CUDA学习（归约算法）

parallel reduction

Nvidia官网归约ppt

中文ppt

可以理解为将数组中所有数求和的过程并行化

CUDA本身并不支持全局同步，将每一层归约作为一个kernel重复递归调用

(1)Interleaved Addressing 交错寻址

数组中相邻地址元素求和

__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] = (i < n) ? g_idata[i] : 0;

  __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];
}

（1）

  // 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] = (i < n) ? g_idata[i] : 0;

将数组元素从global memory拷贝到shared memory，由于shared memory 只在同一个block中共享，

所以shared memory的索引为每个block中threads的Idx，每一个block保存一部分数组元素，而global memory的索引为正常全局索引，

由blockDim和blockIdx及threadIdx计算而得。

同时保证当索引 i >= n 时，将对应shared memory赋值为0

由于数组元素不在同一个warp内，需要利用__syncthreads()进行同步，保证所有threads执行到同一处。

(2）

在shared memory中进行reduction求和

归约分为多步，数组D

第一步

D[0] = D[0] + D[1];
D[2] = D[2] + D[3];
D[4] = D[4] + D[5];
...
D[2*i] += D[2*i+1]  (i = 0, 1, ..., (blockDim.x-1)/2)

只改变thread满足

threadIdx.x%2 == 0

的数组元素

求和元素步长stride=1

第二步

D[0] = D[0] + D[2];
D[4] = D[4] + D[6];
D[8] = D[8] + D[10];
...
D[2*2*i] += D[2*(2*i+1)]  (i = 0, 1, ..., （blockDim.x-2)/4)

只改变thread满足

threadIdx.x % (2*2) == 0

的元素

求和元素步长stride=2

第三步，

D[0] = D[0] + D[4];
D[8] = D[8] + D[12];
D[16] = D[16] + D[20];
...
D[2*2*2*i] += D[2*2*(2*i+1)]  (i = 0, 1, ..., (blockDim.x-4)/8)

只改变thread满足

threadIdx.x % (2*2*2) == 0

的元素

求和元素步长stride=4

第四步，

只改变thread满足

threadIdx.x % (2*2*2*2) == 0

的元素

求和元素步长stride = 8

以此类推

求和元素步长stride：1 2 4 8 ... 2^(i-1) 其中i=0，1，2，...，且2^(i-1)<blockDim.x

满足条件的thread为：2 4 8 16 ... 2*stride

利用同步函数等待所有线程执行结束，

__syncthreads();

注意由于步长stride<blockDim.x，所以元素下标索引不会超过thread的下标范围

(3)

将所有blocks中的shared memory 第一个元素求和