Optimizing Parallel Reduction in CUDA

Mark Harris
NVIDIA Developer Technology

并行归约 (Parallel Reduction)

并行归约是一种常见且重要的数据并行基本操作。
在 CUDA 中实现起来很容易，但要做到正确（即高效）则更具挑战性。
它是一个很好的优化示例：
- 我们将逐步介绍 7 个不同版本的实现。
- 展示几种重要的优化策略。

树状归约方法

并行归约在每个线程块（thread block）内部采用树状方法进行计算。

为了处理非常大的数组并保持 GPU 上所有多处理器（multiprocessors）的繁忙状态，需要能够使用多个线程块。每个线程块对数组的一部分进行归约，但这引出了一个问题：如何在线程块之间通信部分归约结果？

问题：全局同步 (Global Synchronization)

如果可以在所有线程块之间进行同步，那么归约大型数组就会变得很简单：

在每个块产生其结果后进行全局同步。
一旦所有块都达到同步点，就递归地继续进行归约。

然而，CUDA 没有全局同步机制。原因如下：
* 对于具有高处理器数量的 GPU，在硬件中构建全局同步机制的成本非常高昂。
* 这会迫使程序员运行更少的线程块（不超过多处理器数量 * 每个多处理器的驻留块数），以避免死锁，但这可能会降低整体效率。

解决方案：将计算分解为多个核函数（kernel）。
* 核函数的启动本身就是一个全局同步点。
* 核函数启动的硬件开销可以忽略不计，软件开销也很低。

解决方案：核函数分解 (Kernel Decomposition)

通过将计算分解为多次核函数调用来避免全局同步。

Page 5 展示了通过两级核函数调用来完成大规模归约的过程。第一级（Level 0）由8个块并行处理，第二级（Level 1）由1个块汇总第一级的结果。

在归约的场景中，所有层级的代码都是相同的，因此可以进行递归的核函数调用。

我们的优化目标是什么？

我们的目标是达到 GPU 的峰值性能。为此，需要选择正确的衡量指标：
* GFLOP/s：适用于计算密集型（compute-bound）的核函数。
* 带宽 (Bandwidth)：适用于访存密集型（memory-bound）的核函数。

归约操作的算术强度（arithmetic intensity）非常低（每加载一个元素仅执行一次浮点运算），这是一种带宽优化的场景。因此，我们应该追求峰值带宽。

本示例将使用 G80 GPU：
* 384-bit 内存接口
* 900 MHz DDR
* 理论峰值带宽 = 384 * 1800 / 8 = 86.4 GB/s

优化历程

归约 #1：交错寻址 (Interleaved Addressing)

这是一个基础的并行归约核函数实现。

代码 (reduce0):

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

工作原理图示：
在归约的每一步中，步长（stride）加倍，只有一部分线程（ID为偶数倍步长的线程）处于活动状态。

Page 8 展示了交错寻址的归约过程。随着步长从1增加到8，活跃线程的数量逐渐减少。

问题分析:
if (tid % (2*s) == 0) 这行代码存在两个主要问题：
1. 高度发散的 warp (highly divergent warps)：在同一个 warp 中，线程会根据其 tid 走不同的执行路径，这非常低效。
2. 取模运算符 (%) 非常慢。

性能表现:
对于 4M（2^22 个整数）元素的归约：

注意: 所有测试的块大小（Block Size）均为 128 个线程。

归约 #2：交错寻址（消除分支发散）

通过替换内层循环中的发散分支来优化。

代码对比:
* 原发散分支:

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

工作原理图示:
这种方法虽然解决了分支发散，但引入了新的问题。

Page 12 展示了消除分支发散后的工作流程，并指出了新的问题：共享内存的 bank 冲突。

新问题: 共享内存 bank 冲突 (Shared Memory Bank Conflicts)。

性能表现:
消除分支发散带来了显著的性能提升。

Page 13 对比了内核1和内核2的性能数据，内核2的速度是内核1的2.33倍。

归约 #3：顺序寻址 (Sequential Addressing)

为了解决 bank 冲突，我们采用顺序寻址。线程 tid 读取 tid + s 的值，其中 s 从 blockDim.x/2 开始递减。这种访问模式是无冲突的。

工作原理图示:

Page 14 展示了顺序寻址的归约过程，这种方式可以避免共享内存的bank冲突。

代码实现:
只需替换内层循环中的跨步索引逻辑。

原跨步索引:

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

使用反向循环和基于 threadID 的索引:

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

下图展示了对4M（2^22）个整数进行归约操作时，前三个内核版本的性能对比。

