CUDA 减少了许多大小不等的小型阵列

Question

我想知道是否有人可以建议在 CUDA 中计算大量相对较小但不同大小的数组的均值/标准差的最佳方法？

SDK 中的并行缩减示例适用于单个非常大的数组，看起来大小方便地是每个块的线程数的倍数，但我的情况完全不同：

然而，从概念上讲，我有大量的对象，每个对象都包含两个组件，upper 和 lower，每个组件都有一个 x 和一个 y 坐标。即

upper.x, lower.x, upper.y, lower.y

每个数组的长度约为800，但它在对象之间（不在对象内）有所不同，例如

Object1.lower.x = 1.1, 2.2, 3.3
Object1.lower.y = 4.4, 5.5, 6.6
Object1.upper.x = 7.7, 8.8, 9.9
Object1.upper.y = 1.1, 2.2, 3.3

Object2.lower.x = 1.0,  2.0,  3.0,  4.0, 5.0 
Object2.lower.y = 6.0,  7.0,  8.0,  9.0, 10.0
Object2.upper.x = 11.0, 12.0, 13.0, 14.0, 15.0 
Object2.upper.y = 16.0, 17.0, 18.0, 19.0, 20.0

请注意，以上只是我表示数组的方式，我的数据没有存储在C 结构或类似的东西中：数据可以按我需要的任何方式组织。关键是，对于每个数组，都需要计算平均值、标准偏差和最终的直方图，并且在一个特定对象内，需要计算数组之间的比率和差异。

我应该如何将这些数据发送到 GPU 设备并组织我的线程块层次结构？我的一个想法是对我的所有数组进行零填充，以使它们具有相同的长度，并且有一组块在每个对象上工作，但该方法似乎存在各种问题，如果它可以工作的话。

提前致谢

Answer 1

如果我理解正确，您想将 Object1.lower.x 减少为一个结果，将 Object1.lower.y 减少为另一个结果，依此类推。对于任何给定的对象，有四个要减少的数组，所有数组的长度都相等（对于对象）。

对此有许多可能的方法，其中一个影响因素是系统中的对象总数。我假设这个数字很大。

为了获得最佳性能，您需要一种最佳的内存访问模式，并且您希望避免分歧。由于全等数组的数量是 4，如果您采取简单的方法，即在下面每个线程执行一个数组，那么您不仅会遭受内存访问不佳的困扰，而且硬件还需要检查每次迭代中的线程warp 需要执行循环 - 那些不需要执行的循环将被禁用，这可能效率低下（尤其是如果一个数组比其他数组长得多，例如）。

for (int i = 0 ; i < myarraylength ; i++)
    sum += myarray[i];

相反，如果您让每个 warp 对一个数组求和，那么不仅效率更高，而且您的内存访问模式也会更好，因为相邻线程将读取相邻元素 [1]。

for (int i = tidwithinwarp ; i < warparraylength ; i += warpsize)
{
    mysum += warparray[i];
}
mysum = warpreduce(mysum);

您还应该考虑数组的对齐方式，最好在 64 字节边界上对齐，尽管如果您正在开发 1.2 或更高的计算能力，那么这并不像在旧 GPU 上那么重要。

在此示例中，您将启动每个块的四个 warp，即 128 个线程，以及与对象一样多的块。

[1] 你确实说过你可以选择你喜欢的任何内存排列，通常交错数组很有用，这样 array[0][0] 就在 array[1][0] 旁边，因为这意味着相邻线程可以对相邻数组进行操作并获得合并访问。然而，由于数组的长度不是恒定的，这可能很复杂，需要填充较短的数组。

Answer 2

作为汤姆回答的后续，我想提一下 warp reduction 可以通过CUB 轻松实现。

这是一个有效的例子：

#include <cub/cub.cuh>
#include <cuda.h>

#include "Utilities.cuh"

#include <iostream>

#define WARPSIZE    32
#define BLOCKSIZE   256

const int N = 1024;

/*************************/
/* WARP REDUCTION KERNEL */
/*************************/
__global__ void sum(const float * __restrict__ indata, float * __restrict__ outdata) {

    unsigned int tid = blockIdx.x * blockDim.x + threadIdx.x;

    unsigned int warp_id = threadIdx.x / WARPSIZE;

    // --- Specialize WarpReduce for type float. 
    typedef cub::WarpReduce<float, WARPSIZE> WarpReduce;

    // --- Allocate WarpReduce shared memory for (N / WARPSIZE) warps
    __shared__ typename WarpReduce::TempStorage temp_storage[BLOCKSIZE / WARPSIZE];

    float result;
    if(tid < N) result = WarpReduce(temp_storage[warp_id]).Sum(indata[tid]);

    if(tid % WARPSIZE == 0) outdata[tid / WARPSIZE] = result;
}

/********/
/* MAIN */
/********/
int main() {

    // --- Allocate host side space for 
    float *h_data       = (float *)malloc(N * sizeof(float));
    float *h_result     = (float *)malloc((N / WARPSIZE) * sizeof(float));

    float *d_data;      gpuErrchk(cudaMalloc(&d_data, N * sizeof(float)));
    float *d_result;    gpuErrchk(cudaMalloc(&d_result, (N / WARPSIZE) * sizeof(float)));

    for (int i = 0; i < N; i++) h_data[i] = (float)i;

    gpuErrchk(cudaMemcpy(d_data, h_data, N * sizeof(float), cudaMemcpyHostToDevice));

    sum<<<iDivUp(N, BLOCKSIZE), BLOCKSIZE>>>(d_data, d_result);
    gpuErrchk(cudaPeekAtLastError());
    gpuErrchk(cudaDeviceSynchronize());

    gpuErrchk(cudaMemcpy(h_result, d_result, (N / WARPSIZE) * sizeof(float), cudaMemcpyDeviceToHost));

    std::cout << "output: ";
    for(int i = 0; i < (N / WARPSIZE); i++) std::cout << h_result[i] << " ";
    std::cout << std::endl;

    gpuErrchk(cudaFree(d_data));
    gpuErrchk(cudaFree(d_result));

    return 0;
}

在本例中，创建了一个长度为N 的数组，结果是32 连续元素的总和。所以

result[0] = data[0] + ... + data[31];
result[1] = data[32] + ... + data[63];
....