CUDA - 为啥基于扭曲的并行减少速度较慢？

Posted 2023-05-08

【中文标题】CUDA - 为啥基于扭曲的并行减少速度较慢？

【英文标题】：CUDA - why is warp based parallel reduction slower?CUDA - 为什么基于扭曲的并行减少速度较慢？

【发布时间】：2012-09-25 19:58:56

【问题描述】：

我有一个基于 warp 的并行减少的想法，因为根据定义，warp 的所有线程都是同步的。

所以想法是输入数据可以减少 64 倍（每个线程减少两个元素），而无需任何同步。

与 Mark Harris 的原始实现相同，减少应用于块级，数据位于共享内存上。 http://gpgpu.org/static/sc2007/SC07_CUDA_5_Optimization_Harris.pdf

我创建了一个内核来测试他的版本和我的基于 warp 的版本。 内核本身完全相同地将 BLOCK_SIZE 元素存储在共享内存中，并在输出数组中的唯一块索引处输出其结果。

算法本身运行良好。用一个完整的数组进行测试以测试“计数”。

实现的函数体：

/**
 * Performs a parallel reduction with operator add 
 * on the given array and writes the result with the thread 0
 * to the given target value
 *
 * @param inValues T* Input float array, length must be a multiple of 2 and equal to blockDim.x
 * @param targetValue float 
 */
__device__ void reductionAddBlockThread_f(float* inValues,
    float &outTargetVar)

    // code of the below functions

1.他的版本的实现：

if (blockDim.x >= 1024 && threadIdx.x < 512)
    inValues[threadIdx.x] += inValues[threadIdx.x + 512];
__syncthreads();
if (blockDim.x >= 512 && threadIdx.x < 256)
    inValues[threadIdx.x] += inValues[threadIdx.x + 256];
__syncthreads();
if (blockDim.x >= 256 && threadIdx.x < 128)
    inValues[threadIdx.x] += inValues[threadIdx.x + 128];
__syncthreads();
if (blockDim.x >= 128 && threadIdx.x < 64)
    inValues[threadIdx.x] += inValues[threadIdx.x + 64];
__syncthreads();

//unroll last warp no sync needed
if (threadIdx.x < 32)

    if (blockDim.x >= 64) inValues[threadIdx.x] += inValues[threadIdx.x + 32];
    if (blockDim.x >= 32) inValues[threadIdx.x] += inValues[threadIdx.x + 16];
    if (blockDim.x >= 16) inValues[threadIdx.x] += inValues[threadIdx.x + 8];
    if (blockDim.x >= 8) inValues[threadIdx.x] += inValues[threadIdx.x + 4];
    if (blockDim.x >= 4) inValues[threadIdx.x] += inValues[threadIdx.x + 2];
    if (blockDim.x >= 2) inValues[threadIdx.x] += inValues[threadIdx.x + 1];

    //set final value
    if (threadIdx.x == 0)
        outTargetVar = inValues[0];

资源：

使用了 4 个同步线程 12个if语句使用 11个读+加+写操作 1 次最终写入操作 5 寄存器用法

性能：

5 次测试运行平均：~ 19.54 毫秒

2。基于 Warp 的方法：（与上述相同的函数体）

/*
 * Perform first warp based reduction by factor of 64
 *
 * 32 Threads per Warp -> LOG2(32) = 5
 *
 * 1024 Threads / 32 Threads per Warp = 32 warps
 * 2 elements compared per thread -> 32 * 2 = 64 elements per warp
 *
 * 1024 Threads/elements divided by 64 = 16
 * 
 * Only half the warps/threads are active
 */
if (threadIdx.x < blockDim.x >> 1)

    const unsigned int warpId = threadIdx.x >> 5;
    // alternative threadIdx.x & 31
    const unsigned int threadWarpId = threadIdx.x - (warpId << 5);
    const unsigned int threadWarpOffset = (warpId << 6) + threadWarpId;

    inValues[threadWarpOffset] += inValues[threadWarpOffset + 32];
    inValues[threadWarpOffset] += inValues[threadWarpOffset + 16];
    inValues[threadWarpOffset] += inValues[threadWarpOffset + 8];
    inValues[threadWarpOffset] += inValues[threadWarpOffset + 4];
    inValues[threadWarpOffset] += inValues[threadWarpOffset + 2];
    inValues[threadWarpOffset] += inValues[threadWarpOffset + 1];


// synchronize all warps - the local warp result is stored
// at the index of the warp equals the first thread of the warp
__syncthreads();

// use first warp to reduce the 16 warp results to the final one
if (threadIdx.x < 8)

    // get first element of a warp
    const unsigned int warpIdx = threadIdx.x << 6;

    if (blockDim.x >= 1024) inValues[warpIdx] += inValues[warpIdx + 512];
    if (blockDim.x >= 512) inValues[warpIdx] += inValues[warpIdx + 256];
    if (blockDim.x >= 256) inValues[warpIdx] += inValues[warpIdx + 128];
    if (blockDim.x >= 128) inValues[warpIdx] += inValues[warpIdx + 64];

    //set final value
    if (threadIdx.x == 0)
        outTargetVar = inValues[0];

资源：

使用了 1 个同步线程 7 if 语句 10个读写操作 1 次最终写入操作 5 寄存器用法

5 位移位 1 个添加 1个子

性能：

5 次测试运行平均：~ 20.82 毫秒