Kernel 1: 采用交错寻址（interleaved addressing）和发散分支（divergent branching），耗时8.054 ms，带宽为2.083 GB/s。
Kernel 2: 采用交错寻址但存在bank冲突，通过消除发散分支，性能提升。耗时3.456 ms，带宽4.854 GB/s，实现了2.33倍的单步加速。
Kernel 3: 采用顺序寻址（sequential addressing），消除了bank冲突，性能进一步提升。耗时1.722 ms，带宽9.741 GB/s，实现了2.01倍的单步加速，累积加速达到4.68倍。

归约 #4：加载时执行首次相加 (First Add During Load)

空闲线程问题:
在版本3的归约循环过程中，for循环的每次迭代都会使参与计算的线程数量减半（s>>=1）。这意味着在第一次循环迭代中，就有一半的线程处于空闲状态，这是对计算资源的浪费。

优化策略:
为了解决空闲线程问题并进一步优化，我们采取了“加载时执行首次相加”的策略。我们将块（block）的数量减半，并让每个线程加载两个元素，在加载过程中完成第一次归约相加操作，然后将结果写入共享内存。

原始加载方式：每个线程从全局内存加载一个元素到共享内存。

优化后：每个线程加载两个元素，并执行第一次归约相加。

性能更新:
下表更新了包含Kernel 4的性能数据。

Kernel 4: 在全局加载期间执行首次相加。耗时减少到0.965 ms，带宽提升至17.377 GB/s，相较于Kernel 3实现了1.78倍的单步加速，累积加速达到8.34倍。

归约 #5：展开最后一个Warp (Unrolling the Last Warp)

指令瓶颈 (Instruction Bottleneck):
尽管带宽达到了17 GB/s，但这仍然远未达到硬件的带宽上限。我们知道归约操作的算术强度（arithmetic intensity）很低。因此，一个可能的瓶颈是指令开销，这些开销来自于非加载、存储或核心计算的辅助指令，如地址计算和循环开销。策略是展开循环 (unroll loops)。

展开最后一个Warp:
随着归约的进行，“活跃”线程的数量不断减少。当 s <= 32 时，只剩下一个warp的线程在工作。在一个warp内部，指令是SIMD（单指令多数据）同步的。这意味着当 s <= 32 时：
- 我们不再需要 __syncthreads()，因为warp内的同步是隐式的。
- 我们不再需要 if (tid < s) 判断，因为它不会节省任何工作（warp内的所有线程要么都执行，要么都不执行）。
基于以上分析，我们可以展开内部循环的最后6次迭代。

实现:
我们创建一个 warpReduce 设备函数来处理最后32个元素（一个warp）的归约。
重要提示: 为了保证正确性，必须使用 volatile 关键字。这可以防止编译器将共享内存的读写操作优化到寄存器中，确保数据在线程间的可见性。

然后，在主循环之后，当线程ID小于32时，调用 warpReduce 函数。

注意: 这种优化不仅为最后一个warp节省了不必要的工作，而是为所有warps都节省了工作。在未展开的版本中，所有warps在循环的每次迭代中都会执行 for 循环和 if 语句的开销。

性能更新:
下表更新了包含Kernel 5的性能数据。

Kernel 5: 展开最后一个warp。耗时进一步减少到0.536 ms，带宽达到31.289 GB/s，实现了1.8倍的单步加速，累积加速达到15.01倍。

归约 #6：完全展开 (Complete Unrolling)

完全展开:
如果我们在编译时就知道迭代次数，我们就可以完全展开归约循环。幸运的是，GPU限制了块大小（例如最多512个线程），并且我们通常坚持使用2的幂次的块大小。但问题是，如何为编译时未知的块大小进行通用化展开？

解决方案：使用模板 (Templates)！CUDA支持在设备函数和主机函数上使用C++模板参数。我们可以将块大小（block size）指定为函数模板参数。

实现:
通过使用模板，我们可以编写一个完全展开的归约版本。下图中所有用红色标出的代码都将在编译时进行评估。例如，if (blockSize >= 64) 这样的判断，blockSize 是一个模板参数（编译时常量），因此编译器会直接生成特定大小的代码，移除不必要的分支，产生一个非常高效的内层循环。

调用模板化的内核:
我们仍然需要在运行时确定块大小，那么如何在编译时不知道块大小的情况下调用这些模板化的内核呢？我们只需要一个 switch 语句来处理10种可能的块大小。

最终性能:
下表展示了所有六个内核版本的最终性能对比。

Kernel 6: 完全展开版本。耗时缩短至0.381 ms，带宽高达43.996 GB/s，实现了1.41倍的单步加速，累积加速达到惊人的21.16倍。

并行归约的复杂度分析

