CUDA, 軟體抽象的幻影背後 之三

標籤：cuda   gpu計算   指令吞吐   線程間通訊   divergent   

本文原載於我的首頁：planckscale.info，轉載於此。

著作權聲明：原創作品，歡迎轉載，但轉載請以超連結形式註明文章來源(planckscale.info)、作者資訊和本聲明，否則將追究法律責任。

上一篇中談到了編程模型中的Block等概念如何映射到硬體上執行，以及CUDA如何用並行來掩蓋延遲。這一篇繼續剖析SIMT，談一談控制流程分叉，指令吞吐和線程間通訊機制。
雖然我們說warp中的線程類似於SIMD,但事實上它是真正的線程。warp中的每一個thread都有自己的指令地址寄存器，允許它們各自執行不同的任務（控制流程分叉）。最簡單的，比如一個

[php]if(threadIdx < 10){...}else{...}[/php]

語句，將threadIdx=0...31這一個warp劃分成兩個分支，各自做不同的事情。這個靈活性以效能為代價，當一個warp中控制流程出現分叉時，不同分支的線程會被分組相繼執行，直到各分支執行完畢後，控制流程重新匯聚成一支（上例中即if語句的結束點）。這種情況下執行單元的利用率較低，因為每個分支執行時都需要關閉其他分支的線程，所以這時一些執行單元是用不到的。
為了儘可能高效的計算，需要約束控制流程分叉的出現。除了減少流程式控制制語句外，還需要注意，並不是只要有流程式控制制語句就一定會帶來控制流程分叉。關鍵是，控制流程分叉只是針對同一warp中的線程而言，不同warp的線程原本就是序列化執行的，分叉對其無影響。因此，只有流程式控制制語句的條件在
同一warp內不一致時，才會有控制流程分叉。這樣,諸如

[php]if(threadIdx.x / WARPSIZE < n){...}else{...}[/php]

這樣的語句是不會有分叉的。當然，更寬鬆的條件如

[php]if(blockIdx.x < n){...}else{...}[/php]

也不會有分叉。依賴於輸入資料的條件如

[php]if(globalArray[threadIdx.x] < n){...}else{...}[/php]

則會帶來分叉。

CUDA的指令都是針對一個warp中32個線程的並行指令，因而一條指令需要在每個線程中都被執行完才算執行完畢。對於簡單的指令如32位浮點數的加、乘，32位整數的加減等，通常都可以由CUDA Core在一個刻度內完成，而每個SM中通常都有不少於32個的CUDA Cores，因而對於一個warp中上述類型的簡單指令，就是一個CUDA Core處理一個線程，一周期內就可以執行完畢。而對於一些較複雜的指令，執行單元並不能提供這麼高的吞吐率，此時一個warp中32個操作需要在多個周期內序列化處理。
我們可以用單位周期內進行的運算元目N除以32來計算指令的吞吐率。以GM204為例，它的SM中有32*4 = 128個CUDA Cores，32個SFU（特殊函數單元），因而在計算32位浮點加法時具有最高吞吐，一個周期內完成128次操作，單位周期內指令吞吐為128/32 = 4；而計算如sin/cos等超越函數時線程不再一一分配到CUDA Cores上，而是要在32個SFU上計算，單位周期內只能完成32次操作，指令吞吐為1條指令每周期.
指令的吞吐率資料可參考CUDA C Programming Guide中 5.4.1. Arithmetic Instructions，該小節以單位刻度每SM上能夠進行的運算元的形式給出了各指令的吞吐率。
指令吞吐率是我們進行效能最佳化的有一個重要指標。通常，影響指令吞吐率的因素除了數值計算操作的複雜度、精確度之外，控制流程分叉也是一個貢獻因子。這裡的原因不難理解，控制流程分叉時執行單元的利用率下降，使得單位周期內執行的運算元目下降，從而降低了指令吞吐。

 

到這裡，硬體圖景下線程的執行就基本說完了，只剩下一個留到最後的話題：線程間互動。通常，不存在任何相互作用的線程，它們之間才能夠以任意的順序執行，像block。但對於warp這樣的線程組，是可能與同一block中其他warp通訊或同步的，這時執行順序就不能任意。所幸即便在block之內，線程間的互動仍然是較弱的，因而底層可以將block劃分成warp來分組序列化執行，遇到互動時再另作處理。我們現在來看看這些互動機制。

線程間互動可以細分為通訊和同步兩類。通訊主要由公用儲存地區交換資料來實現，但也不排除像shuffle這樣的特殊方式存在。
從通訊的粒度來看，可以分為warp內部線程間通訊，block內部線程間通訊，block間通訊，更粗的粒度這裡不考慮。block之間的通訊則只能基於global memory，block內部的通訊主要基於shared memory/global memory，warp內部線程間除了可以利用上述所有方式，還有一種特殊的shuffle機制.下面我們以通訊的粒度分類陳述各種通訊的實現方式。

block間通訊通常基於兩次kernel發射，一次將通訊資料寫入global memory，另一次發射讀global memory進行後續處理。這種通訊開銷較大，主要來自於global memory訪存和kernel發射，所以如果有可能，盡量把任務放在一次kernel發射中完成。
或許有人會問，同一個kernel發射中的兩個block具有共同的global memory，是不是也可以利用這個特點來構造同一kernel下block間的通訊呢？通常的答案是no，因為block之間執行順序不定，很難構造有意義的通訊；但如果要較真，答案是yes，我們真的可以構造一些特殊的block間通訊方式。一個例子如下所示，該執行個體來自於CUDA C Programming Guide B.5. Memory Fence Functions：

