【CUDA】Reduce归约求和（上）

前言

学习 UP 主 比飞鸟贵重的多_HKL 的 【CUDA】Reduce规约求和（已完结~） 视频，记录下个人学习笔记，仅供自己参考😄

refer1：【CUDA】Reduce规约求和（已完结~）

refer2：Optimizing Parallel Reduction in CUDA

refer3：https://github.com/BBuf/how-to-optim-algorithm-in-cuda/tree/master/reduce

refer4：深入浅出GPU优化系列：reduce优化

refer5：https://chatgpt.com/

1. 问题介绍

通俗来说，reduce 就是要对一个数组求 sum，min，max，avg 等等。reduce 又称为归约，意思是递归约减，因为一般来说最后获得的输出相比于输入维度上会递减。比如 NVIDIA 博客 上的 reduce sum 归约求和问题，一个长度为 8 的数组求和之后得到的输出只有一个数，从一维数组变成了一个标量，如下图所示：

我们以树形图的方式去执行数据累加，最终得到总和，但由于 GPU 没有针对 global memory 的同步操作，只能针对 block 中的 shared memory 数据进行同步，所以，一般而言我们会将 reduce 分为两个阶段，如下图所示：

在第一阶段中，开启 m 个 block 计算出 m 个小份的 reduce 值；在第二阶段中，使用一个 block 将 m 个小份的 reduce 值再次进行 reduce，得到最终的结果。由于第二阶段本质上是可以调用第一个阶段的 kernel，所以不做单独说明，这里我们只探索第一阶段的优化技巧

Note：硬件环境 NVIDIA RTX3060，CUDA-11.6

本文以 reduce sum 为例来记录 reduce 优化，同时学习 Nsight Compute 等工具的简单使用

2. reduce baseline算法介绍

NVIDIA 博客中提供的 reduce sum 的 baseline 算法实现如下：

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

这里的 g_idata 表示的是输入数据的指针，而 g_odata 则表示输出数据的指针，在核函数中首先会把 global memory 中的数据 load 到 shared memory 当中，然后在 shared memory 中对数据进行 reduce sum 操作，最后将 reduce sum 的结果写回 global memory 中

下面是具体的代码分析

共享内存声明：

extern __shared__ int sdata[];

• 动态共享内存声明：该语句声明了一个 int 整型数组 sdata，存储在 shared memory 共享内存中。这里没有指定数组大小，意味着大小将在 kernel 调用时由主机代码通过 kernel 启动配置中的第三个参数动态指定（<<<gridDim, blockDim, shared memory, stream>>>）
• extern 关键字：表示该数组的大小和存储空间不是在编译时确定的，而是外部定义的（由 kernel 调用时的配置参数决定）
• __shared__ 关键字：指明这个数组位于 device 的 shared memory 区域中，这部分内存在同一个线程块（block）内的所有线程是共享的，并且访问速度远快于全局内存（global memory）

索引计算：

unsigned int tid = threadIdx.x;
unsigned int i = blockIdx.x * blockDim.x + threadIdx.x;

• tid：计算的是每个线程在所属块（block）内部的局部索引，其值从 0 到 blockDim.x - 1
• i：计算的是每个线程在整个线程网格（grid）的全局索引，其值从 0 到 总线程数 - 1

关于局部索引 tid 和全局索引 i 的区别可以参考下图：

Note：图片来自于：https://zhuanlan.zhihu.com/p/645330027

以上图为例，我们设置的 gridDim = Dim3(4096, 1, 1)，blockDim = Dim3(256, 1, 1)，也就是说我们是 1D-layout 的布局，启动的线程块是 4096 个，每个线程块中的线程数量是 256 个，从图中可知红色框线程的索引 tid = 3，i = 515

加载数据到共享内存：

sdata[tid] = g_idata[i];
__syncthreads();

• 根据前面计算的线程索引将全局内存中的数据 g_idata[i] 复制到共享内存 sdata[tid] 中
• 同步：__syncthreads() 确保所有线程都完成了数据加载

归约操作：

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

• 归约思想：每一轮（s 从 1 开始倍增），线程间以二叉树方式将数据对加起来
  • 例如，当 s = 1 时，线程 0、2、4、… 会分别将自己与后面的相邻元素相加
  • 当 s = 2 时，线程 0、4、8、… 再将归约后的结果进行累加
  • 如此反复，直到 s 大于等于 blockDim.x