步骤复杂度 (Step Complexity)：并行归约需要 Log(N) 个并行步骤，每个步骤 s 执行 N/2^s 个独立操作。因此步骤复杂度为 O(log N)。
工作复杂度 (Work Complexity)：对于 N=2^D 的情况，总操作数为 Σ_{s∈[1..D]} 2^(D-s) = N-1。工作复杂度为 O(N)。这是工作高效 (work-efficient) 的，因为它执行的操作数不比串行算法多。
时间复杂度 (Time Complexity)：对于 P 个物理并行处理器，时间复杂度为 O(N/P + log N)。
- 与串行归约的 O(N) 相比有显著提升。
- 在一个线程块内，N=P，所以时间复杂度为 O(log N)。

关于成本 (Cost) 的讨论

并行算法的成本 = 处理器数量 × 时间复杂度
- 如果分配 O(N) 个线程（而不是处理器），时间复杂度为 O(log N)，那么成本为 O(N log N)。这不是成本高效 (cost efficient) 的。
Brent定理 建议使用 O(N/log N) 个线程。
- 每个线程执行 O(log N) 的串行工作。
- 然后 O(N/log N) 个线程协作完成 O(log N) 个步骤。
- 成本 = O(N/log N) * log N = O(N)，这是成本高效的。
这种技术有时被称为算法级联 (algorithm cascading)，在实践中可以带来显著的加速。

算法级联 (Algorithm Cascading)

结合顺序归约和并行归约。
- 每个线程加载多个元素到共享内存并进行求和。
- 在共享内存中进行基于树的归约。
Brent 定理指出，每个线程应该对 O(log n) 个元素求和。
- 例如，每个块处理 1024 或 2048 个元素，而不是 256 个。
根据经验，进一步增加每个线程处理的元素数量是有益的。
- 每个线程承担更多工作可能更好地隐藏延迟。
- 每个块使用更多线程可以减少递归核函数调用的层级。
- 在最后几层，当块数量很少时，核函数启动的开销会很高。
在 G80 架构上，使用 64-256 个块，每块 128 个线程，每个线程处理 1024-4096 个元素时性能最佳。

归约优化 #7: 每线程多次加法

该优化通过让每个线程处理多个元素来进一步提高性能。

将原先加载并相加两个元素的代码：

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

替换为使用 while 循环来根据需要添加任意数量的元素：

unsigned int tid = threadIdx.x;
unsigned int i = blockIdx.x*(blockSize*2) + threadIdx.x;
unsigned int gridSize = blockSize*2*gridDim.x;
sdata[tid] = 0;

while (i < n) {
    sdata[tid] += g_idata[i] + g_idata[i+blockSize];
    i += gridSize;
}
__syncthreads();

保持访存合并

在 while 循环中，步长（stride）被设置为 gridSize，即 blockSize*2*gridDim.x。这是一个关键点，因为它确保了对全局内存的访问是合并的（coalesced），从而保持了高内存带宽。线程块内的连续线程访问连续的内存地址，即使在循环的每次迭代中，这个模式也得以维持。

最终性能总结

400万元素归约的性能

下表总结了七个版本的核函数在处理 4M (2^22) 个整数时的性能演进。从最初的交错寻址与分支分化（Kernel 1）到最终的每线程处理多个元素（Kernel 7），总共实现了 30.04倍 的累积加速。最终版本的带宽达到了 62.671 GB/s。

在处理 32M 元素时，Kernel 7 的带宽更是达到了 73 GB/s。

最终优化核函数代码

下面是最终优化后的 CUDA C++ 核函数代码。它包括一个用于 warp 内部归约的设备函数 warpReduce，以及主核函数 reduce6。reduce6 实现了每线程多次加法的逻辑，并在共享内存中执行并行归约，最后由第一个线程将结果写回全局内存。

性能比较

下图展示了七个不同优化版本的核函数在不同数量元素下的性能表现。Y轴（时间）为对数尺度。可以看出，从版本1到版本7，性能持续提升，最终版本（曲线7）在所有测试的数据规模上都表现出最低的执行时间（即最快的速度）。

优化类型总结

有趣的观察:
算法优化
- 包括寻址方式的改变和算法级联。
- 两者结合带来了 11.84倍 的加速！
代码优化
- 主要是循环展开。
- 带来了 2.54倍 的加速。

这表明高层次的算法级优化比底层的代码级优化能带来更显著的性能提升。

结论

理解 CUDA 性能特征
- 内存合并 (Memory coalescing)
- 分支分化 (Divergent branching)
- 银行冲突 (Bank conflicts)
- 延迟隐藏 (Latency hiding)
使用峰值性能指标来指导优化
理解并行算法的复杂度理论
知道如何识别瓶颈类型
- 例如，是内存、核心计算还是指令开销瓶颈。
先优化算法，然后才展开循环
使用模板参数来生成最优代码

问题联系方式: [email protected]