[php]__device__ unsigned int count = 0;__shared__ bool isLastBlockDone;__global__ void sum(const float* array, unsigned int N,volatile float* result){// Each block sums a subset of the input array.float partialSum = calculatePartialSum(array, N);if (threadIdx.x == 0) {// Thread 0 of each block stores the partial sum// to global memory. The compiler will use// a store operation that bypasses the L1 cache// since the "result" variable is declared as// volatile. This ensures that the threads of// the last block will read the correct partial// sums computed by all other blocks.result[blockIdx.x] = partialSum;// Thread 0 makes sure that the incrementation// of the "count" variable is only performed after// the partial sum has been written to global memory.__threadfence();// Thread 0 signals that it is done.unsigned int value = atomicInc(& count, gridDim.x);// Thread 0 determines if its block is the last// block to be done.isLastBlockDone = (value == (gridDim.x - 1));}// Synchronize to make sure that each thread reads// the correct value of isLastBlockDone.__syncthreads();if (isLastBlockDone) {// The last block sums the partial sums// stored in result[0 .. gridDim.x-1]float totalSum = calculateTotalSum(result);if (threadIdx.x == 0) {// Thread 0 of last block stores the total sum// to global memory and resets the count// varialble, so that the next kernel call// works properly.result[0] = totalSum;count = 0;}}}[/php]

代碼 1. block間通訊實現數組求和
本代碼摘錄自 CUDA C Programming Guide B.5. Memory Fence Functions

該例實現一個數組的求和，首先各個block計算部分和，然後由最後一個完成部分和計算的block再把所有的部分和加和出最終結果。block間通過一個位於global memory的變數count通訊，它記錄了目前已經完成計算的線程數。這樣，最後一個完成部分和計算的block就會發現count的數值為最大線
程id，因此可以判定需要由它自己來完成最後從部分和向總和的計算。
不過，為了更好的軟體結構，最好還是避免同一kernel的block間產生耦合。同一kernel中block的通訊還涉及到CUDA的weakly-ordered記憶體模型問題，一個線程中先後兩次記憶體操作在另一個線程看來未必能夠保持原有順序，這產生了相當大的複雜性。我們在下文還會提到這一問題。

block內的線程通訊機制較為豐富，尤其是線程同屬一個warp時的shuffle機制。shuffle在Kepler後出現，是一種相當快的線程間通訊方式，它允許同屬一個warp的線程間可以互相引用彼此的寄存器，比如下例：

[php]__global__ void bcast(int arg){int laneId = threadIdx.x & 0x1f;int value;if (laneId == 0) // Note unused variable forvalue = arg; // all threads except lane 0value = __shfl(value, 0); // Get "value" from lane 0if (value != arg)printf("Thread %d failed.\n", threadIdx.x);}[/php]

代碼 2. shuffle機制實現一個值向整個warp的廣播
本代碼摘錄自 CUDA C Programming Guide B.14. Warp Shuffle Functions

laneId是warp中線程的一個index，有threadIdx對32取餘得到。__shfl(value, 0)語句使得各線程能夠訪問laneId==0這一線程中value的值。

更常用的通訊機制自然是shared memory和global memory了。其中shared memory更快速，在大多數時候是構建高效能CUDA程式的必由之路。這些常識不再贅述。基於shared/global memory的線程間資料交換，一定要注意線程的同步。block中線程的同步由__syncthreads()實現。線程會等待同block中其他線程都執行到這一點，並且__syncthreads()語句之前的所有shared/global memory操作都塵埃落定，保證block內所有線程在__syncthreads()之後都能看到這些操作的結果。

最後談一下CUDA採用的weakly-ordered記憶體模型。它導致一個線程中相繼執行的兩個儲存空間操作在另一個線程看來未必是一樣的順序。例如：

[php]__device__ int X = 1, Y = 2;//thread 0__device__ void writeXY(){X = 10;Y = 20;}//thread 1__device__ void readXY(){int B = Y;int A = X;}[/php]

代碼 3. weakly-ordered記憶體模型樣本
本代碼摘錄自 CUDA C Programming Guide B.5. Memory Fence Functions

這段代碼可能產生A=1，B=20這樣的結果。原因是有多種可能的，要麼thread 1看到的X、Y的寫入順序被顛倒，要麼thread 1中讀取順序被顛倒。這種看似相當毀三觀的事情確確實實發生在我們的代碼背後。在一個線程裡兩個相繼但無依賴的記憶體操作，其實際完成的順序可能是不確定的。在這個線程
看來這並沒有導致什麼不同，因為兩個操作無依賴，並不會破壞因果關係鏈結；但在另一個線程的眼裡，它就暴露出來了。
忍不住插句嘴，這簡直就是狹義相對論的世界觀在電腦世界的翻版：一個參考系的觀察者所看到兩個類空間隔事件(可以是相繼發生但因距離遙遠而無因果關聯）在另一個參考系中看來是顛倒的，但有因果關聯的兩事件在所有觀察者看來時序都不會改變。好玩吧？

所以，表面的秩序井然背後有著巨大的複雜性怪獸，為了關牢它的笼子，我們需要約束我們的代碼，用合適的機制來實現線程間通訊。要保證另一個線程看起來，兩組儲存空間操作具有我們所希望的順序，需要用 Memory Fence Function. 這裡不再涉及，對更多細節感興趣的同學，請參考CUDA C Programming Guide B.5. Memory Fence Functions等章節。
（未完待續）

CUDA, 軟體抽象的幻影背後 之三