CUDA 6D 循环计算

Question

我想在 CUDA（Pascal 架构）中并行化以下 6D 嵌套 for 循环。

const int NX = 250, NY = 250, NZ = 250, NA = 100, NB = 100, NC = 100;
float data_out[NX * NY * NZ];
float data_in[NA * NB * NC];
float datax[NX];
float datay[NY];
float dataz[NZ];

for (int ix = 0; ix < Nx; ix++)
{
    for (int iy = 0; iy < Ny; iy++)
    {
        for (int iz = 0; iz < Nz; iz++)
        {
            float result = 0.0f;
            for (int ia = 0; ia < NA; ia++)
            {
                for (int ib = 0; ib < NB; ib++)
                {
                    for (int ic = 0; ic < NC; ic++)
                    {
                        // some exemplary computation (see kernel)
                    }
                }
            }
            data_out[iz + iy * NZ + ix * (NZ * NY)] = result;
        }
    }
}

目前，我实现了一个执行内部 3D 嵌套 for 循环（循环变量ia、ib、ic）的内核，也就是说，到目前为止我还没有使用并行归约。因此，每个内核计算NA * NB * NC = 1000000值的总和。

编辑：for 循环中的计算已更新以考虑值的任何非线性组合，即不能在 for 循环外计算值

__global__ void testKernel
(
    float *data_out,
    const float *data_in,
    const float *datax,
    const float *datay,
    const float *dataz,
    const int NX,
    const int NY,
    const int NZ,
    const int NA,
    const int NB,
    const int NC
)
{
    int ix = threadIdx.x + blockIdx.x*blockDim.x;
    int iy = threadIdx.y + blockIdx.y*blockDim.y;
    int iz = threadIdx.z + blockIdx.z*blockDim.z;

    if (ix >= NX || iy >= NY || iz >= NZ)
        return;

    float3 xyz = make_float3(datax[ix], datay[iy], dataz[iz]);
    float result = 0.0f;
    for (int ia = 0; ia < NA; ia++)
    {
        for (int ib = 0; ib < NB; ib++)
        {
            for (int ic = 0; ic < NC; ic++)
            {
                // some exemplary nonlinear computation to show memory access
                result += nonlinear_combination(data_in[ic + ib * NC + ia * (NC * NB)], xyz, ia, ib, ic);
            }
        }
    }
    data_out[iz + iy * NZ + ix * (NZ * NY)] = result;
}

int main()
{
    cudaEvent_t start, stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);

    const int NX = 250, NY = 250, NZ = 250, NA = 100, NB = 100, NC = 100;

    float *d_data_out, *d_data_in, *d_datax, *d_datay, *d_dataz;

    cudaMalloc((void**)&d_data_out, NX * NY * NZ * sizeof(float));
    cudaMalloc((void**)&d_data_in, NA * NB * NC * sizeof(float));
    cudaMalloc((void**)&d_datax, NX * sizeof(float));
    cudaMalloc((void**)&d_datay, NY * sizeof(float));
    cudaMalloc((void**)&d_dataz, NZ * sizeof(float));

    dim3 blockSize(8, 8, 8);
    dim3 gridSize(128, 128, 64);

    cudaEventRecord(start);
    testKernel<<<gridSize, blockSize>>>(d_data_out, d_data_in, d_datax, d_datay, d_dataz, NX, NY, NZ, NA, NB, NC);
    cudaEventRecord(stop);

    cudaEventSynchronize(stop);
    float milliseconds = 0;
    cudaEventElapsedTime(&milliseconds, start, stop);

    printf("Elapsed time: %.2f ms\n", milliseconds);

    cudaFree(d_data_out);
    cudaFree(d_data_in);
    cudaFree(d_datax);
    cudaFree(d_datay);
    cudaFree(d_dataz);

    return 0;
}

由于外部 for 循环 (NX * NY * NZ = 15625000) 的总迭代次数已经高于并行线程的总数，因此并行化内部 for 循环以及使用并行缩减是否有任何好处？

另外，如何优化内存访问？确保块中的每个线程都访问数据的相同部分并将这部分数据复制到共享内存可能是有益的，对吗？

Answer 1

一般来说，你的方法看起来是正确的。 15625000 个线程已经非常多了，即使对于具有 10000 个内核的最新 GPU 也是如此。对于他们来说，大约 250000 个线程是可取的。虽然你的块线程部门会浪费很多运行的线程。因为 x * 8 的 128 个线程按 x = 1024 阻塞并且比 NX = 250 少得多。等等。

而且 CUDA 也不允许您在一个块中运行超过 1024 个线程。您可以使用像 (NX, 1, 1) 这样的块大小和像 (1, NY, NZ) 这样的网格大小来节省一些计算。只有理想的块线程大小才能被 32 整除。

对于合并的内存访问，请确保相邻线程访问相邻的内存单元，并且（最好）块对齐到（大约）64 字节。变化最快的线程索引是 x，因此例如在第一个经纱中，线程将有 y 和 z = 0 以及 x = 0，... 31。

您通过对局部变量求和并仅写入一次结果来做对了。

至于减少线程数，这可以节省您的线程初始化。 15625000 次 block int ix = threadIdx.x + blockIdx.x*blockDim.x 的执行；和以下或更少。由于您的内部 3D 循环很大，因此这保证了非常小的收益。

是的，一些额外的计算重组可以避免多次读取 data_in 数组。查看 GPU 示例中的经典矩阵乘法。

我也会尝试循环展开。当然，在确保合并和最小内存访问之后（但也许你不需要共享内存，因为 GPU 会自动使用 SM 内存作为缓存）。在获得第一个工作版本后，您将能够通过 NSight Compute 深入了解代码的有效性。