4.2 CUDA Reduction 一步一步優化

Reduction並行分析:

每個線程是基於一個樹狀的訪問模型，從上至下，上一層讀取數據相加得到下一層的數據．不停的迭代，直到訪問完所有的數據．

利用這么多的線程塊(thread block)我們需要做的事情如下：

１．　處理非常大的數組

２．　讓GPU的每個處理器保持忙碌

３．　每個thread block迭代減少數組的區域．　比如這個圖，第一次是８個數據，第二次是４個．

但是碰到一個問題，在thread block中的線程可以利用同步，但是每個thread block都處理完了，CUDA中並不能提供block級別的同步機制．為什么CUDA不支持全局同步呢?由兩個原因：

1.　　打造高性能GPU處理器的硬件個數是非常昂貴的，處理器越多越貴．

2.　　這就強制程序員盡可能少的使用block個數以避免產生死鎖，（此處還為弄明白：block個數不能大於處理器個數* ）

這個問題該怎么處理呢，全局同步問題？

利用多個kernel來解決這個問題：

cuda kernel lanuch可以當做全局同步點．

cuda kernel lanuch硬件方面的消耗幾乎可以忽略，軟件消耗非常底.

 Level0 是第一個kernel，level1 是第二個kernel．

 

我們的優化目標是？

1.　　努力達到GPU性能極限．

2.　　選擇合適的度量，有兩種：

　　　GFLOP/s:  (FLOPS是Floating-point Operations Per Second每秒所執行的浮點運算次數的英文縮寫)用於分析計算kernel的計算性能．

　　　Bandwidth:用於分析kernel的內存使用情況.

3. reduction是算數密集度非常低的,每個元素一個FLOP.所以我們需要優化極限帶寬來提高信能.

4.以Nvida G80型號的GPU為例：

　．384bit　存儲接口寬度，900MHZ DDR.  384*1800(DDR 是doubel rate)/8 = 86.4GB/s

 

Reduction1: Interleaved Addressing

kernel代碼：

__global__ void reduce0(int *g_idata, int *g_odata) {
extern __shared__ int sdata[];
    // each thread loads one element from global to shared mem
    unsigned int tid = threadIdx.x;
    unsigned int i = blockIdx.x*blockDim.x + threadIdx.x;
    sdata[tid] = g_idata[i];
    __syncthreads();
    // do reduction in shared mem
    for(unsigned int s=1; s < blockDim.x; s *= 2) {
        if (tid % (2*s) == 0) {
            sdata[tid] += sdata[tid + s];
        }
    __syncthreads();
    }
    // write result for this block to global mem
    if (tid == 0) g_odata[blockIdx.x] = sdata[0];
}    

從圖中可以看出，尋址的方式並不是聯系的，而是交叉的. 所以稱這個kernel為interleaved addressing.

分析kernel代碼:

使用了sharedMemory，這里的大小是16,

unsigned int i = blockIdx.x*blockDim.x + threadIdx.x;
sdata[tid] = g_idata[i];

上面兩行可以看出block的size也是16. 這個kernel是一維的．

看for循環的代碼,s 是訪問內存的步長stride，隨着迭代的深入，每一下一層的stride都會變成2倍,第一次是1，然后一次是2,4,8. 每個block的結果最終加到tid=0的線程里面.這是訪問的存儲的部分．計算部分：

第一次是tid[0]計算0,1兩個元素，tid[2]計算2,3元素．．．．．t[14]計算14,15元素．

所以第二次需要把tid[0],tid[2]．．．tid[14]的結果相加. 得到tid[0],tid[3]．．．tid[11].可以看到stride變化．

但不是每個線程都需要執行計算任務，只有每次只有一半的任務執行，第一次是8,第二次是4個最后是一個，用if(tid %(2*s) ==0)來約束．

在for循環體中有:

if(tid % (2*s) == 0)｛

　sdata[tid] += sdata[tid + s];

   __syncthreads();   

｝

if(tid ==0)g_odata[blockId.x] = sdata[0]. 最終每個block的結果保存在sdata[0]中並賦值給global outputdata. 當然這並沒計算完最終的結果，最終結果需要在host端把這些global的結果累加得到.

 

所以上面的kernel可以看出：

1.　　有很多線程並不執行計算．

2. 　　內存訪問並不連續，而是交叉的．

3.　　if(tid % (2*s) == 0) 會導致大量的warp divergence.降低性能.

性能：

 

Reduction2: 消除warp divergence

把reduction１中的for循環體變成:

for (unsigned int s=1; s < blockDim.x; s *= 2) {
　　int index = 2 * s * tid;
　　if (index < blockDim.x) {
　　　　sdata[index] += sdata[index + s];
　　}
　　__syncthreads();
}

