多阵列共享内存的 CUDA 减少

Question

我目前正在使用以下归约函数将数组中的所有元素与 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];
    }
}

但是，正如您从函数参数中看到的那样，我希望能够在一个归约函数内对三个单独的数组求和。现在显然一个简单的方法是启动内核三次并每次传递一个不同的数组，这当然可以正常工作。我现在只是把它写成一个测试内核，真正的内核最终会接受一个结构数组，我需要对每个结构的所有 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));

我不确定我应该为这种内核使用什么 Grid 和 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

或者有没有更好的方法，只有一个归约函数，但能够返回 Struct.X、Struct.Y 或 struct.Z 的总和？

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

这样的缩减没有太多的算术/运算强度，因此我们对性能的关注主要围绕数据存储以实现高效访问。在访问全局内存时，GPU 通常会以大块的形式进行 - 32 字节或 128 字节块。为了有效利用内存子系统，我们希望在每次请求时都使用所请求的所有 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 个数组，因此网格和块的尺寸应该完全相同，就像您在单数组情况下用于参考内核的尺寸一样。每个块的共享内存存储需要增加（三倍）。

Answer 2

Robert Crovella 给出了一个很好的答案，强调了 AoS 的重要性 -> SoA 布局转换通常可以提高 GPU 的性能，我只想提出一个可能更方便的中间立场。 CUDA 语言提供了一些向量类型，仅用于您描述的目的（请参阅this section of the CUDA programming guide）。

例如，CUDA 定义了 int3，一种存储 3 个整数的数据类型。

 struct int3
 {
    int x; int y; int z;
 };

浮点数、字符、双精度数等存在类似的类型。这些数据类型的优点在于它们可以用一条指令加载，这可能会给您带来小的性能提升。有关此问题的讨论，请参阅 this NVIDIA blog post。对于这种情况，它也是一种更“自然”的数据类型，它可能会使代码的其他部分更容易使用。例如，您可以定义：

struct planets {
    float3 position[N];
    float3 velocity[N];
    int mass[N];
};

使用这种数据类型的归约内核可能看起来像这样（改编自 Robert 的）。

__inline__ __device__ void SumInt3(int3 const & input1, int3 const & input2, int3 & result)
{
    result.x = input1.x + input2.x;
    result.y = input1.y + input2.y;
    result.z = input1.z + input2.z;
}

__inline__ __device__ void WarpReduceInt3(int3 const & input, int3 & output, unsigned int const tID)
{
    output.x = WarpReduce(input.x, tID);
    output.y = WarpReduce(input.y, tID);
    output.z = WarpReduce(input.z, tID);    
}

__global__ void reduceSum(int3 * inputData, int3 * output, int size){
    extern __shared__ int3 sdata[];

    int3 temp;

    unsigned int tID = threadIdx.x;
    unsigned int i = tID + blockIdx.x * (blockDim.x * 2);

    // Load and sum two integer triplets, store the answer in temp.
    SumInt3(input[i], input[i + blockDim.x], temp);

    // Write the temporary answer to shared memory.
    sData[tID] = temp;

    __syncthreads();

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

    // Sum the intermediate results accross a warp.
    // No need to write the answer to shared memory,
    // as only the contribution from tID == 0 will matter.
    if (tID < 32)
    {
        WarpReduceInt3(sdata[tID], tID, temp);
    }

    if (tID == 0)
    {
        output[blockIdx.x] = temp;
    }
}