在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 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]);
}

1. 正如Pavan所说,您需要初始化共享内存阵列.如果gridDim.x*blockDim.x大于,则启动的最后一个块可能不是"完整"块elements.
2. 请注意,在这一行,即使我们检查线程操作(gid)小于elements,当我们添加s到gid索引到共享内存中,我们依然可以复制到共享内存中的最后一块合法值以外的指标.因此,我们需要注释1中指示的共享内存初始化.
3. 正如您已经发现的那样,您的最后一行不正确.每个块产生它自己的结果,我们必须以某种方式组合它们.如果启动的块数很少(稍后会详细介绍),可以考虑使用原子的一种方法.通常我们引导人们远离使用原子,因为它们在执行时间方面"代价高昂".但是,我们面临的另一个选择是将块结果保存在全局内存中,完成内核,然后可能启动另一个内核来组合各个块结果.如果我最初启动了大量的块(例如超过1024个),那么如果我遵循这种方法,我可能最终会启动两个额外的内核.因此考虑原子.如上所述,没有本地人atomicMax浮点数的函数,但是如文档中所示,您可以使用atomicCAS生成任意原子函数,并且我提供了一个示例,atomicMaxf其中提供了原子最大值float.

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

在启动线程块的内核时,我们实际上只需要启动足够的线程块来保持机器忙碌.根据经验,我们希望每个SM至少运行4-8个warp,而更多可能是一个好主意.但是,从最初启动数千个线程块的机器利用率角度来看,并没有特别的好处.如果我们选择一个像每个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变量重新设置为它的初始值.这实际上是不必要的,如果我们保证grid(gridDim.x*blockDim.x)的大小小于(elements在内核启动时不难),则可以进一步简化块减少循环代码.

请注意,使用此原子方法时,必须将结果(*d_max在本例中)初始化为适当的值,例如-FLOAT_MAX.

同样,我们通常会引导人们从原子使用方式出发,但在这种情况下,如果我们仔细管理它,则值得考虑,它允许我们节省额外内核启动的开销.

有关如何进行快速平行缩减的忍者级分析,请查看Mark Harris的优秀白皮书,该白皮书可与相关的CUDA样本一起使用.