【并行計算-CUDA開發】CUDA ---- Warp解析

Warp

邏輯上，所有thread是并行的，但是，從硬體的角度來說，實際上并不是所有的thread能夠在同一時刻執行，接下來我們将解釋有關warp的一些本質。

Warps and Thread Blocks

warp是SM的基本執行單元。一個warp包含32個并行thread，這32個thread執行于SMIT模式。也就是說所有thread執行同一條指令，并且每個thread會使用各自的data執行該指令。

block可以是一維二維或者三維的，但是，從硬體角度看，所有的thread都被組織成一維，每個thread都有個唯一的ID(ID的計算可以在之前的博文檢視)。

每個block的warp數量可以由下面的公式計算獲得：

一個warp中的線程必然在同一個block中，如果block所含線程數目不是warp大小的整數倍，那麼多出的那些thread所在的warp中，會剩餘一些inactive的thread，也就是說，即使湊不夠warp整數倍的thread，硬體也會為warp湊足，隻不過那些thread是inactive狀态，需要注意的是，即使這部分thread是inactive的，也會消耗SM資源。

Warp Divergence

控制流語句普遍存在于各種程式設計語言中，GPU支援傳統的，C-style，顯式控制流結構，例如if…else,for,while等等。

CPU有複雜的硬體設計可以很好的做分支預測，即預測應用程式會走哪個path。如果預測正确，那麼CPU隻會有很小的消耗。和CPU對比來說，GPU就沒那麼複雜的分支預測了（CPU和GPU這方面的差異的原因不是我們關心的，了解就好，我們關心的是由這差異引起的問題）。

這樣我們的問題就來了，因為所有同一個warp中的thread必須執行相同的指令，那麼如果這些線程在遇到控制流語句時，如果進入不同的分支，那麼同一時刻除了正在執行的分之外，其餘分支都被阻塞了，十分影響性能。這類問題就是warp divergence。

請注意，warp divergence問題隻會發生在同一個warp中。

下圖展示了warp divergence問題：

為了獲得最好的性能，就需要避免同一個warp存在不同的執行路徑。避免該問題的方法很多，比如這樣一個情形，假設有兩個分支，分支的決定條件是thread的唯一ID的奇偶性：

__global__ void mathKernel1(float *c) {
    int tid = blockIdx.x * blockDim.x + threadIdx.x;
    float a, b;
    a = b = 0.0f;
    if (tid % 2 == 0) {
        a = 100.0f;
    } else {
        b = 200.0f;
    }
    c[tid] = a + b;
}

一種方法是，将條件改為以warp大小為步調，然後取奇偶，如下：

__global__ void mathKernel2(void) {
    int tid = blockIdx.x * blockDim.x + threadIdx.x;
    float a, b;
    a = b = 0.0f;
    if ((tid / warpSize) % 2 == 0) {
        a = 100.0f;
    } else {
        b = 200.0f;
    }
    c[tid] = a + b;
}

代碼：

【并行計算-CUDA開發】CUDA ---- Warp解析

int main(int argc, char **argv) {
// set up device
int dev = 0;
cudaDeviceProp deviceProp;
cudaGetDeviceProperties(&deviceProp, dev);
printf("%s using Device %d: %s\n", argv[0],dev, deviceProp.name);
// set up data size
int size = 64;
int blocksize = 64;
if(argc > 1) blocksize = atoi(argv[1]);
if(argc > 2) size = atoi(argv[2]);
printf("Data size %d ", size);
// set up execution configuration
dim3 block (blocksize,1);
dim3 grid ((size+block.x-1)/block.x,1);
printf("Execution Configure (block %d grid %d)\n",block.x, grid.x);
// allocate gpu memory
float *d_C;
size_t nBytes = size * sizeof(float);
cudaMalloc((float**)&d_C, nBytes);
// run a warmup kernel to remove overhead
size_t iStart,iElaps;
cudaDeviceSynchronize();
iStart = seconds();
warmingup<<<grid, block>>> (d_C);
cudaDeviceSynchronize();
iElaps = seconds() - iStart;
printf("warmup <<< %4d %4d >>> elapsed %d sec \n",grid.x,block.x, iElaps );
// run kernel 1
iStart = seconds();
mathKernel1<<<grid, block>>>(d_C);
cudaDeviceSynchronize();
iElaps = seconds() - iStart;
printf("mathKernel1 <<< %4d %4d >>> elapsed %d sec \n",grid.x,block.x,iElaps );
// run kernel 3
iStart = seconds();
mathKernel2<<<grid, block>>>(d_C);
cudaDeviceSynchronize();
iElaps = seconds () - iStart;
printf("mathKernel2 <<< %4d %4d >>> elapsed %d sec \n",grid.x,block.x,iElaps );
// run kernel 3
iStart = seconds ();
mathKernel3<<<grid, block>>>(d_C);
cudaDeviceSynchronize();
iElaps = seconds () - iStart;
printf("mathKernel3 <<< %4d %4d >>> elapsed %d sec \n",grid.x,block.x,iElaps);
// run kernel 4
iStart = seconds ();
mathKernel4<<<grid, block>>>(d_C);
cudaDeviceSynchronize();
iElaps = seconds () - iStart;
printf("mathKernel4 <<< %4d %4d >>> elapsed %d sec \n",grid.x,block.x,iElaps);
// free gpu memory and reset divece
cudaFree(d_C);
cudaDeviceReset();
return EXIT_SUCCESS;
}

