使用共享内存对三维数组的一维求和

Posted 2023-02-16

【中文标题】使用共享内存对三维数组的一维求和

【英文标题】：Summation over one dimension of a three dimensional array using shared memory

【发布时间】：2012-04-01 12:47:05

【问题描述】：

我需要进行如下计算：A[x][y] = sum从 z=0 到 z=nB[x][y][z]+C[x][y][z ]，其中矩阵 A 的维度为 [height][width]，矩阵 B,C 的维度为 [height][width][n]。

值通过以下方式映射到内存：

index = 0;
for (z = 0; z<n; ++z)
    for(y = 0; y<width; ++y)
        for(x = 0; x<height; ++x) 
            matrix[index] = value;
            index++;
        

我想每个块计算一个总和，因为每个块都有自己的共享内存。为了避免数据竞争，我使用 atomicAdd，如下所示：

全局内存中的部分代码：

dim3 block (n, 1, 1);
dim grid (height, width, 1);

内核：

atomicAdd( &(A[blockIdx.x + blockIdx.y*gridDim.y]), 
           B[blockIdx.x + blockIdx.y*gridDim.y+threadIdx.x*blockDim.x*blockDim.y] 
           + C[blockIdx.x + blockIdx.y*gridDim.y+threadIdx.x*blockDim.x*blockDim.y] );

我想使用共享内存来计算总和，然后将此结果复制到全局内存。

我不确定如何使用共享内存来做这部分。在每个块的共享内存中将只存储一个数字（求和结果）。我应该如何将此数字复制到全局内存中 A 矩阵的正确位置？

【参考方案1】：

您可能不需要共享内存或原子内存访问来进行您所询问的求和。如果我理解正确，您的数据按列主要顺序排列，因此逻辑操作是在输出矩阵中的每个矩阵条目都有一个线程，并让每个线程遍历输入矩阵的 z 轴，并在进行时求和。用于此的内核可能类似于：

__global__ void kernel(float *A, const float *B, const float *C, 
        const int width, const int height, const int n)

    int tidx = threadIdx.x + blockDim.x * blockIdx.x;
    int tidy = threadIdx.y + blockDim.y * blockIdx.y;

    if ( (tidx < height) && (tidy < width) ) 
        int stride = width * height;
        int ipos = tidx + tidy * height;

        float * oval = A + ipos;
        float sum = 0.f;
        for(int z=0; z<n; z++, ipos+=stride) 
            sum += B[ipos] + C[ipos];
        
        *oval = sum;
    

这种方法应该最适合具有width * height >= n 的列主要数据。为此使用共享内存没有性能优势，也不需要使用原子内存操作。如果您遇到width * height << n 的问题，尝试对每个求和进行逐块并行缩减可能是有意义的。但是你没有指出问题的典型维度是什么。如果您的问题更像后者，请发表评论，我可以在答案中添加基于缩减的示例内核。

【讨论】：

但是当每个线程遍历输入矩阵的 z 轴时，计算是序列化的，不是吗？

每个线程计算 z=0 中的值，然后移动到 z=1，依此类推。这是序列化的。并且对于每个 z 都必须访问全局内存。我想避免它。

在内核中，A 中的每个输出条目只有一次全局内存写入。每个线程从B 和C 执行n 读取。 B 和 C 中的每个条目都被读取一次。没有更有效的方法来执行此操作。

如果您有 512x512xn 输入，那么我的答案中的代码是最好的方法。如果你有8x8x10^6，那么基于减少的策略是有意义的。

float * oval = A + ipos; ... ; *oval=sum 与A[ipos]=sum 具有相同的效果，只是我们不需要 ipos 的初始值 - 代码只是预先计算了总和需要的地址被存储以便ipos 可以在循环中更改。