CUDA Reduce：从基础写法到工程优化

Reduce 的目标，是把一段数据逐步归约成一个标量。下面按版本梳理 CUDA 中常见的优化路径，从最基础的写法一路走到工程里更常见的实现。

默认约定：

BLOCK_SIZE 是 2 的幂，工程里通常也会取成 32 的倍数
每个 block 先算出一个部分和，再写到 g_odata[blockIdx.x]
如果还有多个 block 的结果，就继续 launch 下一轮 kernel，直到只剩 1 个值

优化思想汇总：gm加载→共享内存→warp分化→最后几次同步等待太严重→shuffle

v0：最基础的版本

这个版本直接在 global memory 上做归约，主要用来理解整体流程。

__global__ void reduce_v0(float* g_idata, float* g_odata) {
    unsigned int tid = threadIdx.x;
    unsigned int idx = blockIdx.x * blockDim.x + tid;
 
    for (unsigned int s = 1; s < blockDim.x; s *= 2) {
        if (tid % (2 * s) == 0) {
            g_idata[idx] += g_idata[idx + s];
        }
        __syncthreads();
    }
 
    if (tid == 0) {
        g_odata[blockIdx.x] = g_idata[blockIdx.x * blockDim.x];
    }
}

特点：

写法直观，容易看清楚归约过程。
s = 1, 2, 4, 8, ... 的推进方式很容易理解。

问题：

直接读写 global memory，开销大。
tid % (2 * s) == 0 会导致明显的 warp divergence。
线程利用率低。

v1：Shared Memory + 交错归约

思路：

每个线程先读 1 个元素到 shared memory。
再按 s = 1, 2, 4, 8, ... 做归约。
用 tid % (2 * s) == 0 决定谁参与当前轮次。

__global__ void reduce_v1(float* g_idata, float* g_odata) {
    __shared__ float sdata[BLOCK_SIZE];
 
    unsigned int tid = threadIdx.x;
    unsigned int i   = blockIdx.x * blockDim.x + tid;
 
    sdata[tid] = g_idata[i];
    __syncthreads();
 
    for (unsigned int s = 1; s < blockDim.x; s *= 2) {
        // s = 1 时，参与的是 tid = 0, 2, 4, 6, ...
        // tid = 0: sdata[0] += sdata[1]
        // tid = 2: sdata[2] += sdata[3]
        // tid = 4: sdata[4] += sdata[5]
        // s = 2 时，参与的是 tid = 0, 4, 8, 12, ...
        // tid = 0: sdata[0] += sdata[2]
        // tid = 4: sdata[4] += sdata[6]
        if (tid % (2 * s) == 0) {
            sdata[tid] += sdata[tid + s];
        }
        __syncthreads();
    }
 
    if (tid == 0) {
        g_odata[blockIdx.x] = sdata[0];
    }
}

问题：

线程闲置多。
warp 分化严重。
访问模式不够友好。

v2：每个线程先处理 2 个元素

目的：

让每个线程在加载阶段先多做一点事。
减少一部分线程空转。

这一版主要优化的是加载阶段，归约阶段还是 v1 的写法。

__global__ void reduce_v2(float* g_idata, float* g_odata) {
    __shared__ float sdata[BLOCK_SIZE];
 
    unsigned int tid = threadIdx.x;
    // 一个 block 覆盖 blockDim.x * 2 个元素
    unsigned int i   = blockIdx.x * (blockDim.x * 2) + tid;
 
    sdata[tid] = g_idata[i] + g_idata[i + blockDim.x];
    __syncthreads();
 
    for (unsigned int s = 1; s < blockDim.x; s *= 2) {
        if (tid % (2 * s) == 0) {
            sdata[tid] += sdata[tid + s];
        }
        __syncthreads();
    }
 
    if (tid == 0) {
        g_odata[blockIdx.x] = sdata[0];
    }
}

收益：

每个线程一开始就处理 2 个元素。
block 数量会相应减少一些。

遗留问题：

归约阶段依然存在较强的 warp divergence。

v3：顺序归约

核心变化：

不再让“隔一个线程工作一次”。
改成“前半部分线程工作，后半部分线程休息”。
每一轮都让活跃线程连续分布，减少分支分化。

__global__ void reduce_v3(float* g_idata, float* g_odata) {
    __shared__ float sdata[BLOCK_SIZE];
 
    unsigned int tid = threadIdx.x;
    unsigned int i   = blockIdx.x * (blockDim.x * 2) + tid;
 
    sdata[tid] = g_idata[i] + g_idata[i + blockDim.x];
    __syncthreads();
 
    // 每轮把参与归约的范围缩小一半
    for (unsigned int s = blockDim.x / 2; s > 0; s >>= 1) {
        if (tid < s) {
            sdata[tid] += sdata[tid + s];
        }
        __syncthreads();
    }
 
    if (tid == 0) {
        g_odata[blockIdx.x] = sdata[0];
    }
}

收益：

活跃线程集中在前半段，warp divergence 明显减轻。
shared memory 的访问模式也更整齐。

v4：展开最后一个 Warp

思路：