• 线程条件：if(tid % (2 * s) == 0) 的作用是确保在每一轮中只有特定线程进行加法，避免多个线程同时操作相同的数据
• 同步：每次加法后都调用 __syncthreads()，以确保本轮所有线程都完成加法运算，再进入下一轮，保证共享内存中数据的一致性

关于归约操作我们可以看 NVIDIA 博客中的图来理解：

前面我们说了本文只关注 reduce 第一阶段的优化技巧，也就是开启 m 个 block 计算出 m 个小份的 reduce 值，因此我们来关注其中的某一个 block 是如何做归约的就行，在上图中，假设我们数组大小是 16，同时 block 中启动了 16 个线程：

• s = 1 时，线程 0、2、4、6、8、10、12、14 会分别将自己与后面的相邻线程 1、3、5、7、9、11、13、15 中的元素相加
• s = 2 时，线程 0、4、8、12 会分别将自己与线程 2、6、10、14 中的元素相加
• s = 3 时，线程 0、8 会分别将自己与线程 4、12 中的元素相加
• s = 4 时，线程 0 会将自己与线程 8 中的元素相加

写回全局内存：

if(tid == 0){
    g_odata[blockIdx.x] = sdata[0];
}

补充下完整的代码：

#include <cuda.h>
#include <cuda_runtime.h>
#include <time.h>
#include <stdio.h>

#define N 32*1024*1024
#define BLOCK_SIZE 256

__global__ void reduce_v0(float* g_idata, float* g_odata){
    extern __shared__ float 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];
    }
}

int main(){

    // input host
    float* input_host = (float*)malloc(N * sizeof(float));
    for(int i = 0; i < N; ++i){
        input_host[i] = 2.0;
    }

    // input device
    float* input_device;
    cudaMalloc((void**)&input_device, N * sizeof(float));

    // copy data from host to device
    cudaMemcpy(input_device, input_host, N * sizeof(float), cudaMemcpyHostToDevice);
    
    // output host
    int32_t block_num = (N + BLOCK_SIZE - 1) / BLOCK_SIZE;
    float* output_host = (float*)malloc((N / BLOCK_SIZE) * sizeof(float));

    // output device
    float* output_device;
    cudaMalloc((void**)&output_device, (N / BLOCK_SIZE) * sizeof(float));

    // kernel launch
    dim3 grid(N / BLOCK_SIZE, 1);
    dim3 block(BLOCK_SIZE, 1);
    reduce_v0<<<grid, block, BLOCK_SIZE>>>(input_device, output_device);

    // copy result from device to host
    cudaMemcpy(output_host, output_device, block_num * sizeof(float), cudaMemcpyDeviceToHost);

    // print some result
    for(int i = 0; i < 5; ++i){
        float data = output_host[i];
        printf("output_host[%d] = %.2f\n", i, data);
    }

    return 0;
}

为了测试 kernel 的性能，这里我们设定输入数据长度是 32*1024*1024，也就是 4M 的 float 数据，每个 block 的线程数设定为 256（BLOCK_SIZE），启动的 block 总数为输入数据的总长度除以 256 即 dim3 grid(N / BLOCK_SIZE)，每个 block 都会输出其 reduce sum 的结果，因此输出数据的大小是 (N / BLOCK_SIZE) * sizeof(float)

Note：在之后的 reduce 优化中，我们都将采用上面提到的相同大小的输入数据进行测试

执行指令如下：

mkdir bin
/usr/local/cuda/bin/nvcc -o bin/reduce_v0 reduce_v0_baseline.cu
./bin/reduce_v0

执行完成后其输出如下：

output_host[0] = 512.00
output_host[1] = 512.00
output_host[2] = 512.00
output_host[3] = 512.00
output_host[4] = 512.00

利用 Nsight Compute 去 profile 一下，性能和带宽的测试情况如下：

优化手段	耗时(us)	Memory Throughout(%)	DRAM Throughout(%)	加速比
reduce_baseline	2490us	62.01%	15.44%	~

Note：关于 Nsight Compute 的使用大家可以看下篇文章的拓展小节

OK，下面我们就来仔细看看上面实现的这个 baseline 的 reduce sum 核函数存在着什么问题，有哪些优化的空间呢？🤔

NVIDIA 博客中提到存在着以下两个问题：

• highly divergent warps are very inefficient
• % operator is very slow

