CUDA归约计算结果异常：4M元素求和结果不符求助

问题现象 #

当使用128线程块对N=1<<22（4194304个元素）的全1数组求和时，得到结果4210432，与预期值4194304不符；将线程块大小改为1024后结果正确。较小的N（1<<10到1<<21）使用128线程块时结果正常。

问题根源 #

归约的后续轮次中，当当前待求和数组的大小current_size小于线程块大小THREADS_PER_BLOCK时，核函数中所有线程都会尝试访问全局内存：

• 对于tid >= current_size的线程，i = blockIdx.x * blockDim.x + threadIdx.x会超出输入数组的有效范围，导致全局内存越界访问，读取到未初始化的垃圾值。
• 这些垃圾值会被写入共享内存，并在后续归约步骤中被累加，最终导致求和结果错误。

以N=1<<22、THREADS_PER_BLOCK=128为例：

1. 第一次归约后，current_size = 4194304 / 128 = 32768
2. 第二次归约后，current_size = 32768 / 128 = 256
3. 第三次归约后，current_size = 256 / 128 = 2
4. 第四次归约时，blocks = (2 + 127)/128 = 1，此时128个线程中，tid >=2的线程访问全局内存越界，引入垃圾值并参与求和，导致结果错误。

而THREADS_PER_BLOCK=1024时，后续轮次的越界访问恰好读取到0值（属于巧合），因此结果正确，但代码本身仍存在隐患。

修复方案 #

在核函数中添加判断：当线程对应的全局索引i超出当前输入数组的有效范围时，将共享内存对应位置设为0，避免垃圾值参与求和。同时，需要将当前数组大小作为参数传入核函数：

修改后的代码 #

#include <cuda_runtime.h>
#include <iostream>
#include <cstdlib>

#define N (1 << 22)  // 4M ints
#define THREADS_PER_BLOCK 128

__global__ void reduce0(int* g_idata, int* g_odata, int current_size) {
    extern __shared__ int sdata[];

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

    // 处理越界线程，赋值为0避免引入垃圾值
    if (i < current_size)
        sdata[tid] = g_idata[i];
    else
        sdata[tid] = 0;
    __syncthreads();

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

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

int main() {
    const int size = N * sizeof(int);
    int* h_data = (int*)malloc(size);

    for (int i = 0; i < N; ++i)
        h_data[i] = 1;
    
    int* d_idata = nullptr;
    cudaMalloc(&d_idata, size);
    cudaMemcpy(d_idata, h_data, size, cudaMemcpyHostToDevice);

    int* d_odata = nullptr;

    int current_size = N;
    int blocks = 0;
    int* input = d_idata;
    int* output = nullptr;

    while (current_size > 1) {
        blocks = (current_size + THREADS_PER_BLOCK - 1) / THREADS_PER_BLOCK;

        cudaMalloc(&d_odata, blocks * sizeof(int));

        // 传入current_size参数给核函数
        reduce0<<<blocks, THREADS_PER_BLOCK, THREADS_PER_BLOCK * sizeof(int)>>>(input, d_odata, current_size);

        if (input != d_idata) cudaFree(input);
        input = d_odata;
        current_size = blocks;
    }

    int h_result = 0;
    cudaMemcpy(&h_result, input, sizeof(int), cudaMemcpyDeviceToHost);

    std::cout << "GPU total sum = " << h_result << std::endl;

    cudaFree(input);
    if (d_idata != input) cudaFree(d_idata);
    free(h_data);

    return 0;
}

额外优化建议 #

原归约核函数使用的tid % (s*2) == 0的分支方式会导致线程束分化，降低执行效率。可以改为跨步归约的方式，让每个线程处理不重叠的元素，避免分支：

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

这种方式的线程执行更均匀，能显著提升归约性能。

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