4.4 CUDA prefix sum一步一步優化

1. Prefix Sum

前綴求和由一個二元操作符和一個輸入向量組成，雖然名字叫求和，但操作符不一定是加法。先解釋一下，以加法為例：

第一行是輸入，第二行是對應的輸出。可以看到，Output[1] = Input[0] + Input[1]，而Output[length - 1]就是整個輸入向量元素之和。

為什么要使用並行計算？假如用串行計算來計算輸出向量，那么向量中越靠后的元素需要等待的時間越長。最后一個元素需要等待0 + 1 + 2 + ... + (n - 2) = (n - 1)(n - 2) / 2次加法計算后，才開始計算該元素。在很多情況下是不容許等待這么長時間的。因此需要並行計算。

一個很簡單的想法是為輸出向量每一個元素分配一個線程，該線程計算對應元素的值。這個方法有兩個明顯的問題。首先，這樣產生的線程的計算負載非常不均衡，T1只運行一次加法，而T[length - 1]需要運行n - 1次加法。整個並行算法的時間由負載最高的線程決定，也就是O(N)。其次，輸入向量中，第一個元素會被所有線程讀取，第二個被n - 1個線程讀取，總共下來，會有O(N * N)次讀取。前面也提到過，內存讀取是一種比較慢的操作，過多的內存讀取會影響算法的性能。

上面的想法也說明一個現象：並行計算非常的簡單，只要不考慮計算的性能。

 

Exclusive 和 Inclusive區別如下圖

Exclusive 和　Inclusive是可以互相轉換的,　Exclusive把Inclusive添加0作為首元素，去除最后一個元素得到的新的數組就是Exclusive.

 

2. Prefix Sum並行算法一

下圖是一種改進后的算法：

從隊首開始，首先與相鄰元素相加，然后與距離2，距離4的元素相加，以此類推。可以證明其正確性。顯而易見，此算法需要O(logN)步，和O(NlogN)次的加法操作。

以Exclusive為例構造kernel:

__global__ void work_inefficient_scan_kernel(float *X, float *Y, int InputSize) {
　　_shared__ float XY[BLOCK_SIZE];
　　int i = blockIdx.x*blockDim.x + threadIdx.x;
　　if (i < InputSize) {XY[threadIdx.x] = X[i];}
  // the code below performs iterative scan on XY　　
　　for (unsigned int stride = 1; stride <= threadIdx.x; stride *= 2) {
　　　　__syncthreads();　　　　
    　 float in1 = XY[threadIdx.x-stride];
       __syncthreads();
       XY[threadIdx.x] += in1;
  }
}

這個kernel的 block_size = InputSize = sharedMemory size, 通過一個threadBlock完成scan,這個情況搜限制於blocksize的大小，一般是1024，所以在數據量不大的時候（即logN不大），這個算法比較快，可以考慮使用。

 

3. Prefix Sum並行算法二

4.2 CUDA Reduction 一步一步優化里面介紹的思路可以優化Prefix Sum算法.可以分成兩個步驟來執行：

1. Parallel Scan - Reduction Phase

reduction phrase kernel 代碼:

　　for (unsigned int stride = 1; stride < BLOCK_SIZE; stride *= 2) {
　　  　__syncthreads();　
　　  int index = (threadIdx.x+1)*stride*2 - 1; 
    if(index < 2*BLOCK_SIZE)
        XY[index] += XY[index - stride];//index is alway bigger than stride
    __syncthreads();

  }

假設BLOCK_SIZE＝8, 2*BLOCK_SIZE = 16

stride =1,　　threadIdx.x+1 = 1,2,3,4...8, stridek index = 1,3,5,7...15,  第一層thread0 ...thread7都執行計算

stride =2,　　threadIdx.x+1 = 1,2,3,4...8, stridek index = 3,7,11.....31, 第二層thread0 ...thread3都執行計算

stride =3,　　threadIdx.x+1 = 1,2,3,4...8, stridek index = 5,11,19...47, 第三層thread0 ,thread1都執行計算

stride =4,　　threadIdx.x+1 = 1,2,3,4...8, stridek index = 7,23,11.....31, 第四層thread0 都執行計算.

 

我們用Reduction方法得到了一個輸出向量，其元素為圖中黃色的元素。這個輸出向量只包含了部分的最終結果，這8個元素中元素0,1,3, 7計算結束，因此還需要對其進行處理。未得到最終結果的原因是因為有部分的元素對始終沒有得到相加的機會，因此接下來需要將缺失的加法計算補充進去。我們這里使用一個逆Reduction方法，也就是第二個步驟:

2.Reduction Reverse Phase

reduction reverse phase kernel代碼:

 for (unsigned int stride = BLOCK_SIZE/2; stride > 0 ; stride /= 2) {
    __syncthreads();　
　　  int index = (threadIdx.x+1)*stride*2 - 1;
    if((index+stride)< 2*BLOCK_SIZE)
        XY[index + stride] += XY[index];

  }

假設BLOCK_SIZE＝8, 2*BLOCK_SIZE = 16

stride =4,　　threadIdx.x+1 = 1,2,3,4...8, stridek index = 7,23,11.....31,  第一層thread0執行計算,計算element5

stride =2,　　threadIdx.x+1 = 1,2,3,4...8, stridek index = 3,7,11.....31, 第二層thread0 ...thread2都執行計算

stride =1,　　threadIdx.x+1 = 1,2,3,4...8, stridek index = 1,3,5,7...15,, 第三層thread0 ,thread4都執行計算

這一步驟可以計算出:

把這兩個步驟合並到一起來看：

 

注意上圖中已經是16元素的向量了，而不是8元素。從上圖可以看到，逆Reduction第一步補充了Reduction中倒數第二步，依此類推，逆Reduction最后一步補充了Reduction中第一步，至此，所有的元素都得到了相加。

容易看出，上述算法總共需要2 * O(logN)步，以及O(N)次加法操作。相比之前的並行算法，得到了很大的優化。

kernel代碼

__global__ void work_efficient_scan_kernel(float *X, float *Y, int InputSize) {
// XY[2*BLOCK_SIZE] is in shared memory
  _shared__ float XY[BLOCK_SIZE * 2];
　　int i = blockIdx.x*blockDim.x + threadIdx.x;
　　if (i < InputSize) {XY[threadIdx.x] = X[i];}
　　
  // the code below performs iterative scan on XY　　
　　for (unsigned int stride = 1; stride <= BLOCK_SIZE; stride *= 2) {
　　  　__syncthreads();　
　　  int index = (threadIdx.x+1)*stride*2 - 1; 
    if(index < 2*BLOCK_SIZE)
        XY[index] += XY[index - stride];//index is alway bigger than stride
    __syncthreads();
  }
  // threadIdx.x+1 = 1,2,3,4....
  // stridek index = 1,3,5,7...
  
  
  for (unsigned int stride = BLOCK_SIZE/2; stride > 0 ; stride /= 2) {
    __syncthreads();　
　　  int index = (threadIdx.x+1)*stride*2 - 1;
    if(index < 2*BLOCK_SIZE)
        XY[index + stride] += XY[index];
    
  }
  __syncthreads();
  if (i < InputSize) Y[i] = XY[threadIdx.x];
}