一个是 warp divergent，另一个是 % 取余操作非常的慢。warp divergent（线程分歧）是指同一个 warp 内的线程在执行分支指令（如 if/else）时，因条件判断结果不同而走不同执行路径的情况。

我们都知道，在 CUDA 中线程执行时是以 warp 线程束为最小单位去执行的，一个 warp 通常包含 32 个线程，这些线程应同时执行同一指令。当出现分支时，如果部分线程走 if 分支而另一部分线程走 else 分支，就需要序列化这两个分支，这会导致部分线程空闲，从而降低并行效率和整体性能

我们拿之前的图为例，16 个线程在每一轮迭代时实际上都存在大量线程是空闲的，例如 step 1 时线程 1、3、5、7、9、11、13、15 线程都是空闲的，没有最大程度的利用 GPU 硬件

另外关于 % 取余操作的速度问题，在 GPU 架构上，整数除法和取余运算并没有像加法、乘法那样得到高效的硬件加速支持，这些操作往往需要较多的指令周期，因此在性能敏感的 CUDA 核函数中，使用 % 操作会显得比较慢

3. 优化技巧1：交错寻址（Interleaved Addressing）

针对前面提到的 warp divergent 问题，解决方式也比较明了，就是尽可能地让同一个 warp 内的所有线程走到同一个分支里面

NVIDIA 博客中给出的优化技巧 1 的代码如下：

__global__ void reduce_v1(float* g_idata, float* g_odata){
    extern __shared__ float 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){
        int index = 2 * s * tid;
        if(index < blockDim.x){
            sdata[index] += sdata[index + s];
        }
        __syncthreads();
    }

    // write result for this block to global mem
    if(tid == 0){
        g_odata[blockIdx.x] = sdata[0];
    }
}

相比于 baseline 的做法 reduce_v1 主要是将原来的：

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

替换为了：

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();
}

其他部分的代码保持不变，下面我们来具体分析下 NVIDIA 的优化思路

在原始的 reduce_v0 版本中，每个线程都需要判断：

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

这一步骤来决定当前线程是否参与当前步的归约计算，前面我们说这种写法存在两个问题：

• 模运算（% operator）：模运算在 CUDA 中比较慢，因为 GPU 对整数除法和取余的硬件支持有限，会消耗更多指令周期
• warp divergent（线程分歧）：当同一个 warp（通常 32 个线程）内的线程因为条件判断结果不同而分成两组分别执行 if 分支和 else 分支时，会导致分支序列化执行，降低并行效率

为了解决上面提到的两个问题，NVIDIA 提出了交错寻址 interleaved addressing 的方法，它通过下面这段代码实现：

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

NVIDIA 的优化思路包括：

1. 消除模运算

• 原代码中使用 tid % (2 * s) 来判断是否执行加法，而 interleaved addressing 直接通过乘法计算出参与计算的全局索引。乘法运算在 CUDA 上比模运算高效很多，因此这一改动减少了计算开销

2. 规整线程工作分布，减少分支分歧

• 在原始实现中，每个线程都独立判断条件 tid % (2 * s) == 0，不同线程得到不同的判断结果，可能导致同一 warp 内有些线程执行加法而有些不执行，进而引起分支发散
• 优化后，每个线程计算出一个 index = 2 * s * tid，这样，参与计算的线程自然地就是那些满足 index < blockDim.x 的线程，而这一条件下，大部分（甚至整个 warp）的线程在同一轮中都执行相同的操作，从而避免了分支的发散问题

3. 计算模式的变化

• 对于每个归约步（每轮 s 值翻倍），优化后的代码让线程编号 tid 直接映射到它们需要操作的元素位置 index，而不是依赖于线程编号本身进行模运算
• 例如，当 s = 1 时，index = 2 * tid，这使得 tid 为 0、1、2、… 的线程分别操作位置 0、2、4、…；当 s = 2 时，index = 4 * tid，这时线程 0、1、2、… 分别操作位置 0、4、8、…

总的来说，interleaved addressing 的优化并不是改变归约计算的基本逻辑（依然是逐步将相邻元素累加），而是通过重新映射线程与数据的关系，消除了原先判断条件中的模运算和由此引起的 warp 分歧，从而提升了 CUDA 核函数的执行效率。

关于 interleaved addressing 我们可以结合 NVIDIA 博客中的图来理解：

