CUDA核函数中__shared__内存引发计算错误的原因排查

在CUDA核函数中将3D数组XY维度值p_x、p_y存入__shared__内存数组size_t size_p[2]后，使用共享内存值计算出现错误，但直接用参数则结果正常。以下是问题详情、根因分析及正确的共享内存使用方案：

出错的核函数代码 #

__global__ void element_wise_add(
    float* p,
    size_t p_x,
    size_t p_y,
    size_t p_z,
    float* B,
    size_t B_x,
    size_t B_y,
    size_t B_z,
    unsigned int path_x,
    unsigned int path_y,
    unsigned int path_z,
    const float scalar) // try making this in __shared__ memory
{

    int index = blockIdx.x * blockDim.x + threadIdx.x;


    __shared__ size_t size_p[2], size_B[2];

    if (index == 0)
    {
        size_p[0] = p_x;
        size_p[1] = p_y;
        size_B[0] = B_x;
        size_B[1] = B_y;
        
    }
    
    __syncthreads();
    if (index == 100)
        printf("%d == %d == %d == %d", p_x, p_y, size_p[0], size_p[1]);
    if (index < size_p[0] * size_p[1] * p_z)
    {
        //Get ijk indices from each index
        unsigned int k = index / (p_x * p_y);
        index -= k * p_x * p_y;
        unsigned int j = index / p_x; //maybe here yLen 
        index -= j * p_x;
        unsigned int i = index / 1;

         

        
        B[arrayIndex(i+path_x, j+path_y, k+path_z, B_x, B_y)] += scalar*p[arrayIndex(i, j, k, p_x, p_y)];

        //index = arrayIndex(i + path_x, j + path_y, k + path_z, size_B[0], size_B[1]);
        //int index_B = arrayIndex(i, j, k, size_p[0], size_p[1]);

        //atomicAdd((B + index), scalar * p[index_B]); // make arrayIndex function a preprocessor micro for speed
    }
}

正确的核函数代码 #

__global__ void element_wise_add(
    float* p,
    size_t p_x,
    size_t p_y,
    size_t p_z,
    float* B,
    size_t B_x,
    size_t B_y,
    size_t B_z,
    unsigned int path_x,
    unsigned int path_y,
    unsigned int path_z,
    const float scalar) // try making this in __shared__ memory
{
        
    int index = blockIdx.x * blockDim.x + threadIdx.x;
    


    if (index < p_x * p_y * p_z) 
    {
        //Get ijk indices from each index
        unsigned int k = index / (p_x * p_y);
        index -= k * p_x * p_y;
        unsigned int j = index / p_x; //maybe here yLen 
        index -= j * p_x;
        unsigned int i = index / 1;

    

        B[arrayIndex(i+path_x, j+path_y, k+path_z, B_x, B_y)] += scalar*p[arrayIndex(i, j, k, p_x, p_y)];

        
    }
}

最小驱动代码 #

__host__ __device__ int arrayIndex(int x, int y, int z, int height, int width) {
    return x + y * height + z * height * width;
}


void print_3d_serial_array(float* ptr, size_t X, size_t Y, size_t Z);


void kernel_sample_driver_()
{
    const int Nx = 10;
    const int Ny = 10;
    const int Nz = 10;

    const int px = 10;
    const int py = 2;
    const int pz = 2;


    float a[Nx * Ny * Nz], b[px * py * pz];

    for (size_t k = 0; k < Nz; k++)
    {
        for (size_t j = 0; j < Ny; j++)
        {
            for (size_t i = 0; i < Nx; i++)
            {
                a[arrayIndex(i, j, k, Nx, Ny)] = i + j + k;

            }
        }
    }
    for (size_t k = 0; k < pz; k++)
    {
        for (size_t j = 0; j < py; j++)
        {
            for (size_t i = 0; i < px; i++)
            {
                b[arrayIndex(i, j, k, px, py)] = 1000 * (i + j + k + 1);
            }
        }
    }


    print_3d_serial_array(a, Nx, Ny, Nz);
    print_3d_serial_array(b, px, py, pz);


    gpu::dev_array<float> d_a(Nx * Ny * Nz);
    gpu::dev_array<float> d_b(px * py * pz);

    d_a.set(a, Nx * Ny * Nz);
    d_b.set(b, px * py * pz);


    dim3 threadsPerBlock;
    dim3 blocksPerGrid;
    threadsPerBlock.x = Nx * Ny * Nz;
    threadsPerBlock.y = 1;
    blocksPerGrid.x = ceil(((double)(Nx * Ny * Nz)) / (threadsPerBlock.x));

    element_wise_add <<<blocksPerGrid, threadsPerBlock>>> (d_b.getData(), px, py, pz, d_a.getData(), Nx, Ny, Nz, 0, 1, 1, 1);

    cudaDeviceSynchronize();


    d_a.get(a, Nx * Ny * Nz);

    print_3d_serial_array(a, Nx, Ny, Nz);

}