可以看出，改成這樣計算什么都沒變，但是消除了warp divergence,雖然如此，卻引起了新的問題：bank conflict.

如果s =2, 那么thread0,thread2就會有bankconflict.

如果s=4, 那么thread0,thread4就會有bankconflict.關於bank conflict的概念可以參考下面的這篇文章

2.2CUDA-Memory(存儲)和bank-conflict

reduction2的性能:

 

 

Reduction3: sequential addressing(連續尋址,就是合並訪問)

如上圖，如果增加步長的長度,可以起到合並訪問的效果，提高內存訪問速率.

從上圖可以看出stride不在是從小到大變化，而是由大到小的變化,從8到4,2,1

修改for循環體:

for (unsigned int s=blockDim.x/2; s>0; s>>=1) {
　　if (tid < s) {
　　　　sdata[tid] += sdata[tid + s];
　　}
　　__syncthreads();
}

但是bandwidth的性能提升是明顯的,見下圖：

繼續分析，發現有一半的線程是一直沒有在執行計算任務的.這很浪費處理器資源嘛.

 

Reduction4：第一次加載內存時執行相加.

我們把sharedMemory提高到之前的2倍,之前一個block算一個block的和，那么現在一個block計算兩個block大小的和.

for循環題和reduction3的相同,修改加載代碼：

unsigned int tid = threadIdx.x;

unsigned int i = blockIdx.x*(blockDim.x*2)+ threadIdx.x

sdata[tid] = g_idata[i] + g_idata[i+blockDim.x];

這樣sharedMemory,16個數據已經執行了一次加法，把blocksize*2的數據變成了blocksize的數據,后面的循環體和上面的reduction3相同. 所以稱這總叫做:第一次載入執行相加，重點是每個載入的元素已經是一次加法的和，stride步長是blocksize.

reduction4的性能：

 

Reduction5: 修改最后一個warp

這時我們的數據帶寬已經達到了17 GB/s，而我們清楚Reduction的算術強度（arithmetic intensity）很低，因此系統的瓶頸可能是由於Parallel Slowdown，即系統對於指令、調度的花費超過了實際數據處理的花費。在本例中即address arithmetic and loop overhead。

我們的解決辦法是將for循環展開（Unroll the loop）。我們知道，在Reduce的過程中，活動的線程數是越來越少的，當活動的線程數少於32個時，我們將只有一個線程束（Warp）。在單個Warp中，指令的執行遵循SIMD（Single Instruction Multiple Data）模式，也就是說在活動線程數少於32個時，我么不需要進行同步控制，即我們不需要 if (tid < s) 。

修改kernel如下：

首先展開最后一個warp.

for (unsigned int s=blockDim.x/2; s>32; s>>=1)　//步長stride小於32不用進循環，執行下面的動//作，並且不需要同步
{
　　if (tid < s)
　　　　sdata[tid] += sdata[tid + s];
　　__syncthreads();
}

if (tid < 32) {
　　if (blockSize >=64) sdata[tid] += sdata[tid + 32];
　　if (blockSize >=32) sdata[tid] += sdata[tid + 16];
　　if (blockSize >=16) sdata[tid] += sdata[tid + 8];
　　if (blockSize >=8) sdata[tid] += sdata[tid + 4];
　　if (blockSize >=4) sdata[tid] += sdata[tid + 2];
　　if (blockSize >=2) sdata[tid] += sdata[tid + 1];
} // write result for this block to global mem if (tid == 0) g_odata[blockIdx.x] = sdata[0]; }

 

Reduction６:完全展開循環體 

if (blockSize >= 512) {
　　if (tid < 256) { sdata[tid] += sdata[tid + 256]; } __syncthreads();
}
if (blockSize >= 256) {
　　if (tid < 128) { sdata[tid] += sdata[tid + 128]; } __syncthreads();
}
if (blockSize >= 128) {
　　if (tid < 64) { sdata[tid] += sdata[tid + 64]; } __syncthreads();
}
if (tid < 32) {
　　if (blockSize >=64) sdata[tid] += sdata[tid + 32];
　　if (blockSize >=32) sdata[tid] += sdata[tid + 16];
　　if (blockSize >=16) sdata[tid] += sdata[tid + 8];
　　if (blockSize >=8) sdata[tid] += sdata[tid + 4];
　　if (blockSize >=4) sdata[tid] += sdata[tid + 2];
　　if (blockSize >=2) sdata[tid] += sdata[tid + 1];
}

上面代碼所有紅色的部分在編譯階段都會進行優化,結果是一個非常有效率的內循環.
Reduction5,6性能：


到此並行規約算法分析完畢,只要是符合規約類型的運算均符合這種優化思路