我们可以对比着前面归约操作的那张图来看，我们只看 step 1 就能体会到 interleaved addressing 优化的巧妙了，在前面 reduce_v0 的 step 1 中，线程 0、2、4、6、8、… 偶数线程会执行加法操作，而奇数线程 1、3、5、7、9、… 会空闲，而在 reduce_v1 的 step 1 中，线程 0、1、2、…、127 即一个 block 中的 256 个线程的前半部分会执行加法操作，而后半部分的线程 128、129、130、…、255 会空闲

前面我们有提到 kernel 执行时的最小单位是 warp，而 warp 一般是 32 个线程，现在我们假设 0~31 号线程组成第一个 warp，32~63 号线程组成第二个 warp，…，依此类推，我们会发现一个有趣的现象，在 reduce_v0 的 step 1 中由于偶数线程会执行加法操作，而奇数线程会空闲，这导致一个 block 中所有的 warp 内的线程总是有一半在执行有一半在空闲，也就是所有 warp 都在工作，但 warp 内的所有线程中有一半是空闲的

相反 reduce_v1 的 step 2 中是前半部分线程执行加法操作，后半部分线程空闲，意味着 warp0 到 warp7 是工作的，且每个 warp 中的所有线程都在执行加法，而 warp8 到 warp15 是空闲的，且每个 warp 中的所有线程都在空闲

我们可以让线程空闲，但是最好让一个 warp 内的线程同时空闲或者同时执行，因为我们反复在说 kernel 执行时的最小单位是 warp，显然 reduce_v1 的操作更合理

Note：interleaved addressing 这个优化手段也在告诉我们在编写 CUDA Kernel 时应尽量避免使用 % 取模操作，一个是因为昂贵的运算成本，另一个是因为它可能会引起 warp divergent

性能和带宽的测试情况如下：

优化手段	耗时(us)	Memory throughout(%)	DRAM throughout(%)	加速比
reduce_baseline	2490us	62.01%	15.44%	~
reduce_v1_interleaved_addressing	1800us	85.78%	21.35%	1.38

可以看到这个优化还是非常有效的，相比于 baseline 的性能有 1.38 倍的提升

对于 reduce_v1 来说，它其实带来了一个新的问题—bank conflict，要解决这个问题我们必须先搞清楚什么是 bank conflict 以及为什么会产生 bank conflict

关于 bank conflict 我们在韩君老师的课程中有讲过，这里我们再简单过下，大家感兴趣的可以看看：二. CUDA编程入门-共享内存以及Bank Conflict

我们知道在 CUDA 编程中 32 个 thread 组成一个 warp，一般程序在执行的时候是以 warp 为单位去执行的，也就是说每 32 个 thread 一起执行统一指令，比如同时读/取数据。而 NVIDIA 硬件设计者为了让我们能够比较高效的访问 shared memory 把它也分成了 32 个不同的部分，我们称之为 bank，分别对应 warp 中的 32 个线程，之后让每一个线程去访问它们其中的一个部分，如下图所示：

我们假设一个 block 中包含 256 个线程，这 256 个线程访问 shared memory 的 32 个 bank 时的示意图如下所示：

一个理想的情况就是 warp 中的 32 个 thread 分别访问 shared memory 中的 32 个不同的 bank，没有 bank conflict，而 bank conflict 指的就是在同一个 warp 内，有 2 个或者以上的线程访问了同一个 bank 上不同地址的内存，例如线程 0 访问到了 bank1 上线程 33 位置的内存，此时线程 0 和线程 1 就发生了 bank conflict

那为什么会发生 bank conflict 呢？我们回到 reduce_v1 的 kernel 中，我们以 warp0 为例

在第一次迭代中（s = 1），warp0 的 0 号线程需要去加载 shared memory 的 0 号和 1 号地址，然后写回 0 号地址；而 warp0 的 16 号线程需要加载 shared memory 的 32 和 33 号地址并写回 32 号地址，此时，warp0 的 0 号线程和 16 号线程同时访问了 shared memory 中 bank0 的不同内存地址，发生了 2 路的 bank conflict

类似地，在第二次迭代中（s = 2），warp0 的 0 号线程需要去加载 shared memory 的 0 号和 2 号地址，然后写回 0 号地址；而 warp0 的 8 号线程需要加载 shared memory 的 32 和 34 号地址并写回 32 号地址；同理 16 号线程会加载 64 和 68 号地址，24 号线程会加载 96 和 100 号线程；此时，由于 0、32、64、96 号地址都在同一个 bank 即 bank0 中，而 warp0 的 0、8、16、24 号线程同时访问了 bank0，产生了 4 路的 bank conflict