void print_3d_serial_array(float* ptr, size_t X, size_t Y, size_t Z)
{
    for (size_t k = 0; k < Z; k++)
    {
        int len = 0;
        printf("Array( : , : , %02d) =\n\n", k);
        for (size_t j = 0; j < Y; j++)
        {
            for (size_t i = 0; i < X; i++)
            {
                printf("%3.1f , ", ptr[arrayIndex(i, j, k, X, Y)]);
                }
            std::cout << std::endl;
        }
        std::cout << '\n';
        for (size_t l = 0; l < X; l++)
        {
            std::cout << "-";
        }
        std::cout << '\n';
        std::cout << std::endl;
        }
}

错误根本原因 #

1. 类型不匹配与比较逻辑错误

   • index是有符号int类型，而size_p[0] * size_p[1] * p_z是无符号size_t类型。当无符号数值超过INT_MAX时，有符号与无符号的比较会触发整数提升规则，导致index < 无符号值的判断逻辑失效，大量线程错误进入计算分支，引发内存越界。
   • 用%d（有符号整数格式符）打印无符号size_t类型的size_p值，会造成输出误导，看似共享内存值与参数一致，实际存在类型截断风险。

2. 共享内存初始化的线程块范围问题

   • 代码中用index == 0（全局索引）初始化共享内存，但共享内存是线程块级别的资源。当使用多个线程块时，只有第一个线程块的全局索引0线程会初始化共享内存，其他线程块的共享内存会读取未初始化的垃圾值，导致计算错误。

正确使用共享内存的优化方案 #

修改核函数，解决类型匹配问题，确保共享内存按线程块正确初始化，同时优化线程索引计算：

__global__ void element_wise_add(
    float* p,
    size_t p_x,
    size_t p_y,
    size_t p_z,
    float* B,
    size_t B_x,
    size_t B_y,
    size_t B_z,
    unsigned int path_x,
    unsigned int path_y,
    unsigned int path_z,
    const float scalar)
{
    // 使用无符号整数避免符号比较问题
    unsigned int global_idx = blockIdx.x * blockDim.x + threadIdx.x;
    const size_t p_total = p_x * p_y * p_z;
    if (global_idx >= p_total) return;

    __shared__ size_t size_p[2], size_B[2];

    // 每个线程块的第一个线程初始化共享内存
    if (threadIdx.x == 0)
    {
        size_p[0] = p_x;
        size_p[1] = p_y;
        size_B[0] = B_x;
        size_B[1] = B_y;
    }
    __syncthreads();

    // 使用共享内存计算索引，避免重复读取全局参数
    unsigned int k = global_idx / (size_p[0] * size_p[1]);
    unsigned int rem = global_idx % (size_p[0] * size_p[1]);
    unsigned int j = rem / size_p[0];
    unsigned int i = rem % size_p[0];

    // 显式类型转换，匹配arrayIndex的参数类型
    unsigned int b_idx = arrayIndex(static_cast<int>(i+path_x), 
                                   static_cast<int>(j+path_y), 
                                   static_cast<int>(k+path_z), 
                                   static_cast<int>(size_B[0]), 
                                   static_cast<int>(size_B[1]));
    unsigned int p_idx = arrayIndex(static_cast<int>(i), 
                                   static_cast<int>(j), 
                                   static_cast<int>(k), 
                                   static_cast<int>(size_p[0]), 
                                   static_cast<int>(size_p[1]));
    B[b_idx] += scalar * p[p_idx];
}

同时修改驱动代码的线程配置为CUDA推荐的合理值（比如256线程/块）：

threadsPerBlock.x = 256;
blocksPerGrid.x = (px * py * pz + threadsPerBlock.x - 1) / threadsPerBlock.x;

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