CUDA并行数组批量求和的方案选型与核函数参数调优

需求背景 #

• 待处理的连续float数组总规模为M*N（示例参数为M=100、N=2000），需要对每段长度为N的连续子数组独立求和，最终输出长度为M的求和结果数组。
• 当前已验证两类基础实现路径：
  1. 基于cuBLAS的循环调用方案：逐段循环M次，每段单独调用cublasSasum接口完成单段数组的求和计算
  2. 自定义朴素核函数方案：每个线程负责一段数组，通过循环遍历累加段内所有元素得到求和结果

已观测的性能特征 #

• cuBLAS方案的耗时核心由循环次数M决定，单段长度N对总耗时的影响极低
• 自定义朴素核函数的性能受N值影响极强：
  • 当N<5000时，运行速度远高于cuBLAS方案
  • 随N增大耗时快速上升，实测对应耗时如下：

    N值	实测耗时
    4000	254ms
    10000	422ms
    40000	1365ms
    80000	4361ms
    1e6	5399ms
    2e6	10635ms

  • 当N足够大时，自定义朴素核函数的性能会远低于cuBLAS方案

核心问题解答 #

1. 跨场景的方案选型依据 #

• 不需要在每台设备上做全量参数扫测，可通过两类实现的固定成本结构推导通用选型阈值：
  • cuBLAS方案的核心开销是单次函数调用的固定启动成本，单次cublasSasum的主机到设备调用固定开销约为5~10us，总开销近似为M * 单次调用固定开销。库内部已针对大规模规约计算做了访存合并、并行优化，因此N值变化对总耗时的影响很小。
  • 朴素自定义核函数的核心开销是单线程串行累加的访存延迟：每个线程串行读取N个元素，无并行规约优化，当N较大时，单线程需要等待N次全局内存返回，无法通过并行调度覆盖访存延迟。
• 通用选型参考：
  • 当N < 2000~5000（不同GPU的该阈值随内存带宽波动）时，选择自定义核函数，避免M次cuBLAS调用的固定启动开销
  • 当N大于上述阈值后，直接使用cuBLAS方案，可利用库内置的并行优化能力，避免单线程串行累加的性能损耗
• 跨设备适配时，仅需要在首次运行时用10组左右的M/N样本做一次阈值标定，不需要全量扫测所有参数组合。

2. 自定义核函数的gridSize/blockSize配置逻辑 #

• 观测到<<<M,1>>>配置远快于<<<1,M>>>的核心原因和寄存器占用无关，本质是访存效率差异：
  • 当配置为<<<1,M>>>时，所有线程都在同一个block内，每个线程负责的内存地址是in + N*i，相邻线程的地址间隔为N个float（例如N=1e5时，相邻线程地址间隔400KB），完全无法形成合并访存，每个内存请求都要单独发起，访存效率极低。
  • 当配置为<<<M,1>>>时，每个线程在独立的block内，相邻block对应的线程虽然地址间隔同样为N个float，但CUDA内存调度会把不同SM上的访存请求打散，结合L2缓存的预取能力，相比单block内的零散访存效率可提升数倍。
• 针对单线程负责一段数组的朴素核函数，推荐配置规则为：blockSize取32~128区间，gridSize取ceil(M / blockSize)，不需要刻意匹配SM的硬件规模。

3. 大blockSize场景的性能下降原因 #

• blockSize取1024（达到理论满占用率）时性能更差的核心原因是访存延迟覆盖不足：
  • 当blockSize为1024时，受每个SM的硬件资源限制，单个SM上仅能驻留1~2个block。此时SM上的总并行线程数虽然达到了硬件理论满额，但每个线程都是串行访问自己负责的N个元素，线程在等待全局内存返回的间隙，没有足够的其他可调度线程把访存延迟隐藏起来。
  • 当blockSize为10~200时，每个SM上可以驻留更多的block，总并行线程数虽然没有达到理论满值，但每个线程的串行访存间隙可以被同一SM上其他block的线程计算充分覆盖，最终的实际内存带宽利用率反而更高。
• 注意：上述结论仅针对单线程负责整段累加的朴素设计，如果调整为段内多线程并行规约（例如每个block负责一段数组的求和），blockSize的最优值会回到256~512区间。

参考实现代码 #

// cuBLAS实现路径
for (int j = 0;j<M;j++)
    cublasStatus = cublasSasum(cublasHandle,N,d_in+N*j,1,d_out+j);

// 自定义求和核函数
__global__ void getSum(int M, int N, float* in, float * out)
{
    int i = threadIdx.x + blockIdx.x * blockDim.x;
    if(i<M){
        float tmp = 0;
        for(int ii = 0; ii<N; ii++){
            tmp += *(in+N*i+ii);
        }
        out[i] = tmp;
    }
}

// 不同网格配置的性能对比
getSum<<<M,1>>>(M, N, d_in, d_out); // 实测速度更快
getSum<<<1,M>>>(M, N, d_in, d_out); 

// blockSize扫测的计时逻辑
cudaEventRecord(start, 0);
// 测试blockSize取值范围为1~1024
int gridSize = (M + blockSize - 1) / blockSize;
getSum<<<gridSize, blockSize>>>(M, N, d_in, d_out);
cudaEventRecord(stop, 0);
cudaEventSynchronize(stop);
cudaEventElapsedTime(&time, start, stop);

内容的提问来源于stack exchange，提问作者Wei Zhang