編譯運作：

$ nvcc -O3 -arch=sm_20 simpleDivergence.cu -o simpleDivergence
$./simpleDivergence

輸出：

$ ./simpleDivergence using Device 0: Tesla M2070
Data size 64 Execution Configuration (block 64 grid 1)
Warmingup elapsed 0.000040 sec
mathKernel1 elapsed 0.000016 sec
mathKernel2 elapsed 0.000014 sec

我們也可以直接使用nvprof（之後會詳細介紹）這個工具來度量性能：

$ nvprof --metrics branch_efficiency ./simpleDivergence

輸出為：

Kernel: mathKernel1(void)
1 branch_efficiency Branch Efficiency 100.00% 100.00% 100.00%
Kernel: mathKernel2(void)
1 branch_efficiency Branch Efficiency 100.00% 100.00% 100.00%

Branch Efficiency的定義如下：

到這裡你應該在奇怪為什麼二者表現相同呢，實際上當我們的代碼很簡單，可以被預測時，CUDA的編譯器會自動幫助優化我們的代碼。稍微提一下GPU分支預測（了解的有點暈，不過了解下就好），這裡，一個被稱為預測變量的東西會被設定成1或者0，所有分支都會得到執行，但是隻有預測值為1時，才會得到執行。當條件狀态少于某一個門檻值時，編譯器會将一個分支指令替換為預測指令，是以，現在回到自動優化問題，一份較長的代碼就會導緻warp divergence了。

可以使用下面的指令強制編譯器不優化（貌似不怎麼管用）：

$ nvcc -g -G -arch=sm_20 simpleDivergence.cu -o simpleDivergence

Resource Partitioning

一個warp的context包括以下三部分：

Program counter
Register
Shared memory

再次重申，在同一個執行context中切換是沒有消耗的，因為在整個warp的生命期内，SM處理的每個warp的執行context都是on-chip的。

每個SM有一個32位register集合放在register file中，還有固定數量的shared memory，這些資源都被thread瓜分了，由于資源是有限的，是以，如果thread比較多，那麼每個thread占用資源就叫少，thread較少，占用資源就較多，這需要根據自己的要求作出一個平衡。

資源限制了駐留在SM中blcok的數量，不同的device，register和shared memory的數量也不同，就像之前介紹的Fermi和Kepler的差别。如果沒有足夠的資源，kernel的啟動就會失敗。

當一個block或得到足夠的資源時，就成為active block。block中的warp就稱為active warp。active warp又可以被分為下面三類：

Selected warp
Stalled warp
Eligible warp

SM中warp排程器每個cycle會挑選active warp送去執行，一個被選中的warp稱為selected warp，沒被選中，但是已經做好準備被執行的稱為Eligible warp，沒準備好要執行的稱為Stalled warp。warp适合執行需要滿足下面兩個條件：

32個CUDA core有空
所有目前指令的參數都準備就緒

例如，Kepler任何時刻的active warp數目必須少于或等于64個（GPU架構篇有介紹）。selected warp數目必須小于或等于4個（因為scheduler有4個？不确定，至于4個是不是太少則不用擔心，kernel啟動前，會有一個warmup操作，可以使用cudaFree()來實作）。如果一個warp阻塞了，排程器會挑選一個Eligible warp準備去執行。

CUDA程式設計中應該重視對計算資源的配置設定：這些資源限制了active warp的數量。是以，我們必須掌握硬體的一些限制，為了最大化GPU使用率，我們必須最大化active warp的數目。

Latency Hiding

