CUDA编程-ReduceSum优化记录（文字+图解）

ReduceSum

在GPU中，reduce采用了一种树形的计算方式。如下图所示。

从上至下，将数据不断地累加，直到得出最后的结果，即25。但由于GPU没有针对global数据的同步操作，只能针对block的数据进行同步。所以，一般而言将reduce分为两个阶段，其示意图如下：

我们仔细来看看这个事，假设给定一个长度为N的数组，需要计算该数组的所有元素之和。首先需要将数组分为m个小份。而后，在第一阶段中，开启m个block计算出m个小份的reduce值。最后，在第二阶段中，使用一个block将m个小份再次进行reduce，得到最终的结果。由于第二阶段本质上是可以调用第一个阶段的kernel，所以不做单独说明，本文只是探索第一阶段的优化技巧。

所以kernel接口为：

__global__ void reduce(T *input, T* output)

其中，input代表输入的数组，即一个长度为N的数组，output代表输出数组，即第一阶段的结果，即长度为M的数组。随后要开始激动人心的coding阶段，但在CUDA编程中，我们首先需要设置三个参数:

1. BlockNum：即开启的block数量，即上面所说的M，代表需要将数组切分为几份。

2. Thread_per_block:每个block中开启的线程数，一般而言，取128，256，512，1024这几个参数会比较多。

3. Num_per_block:每个block需要进行reduce操作的长度。

其中，BlockNum* Num_per_block=N，三个参数的示意图如下：

（1）Baseline

Baseline算法比较简单，分为三个步骤。第一个步骤是将数据load至shared memory中，第二个步骤是在shared memory中对数据进行reduce操作，第三个步骤是将最后的结果写回global memory中。代码如下：

__global__ void reduce0(float* d_in, float* d_out) {
    __shared__ float sdata[THREAD_PER_BLOCK];
    // 每个线程从全局内存中读取一个数据到共享内存中
    unsigned int tid = threadIdx.x; // 获取每个thread的x索引，因为thread是一维的，所以y始终为1
    unsigned int i = blockIdx.x * blockDim.x + threadIdx.x; // 通过thread的x索引计算出该数字在原始数组中的索引
    sdata[tid] = d_in[i];
    __syncthreads();
    // 在共享内存中计算reducesum
    for (unsigned int s = 1; s < blockDim.x; s *= 2) {
        if (tid % (2 * s) == 0) {
            sdata[tid] += sdata[tid + s]; // 每次循环将(2k + 1)s加到(2k)s上
        }
        __syncthreads(); // 所有线程每个循环都需要等待同步，然后进入下一次循环
    }
    // 所有线程累加到sdata[0]之后将每个block求和后的值传回全局内存
    if (tid == 0) d_out[blockIdx.x] = sdata[0];
}

第一阶段：我们让Num_per_block与Thread_per_block一致，每个block设定为256个线程，一个block负责256个数据的reduce工作。tid代表线程号，i代表在原始数组中的索引号。第tid号线程将第i号的数据从global中取出，放到shared memory的第tid元素中。比如在第0号block中，0号thread将0号元素取出，放到shared memory的第0号位置。示意图见：

从硬件角度来分析一下代码，为了执行代码，GPU需要分配两种资源：

• 存储资源：包括在global memory中分配的一块32M×sizeof(float)的空间以及在shared memory中分配的256×sizeof(float)的空间。需要注意的是，shared memory存在bank冲突的问题，因而需要格外小心。

• 计算资源：一个block中分配256个thread线程，32个线程为一组，绑定在一个SIMD单元。所以256个线程可以简单地理解为分配了8组SIMD单元。

总而言之，在第一个阶段，就是tid号线程将i号数据从global memory中取出，再放进shared memory中，严谨一点的话，中间是走一遍寄存器再到shared memory中的。

第二阶段：进行多轮迭代，在第一轮迭代中，如果tid % 2 == 0, 则第tid号线程将shared memory中第tid号位置的值和第tid + 1号的值进行相加，而后放在第tid号位置。在第二轮迭代中，如果tid % 4 == 0，则第tid号线程将shared memory中第tid号位置的值和第tid + 2号的值进行相加，而后放在第tid号位置。不断迭代，则所有元素都将被累加到第0号位置。其示意图如下。其中，红色的线程代表符合if条件的线程，只有它们有任务，需要干活。

第三阶段：block负责的256个元素之和都放置在shared memory的0号位置，此时，只需要将0号位置的元素写回即可。

优化结果：

（2）warp分歧

问题描述：

目前reduce0存在的最大问题就是线程束分歧（warp divergent）的问题。对于一个block而言，它所有的thread都是执行同一条指令。如果存在if-else这样的分支情况的话，thread会执行所有的分支。只是不满足条件的分支，所产生的结果不会记录下来。可以在上图中看到，在每一轮迭代中都会产生两个分支，分别是红色和橙色的分支。这严重影响了代码执行的效率。

