c++ - 具有多个数组的共享内存的 CUDA 缩减

Question

我目前正在使用以下 Reduction 函数通过 CUDA 对数组中的所有元素求和:

__global__ void reduceSum(int *input, int *input2, int *input3, int *outdata, int size){
    extern __shared__ int sdata[];

    unsigned int tID = threadIdx.x;
    unsigned int i = tID + blockIdx.x * (blockDim.x * 2);
    sdata[tID] = input[i] + input[i + blockDim.x];
    __syncthreads();

    for (unsigned int stride = blockDim.x / 2; stride > 32; stride >>= 1)
    {
        if (tID < stride)
        {
            sdata[tID] += sdata[tID + stride];
        }
        __syncthreads();
    }

    if (tID < 32){ warpReduce(sdata, tID); }

    if (tID == 0)
    {
        outdata[blockIdx.x] = sdata[0];
    }
}

但是，正如您从函数参数中看到的那样，我希望能够在一个缩减函数中对三个单独的数组求和。很明显，一个简单的方法是启动内核 3 次，每次都传递一个不同的数组，这当然可以正常工作。不过我现在只是把它写成一个测试内核，真正的内核最终会采用一个结构数组，我需要对每个结构的所有 X、Y 和 Z 值执行加法，这就是为什么我需要将它们全部汇总到一个内核中。

我已经为所有三个数组初始化并分配了内存

    int test[1000];
    std::fill_n(test, 1000, 1);
    int *d_test;

    int test2[1000];
    std::fill_n(test2, 1000, 2);
    int *d_test2;

    int test3[1000];
    std::fill_n(test3, 1000, 3);
    int *d_test3;

    cudaMalloc((void**)&d_test, 1000 * sizeof(int));
    cudaMalloc((void**)&d_test2, 1000 * sizeof(int));
    cudaMalloc((void**)&d_test3, 1000 * sizeof(int));

我不确定我应该为这种内核使用什么网格和 block 维度，我也不完全确定如何修改缩减循环以按照我的需要放置数据，即输出数组:

Block 1 Result|Block 2 Result|Block 3 Result|Block 4 Result|Block 5 Result|Block 6 Result|

      Test Array 1 Sums              Test Array 2 Sums            Test Array 3 Sums           

我希望这是有道理的。或者有没有更好的方法只有一个归约函数但能够返回 Struct.X、Struct.Y 或 struct.Z 的总和？

结构如下:

template <typename T>
struct planet {
    T x, y, z;
    T vx, vy, vz;
    T mass;
};

我需要将所有 VX 相加并存储，将所有 VY 相加并存储，将所有 VZ 相加并存储。

Answer 1

Or is there a better way to have only one reduction function but be able to return the summation of Struct.X, Struct.Y or struct.Z?

通常加速计算的主要焦点是速度。 GPU 代码的速度(性能)通常在很大程度上取决于数据存储和访问模式。因此，尽管正如您在问题中指出的那样，我们可以通过多种方式实现解决方案，但让我们专注于应该相对较快的事情。

像这样的归约没有太多的算术/操作强度，因此我们对性能的关注主要围绕数据存储以实现高效访问。当访问全局内存时，GPU 通常会以大块的形式进行——32 字节或 128 字节的 block 。为了有效地使用内存子系统，我们希望在每个请求中使用所请求的所有 32 或 128 个字节。

但是你的结构隐含的数据存储模式:

template <typename T>
struct planet {
    T x, y, z;
    T vx, vy, vz;
    T mass;
};

几乎排除了这一点。对于此问题，您关心 vx、vy 和 vz。这 3 个项目在给定的结构(元素)中应该是连续的，但是在这些结构的数组中，它们将被其他结构项目的必要存储分开，至少:

planet0:       T x
               T y
               T z               ---------------
               T vx      <--           ^
               T vy      <--           |
               T vz      <--       32-byte read
               T mass                  |
planet1:       T x                     |
               T y                     v
               T z               ---------------
               T vx      <--
               T vy      <--
               T vz      <--
               T mass
planet2:       T x
               T y
               T z
               T vx      <--
               T vy      <--
               T vz      <--
               T mass

(为了举例，假设T是float)

这指出了 GPU 中结构数组 (AoS) 存储格式的一个主要缺点。由于 GPU 的访问粒度(32 字节)，从连续结构访问相同元素是低效的。在这种情况下，通常的性能建议是将 AoS 存储转换为 SoA(数组结构):

template <typename T>
struct planets {
    T x[N], y[N], z[N];
    T vx[N], vy[N], vz[N];
    T mass[N];
};

以上只是一个可能的示例，可能不是您实际使用的示例，因为该结构没有什么用处，因为我们只有一个结构用于 N 行星。关键是，现在当我访问连续行星的 vx 时，各个 vx 元素在内存中都是相邻的，所以 32 字节的读取给了我 32 字节的 vx 数据，没有浪费或未使用的元素。

经过这样的改造，从代码组织的角度来看，归约问题又变得相对简单了。您可以使用与单个数组缩减代码基本相同的代码，连续调用 3 次，或者直接扩展内核代码以独立处理所有 3 个数组。 “三合一”内核可能看起来像这样:

template <typename T>
__global__ void reduceSum(T *input_vx, T *input_vy, T *input_vz, T *outdata_vx, T *outdata_vy, T *outdata_vz, int size){
    extern __shared__ T sdata[];

    const int VX = 0;
    const int VY = blockDim.x;
    const int VZ = 2*blockDim.x;

    unsigned int tID = threadIdx.x;
    unsigned int i = tID + blockIdx.x * (blockDim.x * 2);
    sdata[tID+VX] = input_vx[i] + input_vx[i + blockDim.x];
    sdata[tID+VY] = input_vy[i] + input_vy[i + blockDim.x];
    sdata[tID+VZ] = input_vz[i] + input_vz[i + blockDim.x];
    __syncthreads();

    for (unsigned int stride = blockDim.x / 2; stride > 32; stride >>= 1)
    {
        if (tID < stride)
        {
            sdata[tID+VX] += sdata[tID+VX + stride];
            sdata[tID+VY] += sdata[tID+VY + stride];
            sdata[tID+VZ] += sdata[tID+VZ + stride];
        }
        __syncthreads();
    }

    if (tID < 32){ warpReduce(sdata+VX, tID); }
    if (tID < 32){ warpReduce(sdata+VY, tID); }
    if (tID < 32){ warpReduce(sdata+VZ, tID); }

    if (tID == 0)
    {
        outdata_vx[blockIdx.x] = sdata[VX];
        outdata_vy[blockIdx.x] = sdata[VY];
        outdata_vz[blockIdx.x] = sdata[VZ];
    }
}

(在浏览器中编码 - 未经测试 - 只是您作为“引用内核”显示的内容的扩展)

上述 AoS -> SoA 数据转换也可能在代码的其他地方带来性能优势。由于建议的内核将同时处理 3 个数组，因此网格和 block 维度应该与您在单数组情况下用于引用内核的维度完全相同。共享内存存储需要每个 block 增加(三倍)。