指令從開始到結束消耗的clock cycle稱為指令的latency。當每個cycle都有eligible warp被排程時，計算資源就會得到充分利用，基于此，我們就可以将每個指令的latency隐藏于issue其它warp的指令的過程中。

和CPU程式設計相比，latency hiding對GPU非常重要。CPU cores被設計成可以最小化一到兩個thread的latency，但是GPU的thread數目可不是一個兩個那麼簡單。

當涉及到指令latency時，指令可以被區分為下面兩種：

Arithmetic instruction
Memory instruction

顧名思義，Arithmetic instruction latency是一個算數操作的始末間隔。另一個則是指load或store的始末間隔。二者的latency大約為：

10-20 cycle for arithmetic operations
400-800 cycles for global memory accesses

下圖是一個簡單的執行流程，當warp0阻塞時，執行其他的warp，當warp變為eligible時從新執行。

你可能想要知道怎樣評估active warps 的數量來hide latency。Little’s Law可以提供一個合理的估計：

對于Arithmetic operations來說，并行性可以表達為用來hide Arithmetic latency的操作的數目。下表顯示了Fermi和Kepler相關資料，這裡是以（a + b * c）作為操作的例子。不同的算數指令，throughput（吞吐）也是不同的。

這裡的throughput定義為每個SM每個cycle的操作數目。由于每個warp執行同一種指令，是以每個warp對應32個操作。是以，對于Fermi來說，每個SM需要640/32=20個warp來保持計算資源的充分利用。這也就意味着，arithmetic operations的并行性可以表達為操作的數目或者warp的數目。二者的關系也對應了兩種方式來增加并行性：

Instruction-level Parallelism（ILP）：同一個thread中更多的獨立指令
Thread-level Parallelism （TLP）：更多并發的eligible threads

對于Memory operations，并行性可以表達為每個cycle的byte數目。

因為memory throughput總是以GB/Sec為機關，我們需要先作相應的轉化。可以通過下面的指令來檢視device的memory frequency：

$ nvidia-smi -a -q -d CLOCK | fgrep -A 3 "Max Clocks" | fgrep "Memory"

以Fermi為例，其memory frequency可能是1.566GHz，Kepler的是1.6GHz。那麼轉化過程為：

乘上這個92可以得到上圖中的74，這裡的數字是針對整個device的，而不是每個SM。

有了這些資料，我們可以做一些計算了，以Fermi為例，假設每個thread的任務是将一個float（4 bytes）類型的資料從global memory移至SM用來計算，你應該需要大約18500個thread，也就是579個warp來隐藏所有的memory latency。

Fermi有16個SM，是以每個SM需要579/16=36個warp來隐藏memory latency。

Occupancy

當一個warp阻塞了，SM會執行另一個eligible warp。理想情況是，每時每刻到保證cores被占用。Occupancy就是每個SM的active warp占最大warp數目的比例：

我們可以使用的device篇提到的方法來擷取warp最大數目：

cudaError_t cudaGetDeviceProperties(struct cudaDeviceProp *prop, int device);

然後用maxThreadsPerMultiProcessor來擷取具體數值。

grid和block的配置準則：

保證block中thrad數目是32的倍數。
避免block太小：每個blcok最少128或256個thread。
根據kernel需要的資源調整block。
保證block的數目遠大于SM的數目。
多做實驗來挖掘出最好的配置。

Occupancy專注于每個SM中可以并行的thread或者warp的數目。不管怎樣，Occupancy不是唯一的性能名額，Occupancy達到當某個值是，再做優化就可能不在有效果了，還有許多其它的名額需要調節，我們會在之後的博文繼續探讨。

Synchronize

同步是并行程式設計的一個普遍的問題。在CUDA的世界裡，有兩種方式實作同步：

System-level：等待所有host和device的工作完成
Block-level：等待device中block的所有thread執行到某個點

因為CUDA API和host代碼是異步的，cudaDeviceSynchronize可以用來停住CUP等待CUDA中的操作完成：

cudaError_t cudaDeviceSynchronize(void);

因為block中的thread執行順序不定，CUDA提供了一個function來同步block中的thread。

__device__ void __syncthreads(void);

當該函數被調用，block中的每個thread都會等待所有其他thread執行到某個點來實作同步。

__global__ void mathKernel1(float *c) {
    int tid = blockIdx.x * blockDim.x + threadIdx.x;
    float a, b;
    a = b = 0.0f;
    if (tid % 2 == 0) {
        a = 100.0f;
    } else {
        b = 200.0f;
    }
    c[tid] = a + b;
}