在具有 256 mb 浮点值的 Geforce 8800 GT 512 mb 上多次测试两个内核。 并以 每块 256 个线程（100 % 占用率）运行内核。

基于 warp 的版本慢了 ~ 1.28 毫秒。

如果未来的卡允许更大的块大小，基于 warp 的方法仍然不需要进一步的同步语句，因为最大值是 4096，它会减少到 64，而最终的 warp 会减少到 1

为什么不快？或者说这个想法，内核的缺陷在哪里？

从资源使用情况来看，warp 方法应该领先吗？

Edit1：更正内核只有一半线程处于活动状态，不会导致超出范围的读取，添加了新的性能数据

【参考方案1】：

我认为您的代码比我的代码慢的原因是在我的代码中，在第一阶段，每个 ADD 的活动扭曲数是其一半。在您的代码中，所有经线在所有第一阶段都处于活动状态。因此，总体而言，您的代码会执行更多的扭曲指令。在 CUDA 中，重要的是要考虑执行的总“warp 指令”，而不仅仅是一个 warp 执行的指令数。

另外，只使用一半的经线是没有意义的。启动 warp 只是为了让它们评估两个分支并退出会产生开销。

另一个想法是使用unsigned char 和short 实际上可能会降低您的性能。我不确定，但它肯定不会为您节省寄存器，因为它们没有打包到单个 32 位变量中。

另外，在我的原始代码中，我用模板参数 BLOCKDIM 替换了 blockDim.x，这意味着它只使用了 5 个运行时 if 语句（第二阶段的 if 被编译器消除）。

顺便说一句，计算threadWarpId 的一种更便宜的方法是

const int threadWarpId = threadIdx.x & 31;

您可以查看this article 了解更多想法。

编辑： 这是另一种基于扭曲的块减少。

template <typename T, int level>
__device__
void sumReduceWarp(volatile T *sdata, const unsigned int tid)

  T t = sdata[tid];
  if (level > 5) sdata[tid] = t = t + sdata[tid + 32];
  if (level > 4) sdata[tid] = t = t + sdata[tid + 16];
  if (level > 3) sdata[tid] = t = t + sdata[tid +  8];
  if (level > 2) sdata[tid] = t = t + sdata[tid +  4];
  if (level > 1) sdata[tid] = t = t + sdata[tid +  2];
  if (level > 0) sdata[tid] = t = t + sdata[tid +  1];


template <typename T>
__device__
void sumReduceBlock(T *output, volatile T *sdata)

  // sdata is a shared array of length 2 * blockDim.x

  const unsigned int warp = threadIdx.x >> 5;
  const unsigned int lane = threadIdx.x & 31;
  const unsigned int tid  = (warp << 6) + lane;

  sumReduceWarp<T, 5>(sdata, tid);
  __syncthreads();

  // lane 0 of each warp now contains the sum of two warp's values
  if (lane == 0) sdata[warp] = sdata[tid];

  __syncthreads();

  if (warp == 0) 
    sumReduceWarp<T, 4>(sdata, threadIdx.x);
    if (lane == 0) *output = sdata[0];
  

这应该会快一点，因为它使用了在第一阶段启动的所有经线，并且在最后阶段没有分支，代价是额外的分支、共享加载/存储和__syncthreads()新的中间阶段。我没有测试过这段代码。如果你运行它，请告诉我它的性能。如果您在原始代码中为 blockDim 使用模板，它可能会更快，但我认为这段代码更简洁。

注意使用临时变量t是因为Fermi和后来的架构使用纯加载/存储架构，所以+=从共享内存到共享内存会导致额外的负载（因为sdata指针必须是易失的）。一次显式加载到临时文件中可以避免这种情况。在 G80 上，它不会对性能产生影响。

【讨论】：

我必须道歉说共享内存数组长度等于 BLOCK_SIZE，但它必须等于双块（这在从全局内存初始化时没有意义）大小或只有一半线程应该是活动的。而且由于所有线程都处于活动状态，因此存在超出范围的读写。我将更正此问题并再次检查。修正后的内核现在绝对是正确的。正如我所说，它总是产生正确的数学结果。奇怪的是为什么这些共享内存的越界访问没有产生错误的结果。

关于 BLOCKDIM 模板参数：这就是为什么我发布了我的算法实现。而且我仍然需要更昂贵的方法来计算warpId，因为我需要它来计算threadWarpOffset。

抱歉，我错过了第三行中的 6 而不是 5。已编辑。我认为我对为什么您的代码较慢的推理可能至少是部分原因。

从所有线程参与的第一个错误内核开始，total warp instructions 明显更高。但是，如果只有一半线程处于活动状态，而其余线程正在等待下一个同步语句，为什么这会相关呢？所以他们也可以做点什么，而不是等待？特别是如果经线中的线程是同步的？

空闲扭曲不使用处理周期。主动经线可以。您的代码总体上具有更多的活动扭曲以获取更多指令，因此总周期更高。我的代码在更少的总周期内执行相同的计算。我认为经线减少是可行的（我自己做过），但我会以不同的方式做。我会尝试挖掘一些代码...

【参考方案2】：

您还应该查看 SDK 中的示例。我记得一个很好的例子，它实现了几种减少方法。其中至少有一个也使用基于扭曲的减少。

（我现在无法查找名称，因为我只在我的另一台机器上安装了它）