在 CUDA 编程中，warp divergence 指的是在一个 warp中的不同线程之间执行不同的代码路径（即分支语句的不同分支），导致 warp中的线程不能同时执行相同的指令，从而影响了性能。

当一个 warp 中的线程执行分支语句时，如果不是所有线程都选择了同一个分支路径，就会发生 warp divergence。这是因为 warp 中的线程以 SIMD（Single Instruction Multiple Data）方式执行，每个线程必须执行相同的指令。如果不同的线程选择了不同的分支路径，就会出现某些线程被阻塞等待其他线程完成它们的操作，这就浪费了处理器时间。

优化方法：

虽然代码依旧存在着if语句，但是却与reduce0代码有所不同。我们继续假定block中存在256个thread，即拥有256/32=8个warp

• 第1次迭代：0-3号warp的index<blockDim.x， 4-7号warp的index>=blockDim.x。对于每个warp而言，都只是进入到一个分支内，所以并不会存在warp divergence的情况。

• 第2次迭代：0、1号两个warp进入计算分支。

• 第3次迭代：只有0号warp进入计算分支。

• 第4次迭代：只有0号warp的前16个线程进入分支。此时开始产生warp divergence。

通过这种方式，我们消除了前3次迭代的warp divergence。

优化代码：

这里有个非常需要注意的点，tid代表的是线程的idx，而计算出来的share_idx则代表当前tid线程进行操作的share memory中的数据索引，因此，第一迭代的时候只用到了[0, 127]的线程idx完成了[0, 255]个share memory中数字的累加操作。

// warp divergence
__global__ void reduce1(float* d_in, float* d_out) {
    __shared__ float sdata[THREAD_PER_BLOCK];
    // 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] = d_in[i];
    __syncthreads();
​
    // do reduction in shared mem
    for (unsigned int s = 1; s < blockDim.x; s *= 2) {
        int share_idx = 2 * s * tid;
        if (share_idx < blockDim.x) {
            sdata[share_idx] += sdata[share_idx + s];
        }
        __syncthreads();
    }
​
    // write result for this block to global mem
    if (tid == 0) d_out[blockIdx.x] = sdata[0];
}

优化结果：

（3）bank冲突

问题描述：

reduce1的最大问题是bank冲突，首先介绍一下CUDA中是如何划分bank的：

• bank划分：先假设共有32个bank，这个数字是与一个warp中的thread数对应的，那么计算bank_idx的公式为：bank_idx = tid % 32，其实也就是说share memory中相邻的数据在不同的bank，bank中相邻的数据间隔为32：

• bank冲突：当同一个warp中的多个线程同时申请访问同一个bank中的不同数据时就会出现bank冲突，此时GPU会强制序列化这些线程的访问，从而导致程序的性能降低。

• 为什么这样划分：之所以将相邻的数据划分到不同的bank其实就是为了避免bank冲突的问题，因为一个warp中的32个thread通常是需要处理连续的数据的，因此如果将连续的32个数据放在同一个bank中就肯定会有bank冲突，并且冲突次数非常多，严重影响性能。

• 最理想访问模式：

通过上面的介绍就可以看出reduce1中存在的问题：

第一次迭代：warp0中thread0访问share memory中的share_idx0和share_idx1，thread16访问share memory中的share_idx32和share_idx33，那么share_idx0和share_idx32都属于bank0就发生了2路的bank冲突。

第二次迭代：warp0中thread0访问share memory中的share_idx0和share_idx2，thread8访问share memory中的share_idx32和share_idx34，thread16访问share memory中的share_idx64和share_idx66，thread24访问share memory中的share_idx96和share_idx98，那么share_idx0和share_idx32和share_idx64和share_idx96都属于bank0就发生了4路的bank冲突。

优化方法：

把for循环逆着来。原来stride从0到256，现在stride从128到0。其伪代码如下：

// bank conflict
__global__ void reduce2(float* d_in, float* d_out) {
    __shared__ float sdata[THREAD_PER_BLOCK];
​
    // 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] = d_in[i];
    __syncthreads();
​
    // do reduction in shared mem
    for (unsigned int s = blockDim.x / 2; s > 0; s >>= 1) {
        if (tid < s) {
            sdata[tid] += sdata[tid + s];
        }
        __syncthreads();
    }
​
    // write result for this block to global mem
    if (tid == 0) d_out[blockIdx.x] = sdata[0];
}

专注以一个warp0之中的计算：

第一次迭代：warp0中thread0访问share memory中的share_idx0和share_idx128，thread1访问share memory中的share_idx1和share_idx129，不同thread访问不同的bank因此没有bank冲突。