__global__ void mathKernel2(void) {
    int tid = blockIdx.x * blockDim.x + threadIdx.x;
    float a, b;
    a = b = 0.0f;
    if ((tid / warpSize) % 2 == 0) {
        a = 100.0f;
    } else {
        b = 200.0f;
    }
    c[tid] = a + b;
}

【并行計算-CUDA開發】CUDA ---- Warp解析

int main(int argc, char **argv) {
// set up device
int dev = 0;
cudaDeviceProp deviceProp;
cudaGetDeviceProperties(&deviceProp, dev);
printf("%s using Device %d: %s\n", argv[0],dev, deviceProp.name);
// set up data size
int size = 64;
int blocksize = 64;
if(argc > 1) blocksize = atoi(argv[1]);
if(argc > 2) size = atoi(argv[2]);
printf("Data size %d ", size);
// set up execution configuration
dim3 block (blocksize,1);
dim3 grid ((size+block.x-1)/block.x,1);
printf("Execution Configure (block %d grid %d)\n",block.x, grid.x);
// allocate gpu memory
float *d_C;
size_t nBytes = size * sizeof(float);
cudaMalloc((float**)&d_C, nBytes);
// run a warmup kernel to remove overhead
size_t iStart,iElaps;
cudaDeviceSynchronize();
iStart = seconds();
warmingup<<<grid, block>>> (d_C);
cudaDeviceSynchronize();
iElaps = seconds() - iStart;
printf("warmup <<< %4d %4d >>> elapsed %d sec \n",grid.x,block.x, iElaps );
// run kernel 1
iStart = seconds();
mathKernel1<<<grid, block>>>(d_C);
cudaDeviceSynchronize();
iElaps = seconds() - iStart;
printf("mathKernel1 <<< %4d %4d >>> elapsed %d sec \n",grid.x,block.x,iElaps );
// run kernel 3
iStart = seconds();
mathKernel2<<<grid, block>>>(d_C);
cudaDeviceSynchronize();
iElaps = seconds () - iStart;
printf("mathKernel2 <<< %4d %4d >>> elapsed %d sec \n",grid.x,block.x,iElaps );
// run kernel 3
iStart = seconds ();
mathKernel3<<<grid, block>>>(d_C);
cudaDeviceSynchronize();
iElaps = seconds () - iStart;
printf("mathKernel3 <<< %4d %4d >>> elapsed %d sec \n",grid.x,block.x,iElaps);
// run kernel 4
iStart = seconds ();
mathKernel4<<<grid, block>>>(d_C);
cudaDeviceSynchronize();
iElaps = seconds () - iStart;
printf("mathKernel4 <<< %4d %4d >>> elapsed %d sec \n",grid.x,block.x,iElaps);
// free gpu memory and reset divece
cudaFree(d_C);
cudaDeviceReset();
return EXIT_SUCCESS;
}

$ nvcc -O3 -arch=sm_20 simpleDivergence.cu -o simpleDivergence
$./simpleDivergence

$ ./simpleDivergence using Device 0: Tesla M2070
Data size 64 Execution Configuration (block 64 grid 1)
Warmingup elapsed 0.000040 sec
mathKernel1 elapsed 0.000016 sec
mathKernel2 elapsed 0.000014 sec

Kernel: mathKernel1(void)
1 branch_efficiency Branch Efficiency 100.00% 100.00% 100.00%
Kernel: mathKernel2(void)
1 branch_efficiency Branch Efficiency 100.00% 100.00% 100.00%

【并行計算-CUDA開發】CUDA ---- Warp解析

Warps and Thread Blocks

Warp Divergence

Resource Partitioning

Latency Hiding

Occupancy

Synchronize

繼續閱讀

Perl與網絡監控

自學Zabbix3.10.2-事件通知Notifications upon events-Actions報警配置點選傳回：自學zabbix集錦

HDU 5678 ztr loves trees

DOG（4）：解析器的部分實作細節先來說說parser一些可能迷惑的地方結果如何傳回?pcd其實是一回事最後的一點說明

拓端tecdat|R語言彈性網絡Elastic Net正則化懲罰回歸模型交叉驗證可視化

二叉樹及其應用--二叉樹建立

PAT (Advanced Level) Practise 1131 Subway Map (30)

ZOJ 3938 Defuse the Bomb

CSU 1565 Word Cloud

ZOJ 3700 Ever Dream

effective c++ 第三版讀書筆記1

可變參數宏， Variadic Macros

ZOJ 1199 Point of Intersection

CSU 1567 Reverse Rot

什麼是BNF範式

詳解STM32單片機的堆棧