计算 CUDA 中嵌套循环的索引

Question

我正在尝试学习 CUDA，但对计算线程索引有点困惑。假设我有一个正在尝试并行化的循环：

...
for(int x = 0; x < DIM_x; x++){
    for(int y = 0; y < DIM_y; y++){
        for(int dx = 0; dx < psize; dx++){
            array[y*DIM_x + x + dx] += 1;
        }
    }
}

在 PyCUDA 中，我设置：

block = (8, 8, 8)
grid = (96, 96, 16)

我见过的大多数并行化循环示例都是这样计算线程索引的：

int x = blockIdx.x * blockDim.x + threadIdx.x;
int y = blockIdx.y * blockDim.y + threadIdx.y;
int dx = blockIdx.z * blockDim.z + threadIdx.z;

if (x >= DIM_x || y >= DIM_y || dx >= psize)
    return;

atomicAdd(&array[y*DIM_x + x + dx], 1)

DIM_x=580，DIM_y=550，psize=50

但是，如果我打印 x，我看到创建了具有相同线程 Id 的多个线程，最终结果是错误的。

相反，如果我使用这个（3D 块的 3D 网格）：

int blockId = blockIdx.x + blockIdx.y * gridDim.x
              + gridDim.x * gridDim.y * blockIdx.z;

int x = blockId * (blockDim.x * blockDim.y * blockDim.z)
        + (threadIdx.z * (blockDim.x * blockDim.y))
        + (threadIdx.y * blockDim.x) + threadIdx.x;

它修复了 x 的多个相同线程 ID 问题，但我不确定如何并行化 y 和 dx。

如果有人能帮助我理解哪里出了问题，并告诉我并行化循环的正确方法，我将不胜感激。

Answer 1

在具有 3D 块的 3D 网格中，线程 ID 为：

    unsigned long blockId = blockIdx.x 
             + blockIdx.y * gridDim.x 
             + gridDim.x * gridDim.y * blockIdx.z; 
    unsigned long threadId = blockId * (blockDim.x * blockDim.y * blockDim.z)
              + (threadIdx.z * (blockDim.x * blockDim.y))
              + (threadIdx.y * blockDim.x)
              + threadIdx.x;

不是您计算的x。 x 只是该 3D 矩阵的 x 索引。

有一个很好的备忘单in this blog

Answer 2

However, if I print x, I see that multiple threads with the same thread Id are created, and the final result is wrong.

您在多维网格中看到具有相同 x 线程 ID 的多个线程是正常的，因为在您的主机代码中观察到具有相同 x 值的循环的多次迭代也是正常的.如果结果错误，则与您显示的任何代码无关，即：

#include <vector>
#include <thrust/device_vector.h>
#include <thrust/copy.h>
#include <assert.h>

void host(int* array, int DIM_x, int DIM_y, int psize)
{
    for(int x = 0; x < DIM_x; x++){
        for(int y = 0; y < DIM_y; y++){
            for(int dx = 0; dx < psize; dx++){
                array[y*DIM_x + x + dx] += 1;
            }
        }
    }
}


__global__
void kernel(int* array, int DIM_x, int DIM_y, int psize)
{
    int x = blockIdx.x * blockDim.x + threadIdx.x;
    int y = blockIdx.y * blockDim.y + threadIdx.y;
    int dx = blockIdx.z * blockDim.z + threadIdx.z;

    if (x >= DIM_x || y >= DIM_y || dx >= psize)
        return;

    atomicAdd(&array[y*DIM_x + x + dx], 1);
}

int main()
{
    dim3 block(8, 8, 8);
    dim3 grid(96, 96, 16);

    int DIM_x = 580, DIM_y = 550, psize = 50;

    std::vector<int> array_h(DIM_x * DIM_y * psize, 0);
    std::vector<int> array_hd(DIM_x * DIM_y * psize, 0);
    thrust::device_vector<int> array_d(DIM_x * DIM_y * psize, 0);

    kernel<<<grid, block>>>(thrust::raw_pointer_cast(array_d.data()), DIM_x, DIM_y, psize);
    host(&array_h[0], DIM_x, DIM_y, psize);

    thrust::copy(array_d.begin(), array_d.end(), array_hd.begin());
    cudaDeviceSynchronize();

    for(int i=0; i<DIM_x * DIM_y * psize; i++) {
        assert( array_h[i] == array_hd[i] ); 
    }

    return 0;
}

编译时运行

$ nvcc -arch=sm_52 -std=c++11 -o looploop loop_the_loop.cu 
$ cuda-memcheck ./looploop 
========= CUDA-MEMCHECK
========= ERROR SUMMARY: 0 errors

未发出任何错误并通过了针对您问题中主机代码的所有元素的检查。

如果您得到不正确的结果，可能是您在运行内核之前初始化设备内存有问题。否则我看不出您所显示的代码如何发出不正确的结果。

一般来说，像您的代码那样执行大量原子内存事务并不是在 GPU 上执行计算的最佳方式。使用非原子事务可能需要依赖其他 a priori 有关问题结构的信息（例如图形分解或问题写入模式的精确描述）。

计算 CUDA 中嵌套循环的索引

Calculating indices for nested loops in CUDA

cuda

pycuda