第二迭代：warp0中thread0访问share memory中的share_idx0和share_idx64，thread1访问share memory中的share_idx1和share_idx65，不同thread访问不同的bank因此没有bank冲突。

第三迭代：warp0中thread0访问share memory中的share_idx0和share_idx32，thread1访问share memory中的share_idx1和share_idx33，不同thread访问不同的bank因此没有bank冲突。

第四迭代：warp0中thread0访问share memory中的share_idx0和share_idx16，thread1访问share memory中的share_idx1和share_idx17，一个thread会访问两个bank，但是不同thread还是访问不同的bank因此没有bank冲突。

优化结果：

（4）idle线程

问题描述：

reduce2最大的问题就是线程的浪费，每个block中启动了256个线程，但实际并没有这么多线程在计算：

第一轮迭代：只有[0, 127]的线程在进行计算，剩余[128, 255]个线程处于空闲状态。

第二轮迭代：只有[0, 63]的线程在进行计算，剩余[64, 255]个线程处于空闲状态。

随着迭代次数增加，参与计算的线程每次都在减半，并且[128, 255]线程除了一开始从global memory读取数据外再也没有参与过计算，一直处于空闲状态。

优化方法：

将block数量减半，num_per_block翻倍，即一个block之前负责求和256个数，现在需要求和512个数。

在每个线程开始阶段从global memory取数据的时候可以直接进行一次从global[256, 511]到share memory[0, 255]的累加，这样的话block中的每个thread都至少进行了一次累加计算。

主要区别在与第一次迭代，所以将修改前和修改后block0中的的第一次迭代进行对比：

从宏观上来看（左侧为2个block每个block负责global memory中256个数字计算，右侧为1个block负责512个数字计算）：

优化代码：

// idle thread
__global__ void reduce3(float* d_in, float* d_out) {
    __shared__ float sdata[THREAD_PER_BLOCK];
​
    // each thread loads one element from global to shared mem
    unsigned int tid = threadIdx.x;
    unsigned int i = blockIdx.x * (blockDim.x * 2) + threadIdx.x;
    sdata[tid] = d_in[i] + d_in[i + blockDim.x]; // 这里相当于从global读取数据的时候先进行累加了一次
    __syncthreads();
​
    // do reduction in shared mem
    for (unsigned int s = blockDim.x / 2; s > 0; s >>= 1) {
        if (tid < s) {
            sdata[tid] += sdata[tid + s];
        }
        __syncthreads();
    }
​
    // write result for this block to global mem
    if (tid == 0) d_out[blockIdx.x] = sdata[0];
}

优化结果：

（5）warp0同步

问题描述：

当进行到最后几轮迭代时，此时的block中只有warp0在干活时，线程还在进行同步操作，这一条语句造成了极大的浪费。

优化方法：

由于一个warp中的32个线程其实是在一个SIMD单元上，这32个线程每次都是执行同一条指令，这天然地保持了同步状态。一个warp内的线程不需要同步；即当执行的线程数小于warpsize时，不需要同步函数。因而当s<=32时，即只有一个SIMD单元在工作时，可以不需要__syncthreads()进行同步。

优化代码：

__device__ void warpReduce(volatile float* cache, unsigned int tid) {
    cache[tid] += cache[tid + 32];
    cache[tid] += cache[tid + 16];
    cache[tid] += cache[tid + 8];
    cache[tid] += cache[tid + 4];
    cache[tid] += cache[tid + 2];
    cache[tid] += cache[tid + 1];
}
​
__global__ void reduce4(float* d_in, float* d_out) {
    __shared__ float sdata[THREAD_PER_BLOCK];
​
    // each thread loads one element from global to shared mem
    unsigned int tid = threadIdx.x;
    unsigned int i = blockIdx.x * (blockDim.x * 2) + threadIdx.x;
    sdata[tid] = d_in[i] + d_in[i + blockDim.x];
    __syncthreads();
    
    //  当调用的线程数小于线程束warp的大小时不需要进行同步，因此终止循环
    for (unsigned int s = blockDim.x / 2; s > 32; s >>= 1) {
        if (tid < s) {
            sdata[tid] += sdata[tid + s];
        }
        __syncthreads();
    }
    
    // 当S<32的时候直接进行[0, 63]、[0, 31]、[0, 15]、[0, 7]、[0, 3]、[0, 1]的6步累加
    if (tid < 32) warpReduce(sdata, tid);
    if (tid == 0) d_out[blockIdx.x] = sdata[0];
}

可以看出此方法下有3个阶段：

• global->share累加阶段：每个线程从global读取数据到share的过程中先进行一次累加：

• 多warp同步阶段：和之前保持一致，避免bank冲突和warp分歧

• warp0直接计算阶段：当s<32的时候直接进行[0, 63]、[0, 31]、[0, 15]、[0, 7]、[0, 3]、[0, 1]的6步累加

优化结果：