cuda - 在 Cuda 中实现 Max Reduce

Question

我一直在学习 Cuda，我仍然在掌握并行性。我目前遇到的问题是对值数组实现最大减少。这是我的内核

__global__ void max_reduce(const float* const d_array,
                     float* d_max,
                     const size_t elements)
{
    extern __shared__ float shared[];

    int tid = threadIdx.x;
    int gid = (blockDim.x * blockIdx.x) + tid;

    if (gid < elements)
        shared[tid] = d_array[gid];
    __syncthreads();

    for (unsigned int s=blockDim.x/2; s>0; s>>=1) 
    {
        if (tid < s && gid < elements)
            shared[tid] = max(shared[tid], shared[tid + s]);
        __syncthreads();
    }

    if (gid == 0)
        *d_max = shared[tid];
}

我已经使用相同的方法(用 min 替换 max 函数)实现了 min reduce，效果很好。

为了测试内核，我使用串行 for 循环找到了最小值和最大值。 min 和 max 值在内核中总是相同的，但只有 min reduce 匹配。

有什么明显的我遗漏/做错了吗？

Answer 1

您在删除的答案中的主要结论是正确的:您发布的内核没有理解这样一个事实，即在内核执行结束时，您已经完成了大量的整体减少，但结果并不完整。每个块的结果必须组合(以某种方式)。正如评论中所指出的，您的代码还有一些其他问题。让我们来看看它的修改版本:

__device__ float atomicMaxf(float* address, float val)
{
    int *address_as_int =(int*)address;
    int old = *address_as_int, assumed;
    while (val > __int_as_float(old)) {
        assumed = old;
        old = atomicCAS(address_as_int, assumed,
                        __float_as_int(val));
        }
    return __int_as_float(old);
}


__global__ void max_reduce(const float* const d_array, float* d_max, 
                                              const size_t elements)
{
    extern __shared__ float shared[];

    int tid = threadIdx.x;
    int gid = (blockDim.x * blockIdx.x) + tid;
    shared[tid] = -FLOAT_MAX;  // 1

    if (gid < elements)
        shared[tid] = d_array[gid];
    __syncthreads();

    for (unsigned int s=blockDim.x/2; s>0; s>>=1) 
    {
        if (tid < s && gid < elements)
            shared[tid] = max(shared[tid], shared[tid + s]);  // 2
        __syncthreads();
    }
    // what to do now?
    // option 1: save block result and launch another kernel
    if (tid == 0)        
        d_max[blockIdx.x] = shared[tid]; // 3
    // option 2: use atomics
    if (tid == 0)
      atomicMaxf(d_max, shared[0]);
}

正如 Pavan 所指出的，您需要初始化您的共享内存阵列。如果 gridDim.x*blockDim.x，启动的最后一个块可能不是“完整”块大于 elements .

请注意，在这一行中，即使我们正在检查正在运行的线程( gid )是否小于 elements ，当我们添加 s至 gid为了索引到共享内存，我们仍然可以索引复制到共享内存中的合法值之外的最后一个块。因此，我们需要注释 1 中指示的共享内存初始化。

正如您已经发现的那样，您的最后一行不正确。每个块产生它自己的结果，我们必须以某种方式组合它们。如果启动的块数量很少(稍后会详细介绍)，您可能会考虑的一种方法是使用 atomics .通常我们会引导人们远离使用原子，因为它们在执行时间方面是“昂贵的”。然而，我们面临的另一个选择是将块结果保存在全局内存中，完成内核，然后可能启动另一个内核来组合单个块结果。如果我最初启动了大量块(比如超过 1024 个)，那么如果我遵循这种方法，我最终可能会启动两个额外的内核。因此考虑原子。如上所述，没有本地 atomicMax浮点数函数，但如 the documentation 所示，您可以使用 atomicCAS生成任意原子函数，我在 atomicMaxf 中提供了一个例子。它为 float 提供了一个原子最大值.

但是运行 1024 个或更多原子函数(每个块一个)是最好的方法吗？可能不是。

当启动线程块的内核时，我们真的只需要启动足够的线程块来保持机器忙碌。根据经验，我们希望每个 SM 至少运行 4-8 个经线，多一些可能是个好主意。但是从机器利用率的角度来看，最初启动数千个线程块并没有什么特别的好处。如果我们选择一个数字，例如每个 SM 8 个线程块，并且我们的 GPU 中最多有 14-16 个 SM，这给了我们相对较少的 8*14 = 112 个线程块。让我们选择 128 (8*16) 作为一个不错的整数。这没有什么神奇之处，它足以让 GPU 保持忙碌。如果我们让这 128 个线程块中的每一个都做额外的工作来解决整个问题，那么我们就可以利用我们对原子的使用，而无需(也许)为此付出太多代价，并避免多次内核启动。那么这看起来如何？:

__device__ float atomicMaxf(float* address, float val)
{
    int *address_as_int =(int*)address;
    int old = *address_as_int, assumed;
    while (val > __int_as_float(old)) {
        assumed = old;
        old = atomicCAS(address_as_int, assumed,
                        __float_as_int(val));
        }
    return __int_as_float(old);
}


__global__ void max_reduce(const float* const d_array, float* d_max, 
                                              const size_t elements)
{
    extern __shared__ float shared[];

    int tid = threadIdx.x;
    int gid = (blockDim.x * blockIdx.x) + tid;
    shared[tid] = -FLOAT_MAX; 

    while (gid < elements) {
        shared[tid] = max(shared[tid], d_array[gid]);
        gid += gridDim.x*blockDim.x;
        }
    __syncthreads();
    gid = (blockDim.x * blockIdx.x) + tid;  // 1
    for (unsigned int s=blockDim.x/2; s>0; s>>=1) 
    {
        if (tid < s && gid < elements)
            shared[tid] = max(shared[tid], shared[tid + s]);
        __syncthreads();
    }

    if (tid == 0)
      atomicMaxf(d_max, shared[0]);
}

使用此修改后的内核，在创建内核启动时，我们不会根据总体数据大小 ( elements ) 决定启动多少线程块。相反，我们启动了固定数量的块(比如 128，你可以修改这个数字来找出运行最快的块)，并让每个线程块(以及整个网格)在内存中循环，计算每个元素的部分最大操作共享内存。然后，在标有注释 1 的行中，我们必须重新设置 gid变量为它的初始值。这实际上是不必要的，如果我们保证网格的大小( gridDim.x*blockDim.x )小于 elements ，则可以进一步简化块减少循环代码。 ，这在内核启动时不难做到。

请注意，在使用这种原子方法时，需要将结果(在本例中为 *d_max)初始化为适当的值，例如 -FLOAT_MAX .

同样，我们通常会引导人们远离原子使用，但在这种情况下，如果我们仔细管理它是值得考虑的，它允许我们节省额外内核启动的开销。

有关如何进行快速并行缩减的忍者级分析，请查看 Mark Harris 的优秀白皮书，该白皮书可通过相关 CUDA sample 获得。 .