以此类推，下一次迭代会产生 8 路的 bank conflict，使得整个 kernel 一直受到 bank conflict 的影响

4. 优化技巧2：顺序寻址（Sequential Addressing）

为了避免 bank conflict 我们可以把循环迭代的顺序修改一下

NVIDIA 博客中给出的优化技巧 2 的代码如下：

__global__ void reduce_v2(float* g_idata, float* g_odata){
    extern __shared__ float 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 = 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){
        g_odata[blockIdx.x] = sdata[0];
    }
}

相比于 reduce_v1 的做法，reduce_v2 主要是将原来的：

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();
}

替换为了：

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

其他部分的代码保持不变，下面我们来具体分析下 NVIDIA 的优化思路

在 reduce_v1 版本中，我们通过：

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

来让部分线程执行累加操作，虽然这种方法消除了 % 运算和 warp divergent 问题，但由于线程访问的是跳跃式（交错式）的地址，会引起共享内存的 bank conflict

为了解决这一问题，NVIDIA 提出了顺序寻址 sequential addressing 方法，其实现如下：

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

NVIDIA 的优化思路包括：

1. 重新组织归约顺序，避免 bank conflict

• 在交错寻址 interleaved addressing 中，每轮归约时线程计算出的索引是 index = 2 * s * tid，这会导致线程访问共享内存的地址交错分布，从而可能让同一个 warp 内的多个线程同时访问同一个 bank 的不同地址，从而导致 bank conflict
• sequential addressing 顺序寻址则直接让线程按照连续地址访问数据，归约操作从 s = blockDim.x / 2 开始，让线程 0 到 s - 1 依次访问 sdata[tid] 和紧邻的 sdata[tid + s]。这种连续的模式大大降低了 bank conflict 的可能性

2. 优化归约迭代方式

• 这里的归约循环采用了从半个 block 大小开始，每次将 s 右移一位（即 s /= 2）的方式
• 每次迭代中，只有前 s 个线程参与加法操作，保证了操作数据是相邻的两个元素，且随着 s 的减少，线程组数目也在按比例减少，形成一个自顶向下的连续归约过程

总的来说，sequential addressing 的优化思路在于调整归约的顺序，让线程按连续顺序访问共享内存，从而避免了 bank conflict 问题

为什么 sequential addressing 顺序寻址就可以避免 bank conflict 呢？我们可以结合 NVIDIA 博客中的图来理解：

我们可以对比着前面交错寻址的那张图来看，以 warp0 为例，在前面 reduce_v1 交错寻址的 step1 中，warp0 的 0 号线程需要去加载 shared memory 的 0 号和 1 号地址，然后写回 0 号地址；而 warp0 的 16 号线程需要加载 shared memory 的 32 和 33 号地址并写回 32 号地址，此时，warp0 的 0 号线程和 16 号线程同时访问了 shared memory 中 bank0 的不同内存地址，发生了 2 路的 bank conflict

而在 reduce_v2 顺序寻址的 step1 中，warp0 的 0 号线程需要去加载 shared memory 的 0 号和 128 号地址，然后写回 0 号地址；而 warp0 的 1 号线程需要去加载 shared memory 的 1 号和 129 号地址，然后写回 1 号地址；依次类推，我们会发现去加载的 0、1、2、…、31 号数据刚好是 shared memory 的一行数据，都处于不同的 bank，不存在 bank conflict，同理 128、129、130、…、159 号数据也是 shared memory 的一行数据，处于不同的 bank

以此类推，接下来的迭代中，warp0 中的线程会去加载 shared memory 中不同 bank 的数据，不会发生 bank conflict

Note：sequential addressing 的优化思路提示我们在编写 CUDA Kernel 时应设计连续、顺序的内存访问模式，避免产生 bank conflict，提高整体执行效率

性能和带宽的测试情况如下：

优化手段	耗时(us)	Memory throughout(%)	DRAM throughout(%)	加速比
reduce_baseline	2490us	62.01%	15.44%	~
reduce_v1_interleaved_addressing	1800us	85.78%	21.35%	1.38
reduce_v2_bank_conflict_free	1730us	89.24%	22.01%	1.44

可以看到相比于 reduce_v1_interleaved_addressing 性能和带宽又提升了一些

那我们继续来思考下还有什么优化的空间呢？🤔