当 block 内只剩最后一个 warp 时，不再每轮都做 __syncthreads()。
直接把最后几步手动展开，减少同步开销。

__device__ void warpReduce(volatile float* sdata, unsigned int tid) {
    sdata[tid] += sdata[tid + 32];
    sdata[tid] += sdata[tid + 16];
    sdata[tid] += sdata[tid + 8];
    sdata[tid] += sdata[tid + 4];
    sdata[tid] += sdata[tid + 2];
    sdata[tid] += sdata[tid + 1];
}
 
__global__ void reduce_v4(float* g_idata, float* g_odata) {
    __shared__ float sdata[BLOCK_SIZE];
 
    unsigned int tid = threadIdx.x;
    unsigned int i   = blockIdx.x * (blockDim.x * 2) + tid;
 
    sdata[tid] = g_idata[i] + g_idata[i + blockDim.x];
    __syncthreads();
 
    // 归约到只剩 64 个值时停下，最后一段交给单个 warp 处理
    for (unsigned int s = blockDim.x / 2; s > 32; s >>= 1) {
        if (tid < s) {
            sdata[tid] += sdata[tid + s];
        }
        __syncthreads();
    }
 
    if (tid < 32) {
        warpReduce(sdata, tid);
    }
 
    if (tid == 0) {
        g_odata[blockIdx.x] = sdata[0];
    }
}

收益：

避免最后几轮重复同步。
进一步减少归约尾部的固定开销。

v5：用 Shuffle 做最后一个 Warp

思路：

v4 的尾归约仍然依赖 shared memory。
可以先把前 32 个线程与后 32 个线程的结果合并，再用 __shfl_down_sync 在一个 warp 内完成剩余归约。
这也是后面工程写法的一个过渡版本。

__device__ __forceinline__ float warpReduceShuffle(float val) {
    constexpr unsigned int FULL_MASK = 0xffffffff;
    for (int offset = 16; offset > 0; offset >>= 1) {
        val += __shfl_down_sync(FULL_MASK, val, offset);
    }
    return val;
}
 
__global__ void reduce_v5(float* g_idata, float* g_odata) {
    __shared__ float sdata[BLOCK_SIZE];
 
    unsigned int tid = threadIdx.x;
    unsigned int i   = blockIdx.x * (blockDim.x * 2) + tid;
 
    sdata[tid] = g_idata[i] + g_idata[i + blockDim.x];
    __syncthreads();
 
    for (unsigned int s = blockDim.x / 2; s > 32; s >>= 1) {
        if (tid < s) {
            sdata[tid] += sdata[tid + s];
        }
        __syncthreads();
    }
 
    if (tid < 32) {
        float val = sdata[tid] + sdata[tid + 32];
        val = warpReduceShuffle(val);
        if (tid == 0) {
            g_odata[blockIdx.x] = val;
        }
    }
}

收益：

最后一个 warp 的归约不再依赖 volatile shared memory。
shared memory 读写更少，写法也更贴近现代 CUDA 的常见实现。

工程里更常见的写法

实际工程里，常见组合通常是：

grid-stride loop
寄存器累加
warp reduce
少量 shared memory 做跨 warp 合并

__device__ __forceinline__ float warp_reduce_sum(float val) {
    constexpr unsigned int FULL_MASK = 0xffffffff;
    for (int offset = 16; offset > 0; offset >>= 1) {
        val += __shfl_down_sync(FULL_MASK, val, offset);
    }
    return val;
}
 
__global__ void reduce_better(const float* __restrict__ g_idata,
                              float* __restrict__ g_odata,
                              int N) {
    __shared__ float sdata[BLOCK_SIZE / 32]; // 每个 warp 一个槽
 
    unsigned int tid = threadIdx.x;
    unsigned int idx = blockIdx.x * blockDim.x * 2 + tid;
    unsigned int stride = blockDim.x * 2 * gridDim.x;
 
    float sum = 0.0f;
 
    // 1. 每个线程通过 grid-stride loop 处理多组数据
    while (idx < N) {
        sum += g_idata[idx];
        if (idx + blockDim.x < N) {
            sum += g_idata[idx + blockDim.x];
        }
        idx += stride;
    }
 
    // 2. 先在 warp 内做一次归约
    sum = warp_reduce_sum(sum);
 
    if ((tid & 31) == 0) {
        sdata[tid >> 5] = sum;
    }
    __syncthreads();
 
    // 3. 第一个 warp 继续把所有 warp 的结果归约掉
    sum = (tid < blockDim.x / 32) ? sdata[tid] : 0.0f;
    if (tid < 32) {
        sum = warp_reduce_sum(sum);
    }
 
    if (tid == 0) {
        g_odata[blockIdx.x] = sum;
    }
}

对应的 kernel launch 可以写成：

constexpr int BLOCK_SIZE = 256;
 
while (n > 1) {
    int blocks = (n + BLOCK_SIZE * 2 - 1) / (BLOCK_SIZE * 2);
    reduce_better<<<blocks, BLOCK_SIZE>>>(src, dst, n);
    n = blocks;
    std::swap(src, dst);
}

循环结束后，结果就在 src[0]

fxlin's Blog

目录

CUDA Reduce