本文深入探讨了一个经典的并行计算算法——并行归约（Parallel Reduction）的性能优化过程，通过七个渐进式的优化步骤，展示了如何将算法性能提升至极致。这项研究基于Mark Harris在NVIDIA网络研讨会中提出的优化方法，在重现这些优化技术的同时，进一步简化了概念阐述以便于理解。配套的GitHub代码库提供了完整的实现细节，为读者深入研究提供了详实的技术支撑。
算法原理分析
并行归约算法是CUDA编程中的一个重要数据并行原语，其核心思想是利用GPU的线程层次结构对向量、矩阵或张量进行并行计算。该算法通过sum()、min()、max()或avg()等操作对数据进行聚合处理。本文将重点使用sum()操作进行数据集归约。尽管这些操作在概念上相对简单，但它们在众多应用场景中发挥着关键作用，因此需要高度优化以避免成为性能瓶颈。
在并行化实现过程中，算法采用基于树的方法，计算任务分布在GPU的各个线程块中。这里面临一个核心技术挑战：如何在线程块之间高效传递部分计算结果？最直观的解决方案是采用全局同步机制——让各个块完成计算后进行全局同步，然后递归继续处理。CUDA架构并不支持全局同步，主要原因是硬件成本过高，且会限制程序员只能使用少量线程块以避免死锁，从而显著降低整体计算效率。
基于树的归约 | 来源：NVIDIA
解决线程块间部分结果通信问题的实用方法是采用内核分解技术。内核分解将大规模的内核任务分解为多个较小的、可管理的子任务，这些子任务可以在不同的线程或块中独立执行。这种方法最大限度地减少了硬件和软件开销，实现了更灵活高效的GPU资源利用，同时降低了同步需求并提升了整体计算性能。
内核分解 | 来源：NVIDIA
性能评估指标体系
算法性能评估基于两个关键维度：执行时间和带宽利用率。这些指标能够准确反映GPU资源利用程度，本质上衡量系统是否达到了峰值性能。我们的优化目标是实现GPU峰值性能，通过计算性能（GFLOP/s）和内存性能（GB/s）两个方面的指标进行量化评估。
为了实现这些指标的优化，需要重点关注数据访问模式和计算瓶颈识别两个核心方面。具体而言，需要评估如何提升数据读写效率，以及如何使计算过程更加快速和高效。在GPU编程中，理想的计算实现不仅要追求高速度，更要确保大部分线程都能有效参与工作。
REDUCE-0：交替寻址基础实现
作为优化的起点，首先实现最基础的并行归约方法。这种朴素的并行化方法需要确定访问存储元素的地址空间模式，检索相应元素，通过求和操作组合这些元素，并在不同线程上递归重复此过程以实现操作的并行化。
交替寻址技术的核心是访问和组合位于当前线程处理段中间位置的地址空间。以包含1024个整数的数组为例，如果每个块使用256个线程，每个线程从不同起点开始，每次跳过256个元素进行处理。线程0将依次处理元素0、256、512和768，每次将当前元素与位于其负责数组段中间位置的另一个元素进行组合。因此，线程0会将元素0与元素128组合，元素256与384组合，元素512与640组合，元素768与896组合。这个过程将递归进行直到获得最终结果。
这种方法在简化线程间同步的同时，确保所有线程都能积极参与并行数据归约，从而实现更加均衡的负载分配和高效的归约计算。
交替寻址 | 来源：NVIDIA
// 归约 0 – 交替寻址 __global__ void reduce0(int *g_in_data, int *g_out_data){ extern __shared__ int sdata[]; // 存储在共享内存中 // 每个线程从全局内存加载一个元素到共享内存 unsigned int tid = threadIdx.x; unsigned int i = blockIdx.x * blockDim.x + threadIdx.x; sdata[tid] = g_in_data[i]; __syncthreads(); // 归约方法 -- 在共享内存中进行，因为sdata存储在那里 for(unsigned int s = 1; s < blockDim.x; s *= 2){ if (tid % (2 * s) == 0) { sdata[tid] += sdata[tid + s]; } __syncthreads(); } if (tid == 0){ g_out_data[blockIdx.x] = sdata[0]; } }实现过程包含六个关键步骤。首先，根据线程ID和块大小为每个线程分配起始索引。接着，每个线程将对应元素从全局内存加载到共享内存。然后同步块内所有线程以确保数据加载完成。在共享内存中执行归约操作时，每个线程将其值与另一个线程在计算偏移量处的值相加，该偏移量在每个后续步骤中减半。每个归约步骤后需要再次同步线程以保证数据完整性。最后，每个块中的第一个线程将归约结果写入全局内存的输出数组。
结果
REDUCE-0 结果
性能瓶颈分析
虽然这种方法为并行编程奠定了良好基础，但仍存在明显的性能瓶颈。从计算和内存两个维度分析，可以识别出以下效率问题。
在计算方面，主要瓶颈来源于模运算符（%）的使用。该运算符在计算上开销很大，因为它涉及除法操作——这是底层硬件上最慢的操作之一。在内核中频繁执行该操作会严重影响性能。此外交替寻址模式导致warp高度发散，因为同一warp内的线程由于条件判断需要执行不同的代码路径。这种路径发散导致warp停滞，等待其他线程完成，严重降低了执行效率。
在内存方面，由于warp发散导致的内存访问模式次优。每个线程访问分布在整个数组中的数据元素，使得内存访问呈现分散且非连续的特点，导致内存带宽利用效率低下和数据检索延迟增高。这种分散的访问模式可能产生多个缓慢的内存事务，而非单个高效的事务，因此无法充分利用GPU的内存带宽能力。
REDUCE-1：改进的交替寻址
针对第一版实现中的计算效率问题，这个版本对寻址方式进行了优化。虽然基本的寻址逻辑保持不变，但在构建归约函数时消除了模运算符和发散条件的使用。通过重构索引计算方式（int index = 2 * s * tid;），REDUCE-1确保每个线程能够一致地执行操作，无需检查相对于步长的位置，从而消除了warp内的发散现象。
这种调整使得warp中的所有线程都遵循相同的执行路径，显著提升了warp执行效率。移除模运算符进一步提升了性能，避免了GPU上运行缓慢的除法相关运算。
// 归约 1 – 无分支发散和%运算的交替寻址 __global__ void reduce1(int *g_in_data, int *g_out_data){ extern __shared__ int sdata[]; // 存储在共享内存中 // 每个线程从全局内存加载一个元素到共享内存 unsigned int tid = threadIdx.x; unsigned int i = blockIdx.x * blockDim.x + threadIdx.x; sdata[tid] = g_in_data[i]; __syncthreads(); // 归约方法 -- 在共享内存中进行 for(unsigned int s = 1; s < blockDim.x; s *= 2){ // 注意步长为 s *= 2：这导致交替寻址 int index = 2 * s * tid; // 现在我们不需要if条件的发散分支 if (index + s < blockDim.x) { sdata[index] += sdata[index + s]; // s用来表示将要合并的偏移量 } __syncthreads(); } if (tid == 0){ g_out_data[blockIdx.x] = sdata[0]; } }结果
REDUCE-1 结果
新问题的出现
尽管REDUCE-1在计算效率和执行一致性方面相比REDUCE-0有显著改进，但引入了一个新的性能问题：共享内存库冲突（Bank Conflicts）。当多个线程同时尝试访问同一内存库的数据时，就会发生这种冲突，导致原本可以并行执行的内存访问被强制串行化。
从REDUCE-0到REDUCE-1的转换过程中，虽然提升了算法的计算效率，但并未解决内存相关问题，反而通过引入步长机制创造了更多内存访问冲突。步长方法使得线程尝试访问相同的共享内存地址。REDUCE-0将线程分散在充当边界的间隔内，将线程访问限制在这些边界内，从而减少了冲突机会。而REDUCE-1依赖步长并移除了这些边界，导致库冲突和进程串行化。
由于每个内存库每个周期只能处理一次访问，当多个访问指向同一库时必须进行串行化处理，这有效降低了内存操作的吞吐量。这种串行化抵消了通过消除warp发散获得的部分性能提升，在较大线程块中可能成为显著的性能瓶颈。
REDUCE-2：顺序寻址优化
为了解决内存访问冲突问题，这个版本采用了更高效的寻址技术。与让线程访问间隔较远元素的交替寻址不同，顺序寻址让每个线程处理连续的数据元素。
在1024元素、每块256线程的示例中，线程0将访问连续的元素0、1、2、3，而非间隔较远的0、256、512、768。线程0依次组合元素0和1、然后处理元素2，以此类推进行递归处理。这种方法充分利用了空间局部性原理，通过提升缓存效率来避免库冲突。该算法具有线性特征，最大限度地减少了会增加等待时间的同步需求。
顺序寻址 | 来源：NVIDIA
这种变化通过将内存访问模式与GPU对连续内存访问的偏好更紧密对齐，显著改善了内存访问效率。通过访问相邻的内存位置，REDUCE-2降低了缓存未命中和内存库冲突的概率，提升了内存带宽利用效率，并改善了归约操作的整体性能。
// 归约 2 – 顺序寻址 __global__ void reduce2(int *g_in_data, int *g_out_data){ extern __shared__ int sdata[]; // 存储在共享内存中 // 每个线程从全局内存加载一个元素到共享内存 unsigned int tid = threadIdx.x; unsigned int i = blockIdx.x * blockDim.x + threadIdx.x; sdata[tid] = g_in_data[i]; __syncthreads(); // 归约方法 -- 在共享内存中进行 for(unsigned int s = blockDim.x/2; s > 0; s >>= 1){ // REDUCE2 -- 查看上面的反向循环 if (tid < s){ // 然后，我们检查线程ID来进行计算 sdata[tid] += sdata[tid + s]; } __syncthreads(); } if (tid == 0){ g_out_data[blockIdx.x] = sdata[0]; } }该方法的主要技术创新包括用反向循环结构结合基于线程ID的索引替换了步长索引机制，从根本上改变了归约过程中的数据处理方式。反向循环从最高可能步长s = blockDim.x / 2开始归约，每次迭代将步长减半。这意味着线程首先处理待求和数据间的最大间隙，快速减少需要处理的数据总量。基于线程ID的索引使每个线程使用其ID来访问连续的数据点对而非分散的数据点，简化了访问模式并最小化了内存延迟。随着步长的减小，线程组合相邻元素，优化了内存使用并提升了数据吞吐量。
结果
REDUCE-2 结果
线程利用率问题
这种方法基本上解决了内存冲突问题。在解决了明显的计算和内存问题后，需要进一步提升算法的智能化程度以获得更好的性能表现。
当前面临的主要问题是在第一个循环迭代中有一半的线程处于空闲状态，这造成了资源浪费并未能充分利用GPU的计算能力。在1024元素的示例中，循环第一次迭代时s=blockDim.x/2（即s=512），条件if (tid < s)将活跃计算限制在块的前512个线程。这意味着虽然这512个线程在积极地对元素对求和（例如sdata[tid]与sdata[tid + 512]），剩余的512个线程却处于空闲状态，对计算没有任何贡献。这种在每个后续迭代中将活跃线程数量减半的模式持续到归约完成，从512减少到256，然后是128、64、32等。这种快速的线程数量衰减导致GPU能力的显著浪费，特别是在初始迭代中只有一小部分可用线程参与工作。
解决方案是在数据加载到共享内存的同时进行第一次计算操作。
REDUCE-3：加载时预归约
为了充分利用空闲线程并提升计算效率，在从全局内存向共享内存加载元素的同时执行第一次计算操作。这种方法能够在加载过程中将两个元素归约为一个，从而将需要处理的数据块数量减半。
具体实现中，在1024元素、256线程的配置下，每个线程将前两个元素的和加载到共享内存中。线程0处理元素0和1，线程1处理元素2和3，以此类推。这样可以将数据块数量和共享内存长度都减半到512。代码的其余部分与REDUCE-2完全相同，这意味着第一次迭代仍然会激活512个线程开始归约操作，因为s=blockDim.x/2 = 512。这种方法让更多线程参与有效工作，避免了计算资源的浪费。
// 归约 3 – 加载时首次加法 __global__ void reduce3(int *g_in_data, int *g_out_data){ extern __shared__ int sdata[]; // 存储在共享内存中 // 每个线程从全局内存加载一个元素到共享内存 unsigned int tid = threadIdx.x; unsigned int i = blockIdx.x*(blockDim.x*2) + threadIdx.x; sdata[tid] = g_in_data[i] + g_in_data[i+blockDim.x]; __syncthreads(); // 归约方法 -- 在共享内存中进行 for(unsigned int s = blockDim.x/2; s > 0; s >>= 1){ // 查看上面的反向循环 if (tid < s){ // 然后，我们检查tid来进行计算 sdata[tid] += sdata[tid + s]; } __syncthreads(); } if (tid == 0){ g_out_data[blockIdx.x] = sdata[0]; } }实现过程中包含三个关键技术变更。首先，在从全局内存加载元素到共享内存时同时进行初始归约步骤：sdata[tid] = g_in_data[i] + g_in_data[i+blockDim.x]。其次，修改索引i的计算方式为unsigned int i = blockId.x * (blockDim.x*2) + threadId.x，因为每个线程现在同时处理两个输入，需要将每个块覆盖的有效索引范围扩大一倍。最后，在主函数中修改内核调用方式，将执行配置设置为int num_blocks = (n + (2*blockSize) - 1 / (2*blockSize)，这样可以将分配给内核的块数量减半，同时保持代码的正确性。
结果
REDUCE-3 结果
指令开销瓶颈识别
当前方法表现良好，但仍有进一步优化的空间。通过分析性能指标发现，在Tesla T4上约41 GB/s的带宽使用率表明我们并未达到或耗尽带宽上限。另一方面，归约操作具有低算术强度的特征，意味着我们也不受计算能力限制。
由于既不受带宽限制也不受计算限制，还存在第三个潜在瓶颈：指令开销。这包括GPU执行的所有辅助指令，这些指令不直接参与数据加载、存储或归约的主要算术操作。具体包括地址算术运算（计算下一个要加载的地址空间）和循环开销（处理循环逻辑、条件判断和迭代控制）。
针对这种瓶颈的优化策略是循环展开技术。
REDUCE-4：Warp级循环展开
首先分析REDUCE-3中的执行模式以理解优化的必要性。在1024元素的示例中，经过初始的元素对加载和相加后，256个线程处理512个元素。此时归约过程中每个线程处理单个元素，线程活跃度呈现递减模式：当s = 256时有256个活跃线程，当s = 128时有128个活跃线程，当s = 64时有64个活跃线程。
关键的优化点出现在s = 32时，此时有32个活跃线程。由于指令在warp内以SIMD方式同步执行，这带来两个重要特性：首先，无需使用__syncthreads()，因为所有线程在同一个warp中以锁步方式工作；其次，无需if (tid < s)条件判断，因为每个线程都需要执行相同的操作。因此可以安全地从这部分代码中移除所有同步命令，显著提升最终归约阶段的执行速度。
// 添加这个函数来帮助展开 __device__ void warpReduce(volatile int* sdata, int tid){ // 目标是让所有warp避免无用的工作 sdata[tid] += sdata[tid + 32]; sdata[tid] += sdata[tid + 16]; sdata[tid] += sdata[tid + 8]; sdata[tid] += sdata[tid + 4]; sdata[tid] += sdata[tid + 2]; sdata[tid] += sdata[tid + 1]; } // 归约 4 – 展开最后的Warp __global__ void reduce4(int *g_in_data, int *g_out_data){ extern __shared__ int sdata[]; // 存储在共享内存中 // 每个线程从全局内存加载一个元素到共享内存 unsigned int tid = threadIdx.x; unsigned int i = blockIdx.x*(blockDim.x*2) + threadIdx.x; sdata[tid] = g_in_data[i] + g_in_data[i+blockDim.x]; __syncthreads(); // 只是将结束限制改为在s = 32之前停止 for(unsigned int s = blockDim.x/2; s > 32; s >>= 1){ // 查看上面的反向循环 if (tid < s){ // 然后，我们检查tid来进行计算 sdata[tid] += sdata[tid + s]; } __syncthreads(); } // 添加这个以在s = 32时使用warpReduce if (tid < 32){ warpReduce(sdata, tid); } if (tid == 0){ g_out_data[blockIdx.x] = sdata[0]; } }实现方法相对简单：在s = 32之前停止主循环，并调用专门的warpReduce函数。该函数包含手写的6次迭代，仅在设备端执行。同时需要使用volatile关键字确保实现的正确性。
结果
REDUCE-4 结果
扩展循环展开策略
这种优化取得了显著的性能提升效果。既然循环展开如此有效，为什么不将这种技术扩展到更多的循环中？
REDUCE-5：完全循环展开
为了进一步扩展展开技术，需要在编译时确定循环的总迭代次数。幸运的是，GPU将线程块大小限制为512个线程，且通常使用2的幂次方配置。因此可以针对固定的块大小进行完全展开，同时保持通用性。CUDA提供的C++模板参数支持使这种优化成为可能。
C++模板技术允许定义具有占位符的函数或类，这些占位符在编译时被具体类型替换。通过使用模板参数来处理blockSize的变化，可以应对不同的展开需求。根据块大小的不同，准备相应的switch case来处理特定的展开要求。完全展开技术消除了大部分归约阶段中不必要的循环和条件判断，最小化了计算开销。
通过编译针对特定块大小（如512、256和128）定制的内核版本，为每个变体优化其特定场景，剥离不必要的操作，最大化内存和计算资源效率。在具体实现中，主函数中将blockSize设置为256以简化方法，同时包含了512、256和128的switch case以展示该方法的灵活性，突出CUDA如何有效利用模板参数来提升不同配置下的性能。
// 添加这个函数来帮助展开并添加模板 template <unsigned int blockSize> __device__ void warpReduce(volatile int* sdata, int tid){ if(blockSize >= 64) sdata[tid] += sdata[tid + 32]; if(blockSize >= 32) sdata[tid] += sdata[tid + 16]; if(blockSize >= 16) sdata[tid] += sdata[tid + 8]; if(blockSize >= 8) sdata[tid] += sdata[tid + 4]; if(blockSize >= 4) sdata[tid] += sdata[tid + 2]; if(blockSize >= 2) sdata[tid] += sdata[tid + 1]; } // 归约 5 – 完全展开 template <unsigned int blockSize> __global__ void reduce5(int *g_in_data, int *g_out_data){ extern __shared__ int sdata[]; // 存储在共享内存中 // 每个线程从全局内存加载一个元素到共享内存 unsigned int tid = threadIdx.x; unsigned int i = blockIdx.x*(blockDim.x*2) + threadIdx.x; sdata[tid] = g_in_data[i] + g_in_data[i+blockDim.x]; __syncthreads(); // 分步执行归约，减少线程同步 if (blockSize >= 512) { if (tid < 256) { sdata[tid] += sdata[tid + 256]; } __syncthreads(); } if (blockSize >= 256) { if (tid < 128) { sdata[tid] += sdata[tid + 128]; } __syncthreads(); } if (blockSize >= 128) { if (tid < 64) { sdata[tid] += sdata[tid + 64]; } __syncthreads(); } if (tid < 32) warpReduce<blockSize>(sdata, tid); if (tid == 0){ g_out_data[blockIdx.x] = sdata[0]; } }相应地，需要修改内核调用方式以支持完全展开：
// 完全展开所需 // 启动内核并同步线程 switch (blockSize) { case 512: reduce6<512><<<num_blocks, 512, 512 * sizeof(int)>>>(dev_input_data, dev_output_data, n); break; case 256: reduce6<256><<<num_blocks, 256, 256 * sizeof(int)>>>(dev_input_data, dev_output_data, n); break; case 128: reduce6<128><<<num_blocks, 128, 128 * sizeof(int)>>>(dev_input_data, dev_output_data, n); break; }实现方式与REDUCE-4相似，主要变化是将blockSize作为编译时确定的模板参数。通过包含条件语句处理不同的blockSize值，以及使用switch语句根据这些值调用相应的内核版本。
结果
REDUCE-5 结果
灵活性与可扩展性的权衡
虽然Reduce5通过为已知块大小完全展开循环来提升效率，但这种方法缺乏灵活性且难以扩展。完全展开技术严重依赖编译时优化，将内核限制为固定的块大小配置。当数据大小与块配置不完全匹配时，可能导致效率降低和GPU资源的次优利用。此外，为每个块大小管理多个内核版本增加了开发复杂度，限制了对变化工作负载的动态适应能力，使其在输入大小变化较大的通用应用中实用性受限。
因此，需要借鉴REDUCE-3中加载时预归约的思想，尝试执行更多的加法操作而非仅限于第一次加法。
REDUCE-6：多重归约与线程级并行
REDUCE-6通过引入"算法级联"的动态方法来解决REDUCE-5中的刚性和可扩展性问题。该方法让每个线程在更广泛的块大小范围内执行多次加法操作，有效减少了对特定块配置的依赖。这种灵活性使算法能够更平滑地适应不同的数据规模，在更广泛的场景中优化资源利用率。
通过结合顺序和并行归约技术，REDUCE-6最小化了延迟并最大化了吞吐量，特别适用于具有高内核启动开销和多样化工作负载的环境。基于Brent定理的工作分配策略确保每个线程在整个归约过程中都能以最优方式贡献计算能力，在与硬件能力有效匹配的同时保持成本效率。
该方法的核心思想是每个线程在同步屏障之前处理多个元素对，从而在更多计算中摊销同步成本，提升整体性能表现。
最终优化内核实现
// 添加这个函数来帮助展开并添加模板 template <unsigned int blockSize> __device__ void warpReduce(volatile int* sdata, unsigned int tid){ if(blockSize >= 64) sdata[tid] += sdata[tid + 32]; if(blockSize >= 32) sdata[tid] += sdata[tid + 16]; if(blockSize >= 16) sdata[tid] += sdata[tid + 8]; if(blockSize >= 8) sdata[tid] += sdata[tid + 4]; if(blockSize >= 4) sdata[tid] += sdata[tid + 2]; if(blockSize >= 2) sdata[tid] += sdata[tid + 1]; } // 归约 6 – 多重加法/线程 template <int blockSize> __global__ void reduce6(int *g_in_data, int *g_out_data, unsigned int n){ extern __shared__ int sdata[]; // 存储在共享内存中 // 每个线程从全局内存加载一个元素到共享内存 unsigned int tid = threadIdx.x; unsigned int i = blockIdx.x*(blockSize*2) + tid; unsigned int gridSize = blockDim.x * 2 * gridDim.x; sdata[tid] = 0; while(i < n) { sdata[tid] += g_in_data[i] + g_in_data[i + blockSize]; i += gridSize; } __syncthreads(); // 分步执行归约，减少线程同步 if (blockSize >= 512) { if (tid < 256) { sdata[tid] += sdata[tid + 256]; } __syncthreads(); } if (blockSize >= 256) { if (tid < 128) { sdata[tid] += sdata[tid + 128]; } __syncthreads(); } if (blockSize >= 128) { if (tid < 64) { sdata[tid] += sdata[tid + 64]; } __syncthreads(); } if (tid < 32) warpReduce<blockSize>(sdata, tid); if (tid == 0){ g_out_data[blockIdx.x] = sdata[0]; } }关键的技术创新体现在while循环中，每个线程直接在共享内存中执行多次加法操作。该循环设计为每次迭代聚合两个数据元素，有效地将必要操作数量和与全局内存的交互频次减半。线程从全局内存加载数据并添加到先前的累积值中，然后按总线程数的两倍向前跳跃，确保在下一次迭代中处理另一对元素。这种模式显著减少了每个线程在任何时刻需要处理的数据量，最大化了可用带宽的利用率并最小化了访问延迟。
性能评估与比较分析
REDUCE-6 结果
所有优化技术的性能对比
与NVIDIA基准的对比分析
本实现与NVIDIA官方实现的主要差异在于硬件平台的不同。NVIDIA的研讨会使用GeForce 8800，而本研究采用Tesla T4。由于Tesla T4具有更优化的架构，使得初始实现就具备了更好的性能基础，但这也意味着性能提升的空间相对有限。虽然无法复现NVIDIA展示的戏剧性加速效果，但成功展示了持续的优化进程和GPU峰值性能的逐步提升。
总结
基于本次优化实践，总结出CUDA内核优化的核心要点如下。
首先，深入理解核心性能特征是优化的基础，包括内存合并访问、分支发散管理、内存库冲突解决以及延迟隐藏技术的应用。其次，充分利用性能指标进行瓶颈识别，通过计算和内存性能指标判断内核是计算受限还是内存受限。第三，系统化地识别瓶颈来源，确定性能限制是由内存访问、计算能力还是指令开销造成的。第四，采用渐进式算法优化策略，先优化基础算法逻辑，再进行循环展开等高级优化。最后，灵活运用模板参数技术进行代码生成的精细调优，确保为不同块大小配置提供最优的实现方案。
作者：Rimika Dhara