目前 reduce_v2 最大的一个问题就是线程的浪费，可以看到对于每个 block 我们启动了 256 个线程，但是经过我们的分析发现在第一轮迭代时只有一半的线程也就是只有 128 个线程在工作，其余的线程空闲；第二轮迭代只有 64 个线程在工作，其余线程空闲；第三轮迭代中只有 32 个线程工作，其余线程空闲；依此类推，在每一轮迭代中都会有大量的线程空闲，太浪费了

5. 优化技巧3：解决Idle线程

为了避免 idle 线程，NVIDIA 博客中给出的优化技巧 3 的代码如下：

__global__ void reduce_v3(float* g_idata, float* g_odata){
    extern __shared__ float sdata[];

    // 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] = g_idata[i] + g_idata[i + blockDim.x];
    __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){
        g_odata[blockIdx.x] = sdata[0];
    }
}

相比于 reduce_v2 的做法，reduce_v3 主要是将原来的：

// 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();

替换为了：

// 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] = g_idata[i] + g_idata[i + blockDim.x];
__syncthreads();

其余部分的代码保持不变，下面我们来具体分析下 NVIDIA 的优化思路

reduce_v3 的优化思路就是我们让每一轮迭代的空闲线程也强行做一点工作，除了从 global memory 中取数之外再额外做一次加法。但需要注意的是，reduce_v3 中我们需要把 block 的数量调成之前的一半，因为这个 kernel 中的每个 block 不再是只管 global memory 中的 256 个元素了，而是管 512 个元素。因此需要修改下 gridDim 的定义，如下所示：

// dim3 grid(N / BLOCK_SIZE, 1);
dim3 grid(N / BLOCK_SIZE / 2, 1);

执行后输出如下：

output_host[0] = 1024.00
output_host[1] = 1024.00
output_host[2] = 1024.00
output_host[3] = 1024.00
output_host[4] = 1024.00

可以看到相比于之前的 kernel，这个 kernel 执行后每个 block 的 reduce sum 结果是 1024，这是因为我们在 kernel 的内部多执行了一次加法操作，每个 block 处理 512 个元素，所以累加起来就是 1024

博主绘制了一个草图来帮助理解：

上图对比了 reduce_v2 和 reduce_v3 执行时不同 block 中 sdata 的计算，从图中我们可以清晰的看到在 reduce_v2 中每个 block 只处理 global memory 中的 256 个元素，而相反 reduce_v3 处理 global memory 中的 512 个元素，reduce_v3 中启动的 block 数量是 65536，恰好是 reduce_v2 中 block 数量的一半

性能和带宽的测试情况如下：

优化手段	耗时(us)	Memory throughout(%)	DRAM throughout(%)	加速比
reduce_baseline	2490us	62.01%	15.44%	~
reduce_v1_interleaved_addressing	1800us	85.78%	21.35%	1.38
reduce_v2_bank_conflict_free	1730us	89.24%	22.01%	1.44
reduce_v3_idle_threads_free	896.19us	89.48%	42.89%	2.78

对于 reduce_v3 来说性能已经算是比较好的了，但是依旧没有达到我们想要的效果，我们再来仔细地看看还有什么可以改进的地方。我们发现，在 reduce_v3 中当进行到最后几轮迭代时，此时 block 中只有 warp0 在工作，而其他 warp 空闲，但此时对于空闲的 warp 还在进行 __syncthreads() 线程同步操作，这造成了极大的浪费

结语

这篇文章我们主要是从三个方面出发优化 reduce，首先 baseline 存在 warp divergent 的问题，我们可以通过交错寻址的方式来避免，但是随之带来了一个新的问题—bank conflict，于是我们将交错寻址方式修改为顺序寻址方式以解决 bank conflict 的问题。之后经过我们的分析发现在 kernel 执行过程中存在着许多空闲的线程，于是我们让空闲的线程额外做了一次加法操作

OK，这就是这篇文章中关于 reduce 优化的全部内容了

篇幅原因，后续的 reduce 优化我们放在下篇文章来讲解，敬请期待😄

下载链接

• Nsight Systems软件包下载链接【提取码:1234】
• Reduce归约求和代码下载链接【提取码:1234】

参考

• 【CUDA】Reduce规约求和（已完结~）
• Optimizing Parallel Reduction in CUDA
• https://github.com/BBuf/how-to-optim-algorithm-in-cuda/tree/master/reduce
• 深入浅出GPU优化系列：reduce优化
• https://chatgpt.com/