前言

第十章学了归约——把 N 个数变成 1 个数。这章学前缀和（Prefix Sum），也叫 Scan——把 N 个数变成 N 个数，每个位置存储之前所有元素的累计结果。看起来计算量更大，但前缀和是并行计算的"瑞士军刀"，能解决很多看似不可并行的问题：流压缩、基数排序、稀疏矩阵运算等。第十一章重点讲解工作效率——如何让并行算法做的总工作量不超过串行算法太多。

📦 配套资源：本系列文章配有完整的 GitHub 仓库，包含每章的练习题解答、CUDA 代码实现和详细注释。所有代码都经过测试，可以直接运行。

前缀和基础

什么是前缀和

前缀和：对于输入数组 [a₀, a₁, a₂, ..., aₙ₋₁]，输出每个位置之前所有元素的累计值。

有两种变体：

包含式扫描（Inclusive Scan）：包含当前元素

1
2
3
输入:  [3, 1, 7, 0, 4, 1, 6, 3]
输出:  [3, 4, 11, 11, 15, 16, 22, 25]
       3  3+1  4+7  11+0  ...

排除式扫描（Exclusive Scan）：不包含当前元素

1
2
3
输入:  [3, 1, 7, 0, 4, 1, 6, 3]
输出:  [0, 3, 4, 11, 11, 15, 16, 22]
       0  前缀  前缀+1  前缀+7  ...

关系：exclusive[i] = inclusive[i-1]，exclusive[0] = 0

数学定义

对于二元结合运算 ⊕：

Inclusive Scan:

yi=⨁j=0ixj=x0⊕x1⊕...⊕xiy_i = \bigoplus_{j=0}^{i} x_j = x_0 \oplus x_1 \oplus ... \oplus x_i

Exclusive Scan:

yi=⨁j=0i−1xj=x0⊕x1⊕...⊕xi−1y_i = \bigoplus_{j=0}^{i-1} x_j = x_0 \oplus x_1 \oplus ... \oplus x_{i-1}

其中 y₀ = identity（单位元）。

应用场景

前缀和用途极广：

应用 说明
流压缩 根据条件筛选元素
基数排序 计算每个桶的起始位置
稀疏矩阵 CSR 格式的行指针数组
多项式求值 Horner 法则的并行化
求解三对角方程 循环归约
直方图均衡化 累积分布函数

核心思想：前缀和把"依赖前面结果"的问题转化为可并行的形式。

串行实现

1
2
3
4
5
6
void inclusive_scan_sequential(float *x, float *y, int n) {
    y[0] = x[0];
    for (int i = 1; i < n; i++) {
        y[i] = y[i-1] + x[i];
    }
}

时间复杂度 O(n)，存在严格的数据依赖：y[i] 依赖 y[i-1]。

看似无法并行？其实可以！

Kogge-Stone 算法

核心思想

利用结合律，不需要按顺序计算：

1
2
y[i] = x[0] + x[1] + ... + x[i]
     = (x[0] + x[1] + ... + x[i-k]) + (x[i-k+1] + ... + x[i])

思路：每一步让每个元素与距离为 2^k 的元素相加。

算法过程

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
初始:     [3, 1, 7, 0, 4, 1, 6, 3]

Step 1 (stride=1):
  y[i] = x[i] + x[i-1]
  结果: [3, 4, 8, 7, 4, 5, 7, 9]

Step 2 (stride=2):
  y[i] = y[i] + y[i-2]
  结果: [3, 4, 11, 11, 12, 12, 11, 14]

Step 3 (stride=4):
  y[i] = y[i] + y[i-4]
  结果: [3, 4, 11, 11, 15, 16, 22, 25]

完成！

每一步，每个元素都并行更新。log₂N 步后完成。

CUDA 实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
__global__ void kogge_stone_scan(float *X, float *Y, int n) {
    __shared__ float XY[SECTION_SIZE];
    
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    
    // 加载到共享内存
    if (i < n) {
        XY[threadIdx.x] = X[i];
    } else {
        XY[threadIdx.x] = 0;
    }
    
    // Kogge-Stone 扫描
    for (unsigned int stride = 1; stride < blockDim.x; stride *= 2) {
        __syncthreads();
        float temp;
        if (threadIdx.x >= stride) {
            temp = XY[threadIdx.x] + XY[threadIdx.x - stride];
        }
        __syncthreads();
        if (threadIdx.x >= stride) {
            XY[threadIdx.x] = temp;
        }
    }
    
    // 写回
    if (i < n) {
        Y[i] = XY[threadIdx.x];
    }
}

为什么需要两次同步？

1
2
3
4
__syncthreads();  // 第一次：确保所有线程读完
temp = XY[...] + XY[...];  // 读取
__syncthreads();  // 第二次：确保所有线程读完再写
XY[...] = temp;  // 写入

如果不用临时变量：

1
2
// 错误！
XY[i] = XY[i] + XY[i-stride];  // 如果邻居还没读完就被覆盖

关键：读和写必须分离，用同步保证。

工作量分析

指标 串行算法 Kogge-Stone
步数 N log₂N
每步操作 1 ~N
总操作 N N·log₂N

问题：Kogge-Stone 做了更多的工作！

对于 N=1024：

  • 串行：1024 次操作
  • Kogge-Stone：1024 × 10 = 10240 次操作

这就是**工作效率（Work Efficiency）**问题。

Brent-Kung 算法

工作效率优化

思路：减少冗余计算，分两个阶段：

  1. 归约阶段（Reduce）：自底向上，构建部分和
  2. 分发阶段（Downsweep）：自顶向下，分发结果

算法过程

阶段 1：归约

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
初始:     [3, 1, 7, 0, 4, 1, 6, 3]

Step 1 (stride=1): 相邻元素配对求和
          [3, 4, 7, 7, 4, 5, 6, 9]
              ↑     ↑     ↑     ↑
             1+3   0+7   1+4   3+6

Step 2 (stride=2): 隔一个配对
          [3, 4, 7, 11, 4, 5, 6, 14]
                   ↑           ↑
                  4+7         5+9

Step 3 (stride=4): 隔三个配对
          [3, 4, 7, 11, 4, 5, 6, 25]
                               ↑
                             11+14

现在 XY[7] = 25 = 全部元素之和。

阶段 2：分发（Downsweep）

1
2
3
4
5
6
7
8
9
10
11
12
13
从右到左，把部分和"传播"下去：

Step 1 (stride=2):
  XY[5] = XY[5] + XY[3] = 5 + 11 = 16
          [3, 4, 7, 11, 4, 16, 6, 25]

Step 2 (stride=1):
  XY[2] = XY[2] + XY[1] = 7 + 4 = 11
  XY[4] = XY[4] + XY[3] = 4 + 11 = 15
  XY[6] = XY[6] + XY[5] = 6 + 16 = 22
          [3, 4, 11, 11, 15, 16, 22, 25]

完成！

CUDA 实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
__global__ void brent_kung_scan(float *X, float *Y, int n) {
    __shared__ float XY[SECTION_SIZE];
    
    int i = 2 * blockIdx.x * blockDim.x + threadIdx.x;
    
    // 加载（每线程两个元素）
    if (i < n) XY[threadIdx.x] = X[i];
    if (i + blockDim.x < n) XY[threadIdx.x + blockDim.x] = X[i + blockDim.x];
    
    // ========== 归约阶段 ==========
    for (unsigned int stride = 1; stride <= blockDim.x; stride *= 2) {
        __syncthreads();
        int index = (threadIdx.x + 1) * 2 * stride - 1;
        if (index < SECTION_SIZE) {
            XY[index] += XY[index - stride];
        }
    }
    
    // ========== 分发阶段 ==========
    for (unsigned int stride = SECTION_SIZE / 4; stride > 0; stride /= 2) {
        __syncthreads();
        int index = (threadIdx.x + 1) * 2 * stride - 1;
        if (index + stride < SECTION_SIZE) {
            XY[index + stride] += XY[index];
        }
    }
    
    __syncthreads();
    
    // 写回
    if (i < n) Y[i] = XY[threadIdx.x];
    if (i + blockDim.x < n) Y[i + blockDim.x] = XY[threadIdx.x + blockDim.x];
}

工作量分析

阶段 步数 每步操作 总操作
归约 log₂N N/2^k N-1
分发 log₂N - 1 N/2^k N-1-log₂N
总计 2N-2-log₂N

对于 N=1024：

  • Kogge-Stone：10240 次操作
  • Brent-Kung：~2046 次操作

工作效率提升 5 倍！

为什么每线程处理两个元素？

1
int i = 2 * blockIdx.x * blockDim.x + threadIdx.x;

如果 Section 有 1024 个元素，只需要 512 个线程。每线程加载两个元素，减少了线程开销。

Kogge-Stone vs Brent-Kung

对比

指标 Kogge-Stone Brent-Kung
总操作数 N·log₂N 2N - 2 - log₂N
步数 log₂N 2·log₂N - 1
并行度（每步） N 递减/递增
控制逻辑 简单 复杂
工作效率 低 高（~串行）

选择策略

场景 推荐算法 原因
小数组（<1K） Kogge-Stone 简单，步数少
大数组 Brent-Kung 工作量小
带宽受限 Kogge-Stone 更多并行隐藏延迟
计算受限 Brent-Kung 总操作少

实际中，混合策略最常见：大部分用 Brent-Kung，最后几层用 Kogge-Stone。

三阶段分层扫描

问题：超出单 Block

当数据量超过共享内存容量时，需要多 Block 协作。

三阶段算法

阶段 1：Block 内扫描

每个 Block 独立扫描自己的 Section，保存最后一个元素（Section Sum）。

1
2
3
4
Block 0: [3,4,11,11]  → sum[0] = 11
Block 1: [4,5,11,14]  → sum[1] = 14
Block 2: [5,6,12,15]  → sum[2] = 15
...

阶段 2：扫描 Section Sums

对所有 Block 的 sum 做扫描：

1
2
sum: [11, 14, 15, ...]
扫描后: [0, 11, 25, 40, ...]  // Exclusive

这给出了每个 Block 的"起始偏移"。

阶段 3：加上偏移

每个 Block 的所有元素加上对应的偏移：

1
2
3
4
Block 0: [3,4,11,11] + 0  → [3,4,11,11]
Block 1: [4,5,11,14] + 11 → [15,16,22,25]
Block 2: [5,6,12,15] + 25 → [30,31,37,40]
...

实现

Kernel 1：Block 内扫描

1
2
3
4
5
6
7
8
9
10
11
12
13
14
__global__ void scan_phase1(float *X, float *Y, float *S, int n) {
    __shared__ float XY[SECTION_SIZE];
    
    // Block 内扫描（Brent-Kung）
    // ...
    
    // 保存 Section Sum
    if (threadIdx.x == blockDim.x - 1) {
        S[blockIdx.x] = XY[SECTION_SIZE - 1];
    }
    
    // 写回
    // ...
}

Kernel 2：扫描 Section Sums

1
2
// 通常只有一个 Block，因为 Section 数量不多
scan_single_block<<<1, numSections>>>(S, S, numSections);

Kernel 3：加偏移

1
2
3
4
5
6
__global__ void scan_phase3(float *Y, float *S, int n) {
    int i = blockIdx.x * SECTION_SIZE + threadIdx.x;
    if (blockIdx.x > 0 && i < n) {
        Y[i] += S[blockIdx.x - 1];  // Exclusive 扫描的结果
    }
}

递归处理

如果 Section 数量也很大，对 S 的扫描本身需要分层。形成递归结构：

1
2
3
4
5
6
7
8
Level 0: 扫描原始数据 → Section Sums
Level 1: 扫描 Section Sums → Super Section Sums
Level 2: 扫描 Super Section Sums
...
Level K: 单 Block 完成

然后逆向加偏移：
Level K → Level K-1 → ... → Level 0

层数：log(N/SECTION_SIZE)

单遍扫描（Single-Pass Scan）

动机

三阶段扫描需要多次 Kernel 启动，有开销。能否一遍完成？

思路

用原子操作 + 内存标志实现 Block 间通信：

  1. 每个 Block 完成扫描后，发布自己的 Section Sum
  2. 等待前一个 Block 的结果
  3. 加上前缀并继续

实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
__global__ void single_pass_scan(float *X, float *Y, float *flags, 
                                  float *sums, int n) {
    __shared__ float XY[SECTION_SIZE];
    __shared__ float prefix;
    
    int bid = blockIdx.x;
    
    // Phase 1: Block 内扫描
    // ... 标准 Brent-Kung ...
    
    // Phase 2: 等待前一个 Block
    if (threadIdx.x == 0) {
        // 发布自己的和
        sums[bid] = XY[SECTION_SIZE - 1];
        __threadfence();  // 确保写入对其他 Block 可见
        atomicExch(&flags[bid], 1);  // 标记完成
        
        // 等待前缀
        if (bid > 0) {
            // 等待 Block bid-1
            while (atomicAdd(&flags[bid - 1], 0) == 0);
            __threadfence();
            
            // 累加所有前缀（可优化为 Kogge-Stone 风格）
            float pre = 0;
            for (int i = 0; i < bid; i++) {
                pre += sums[i];
            }
            prefix = pre;
        } else {
            prefix = 0;
        }
    }
    __syncthreads();
    
    // Phase 3: 加偏移并写回
    int i = bid * SECTION_SIZE + threadIdx.x;
    if (i < n) {
        Y[i] = XY[threadIdx.x] + prefix;
    }
}

问题

顺序依赖：Block 必须按顺序完成，可能导致串行化。

改进：使用 Decoupled Look-back（分离回溯）算法，每个 Block 只查看部分前缀，链式传播。

工作效率的本质

时间 vs 操作

并行时间：所有处理器同时工作需要的时间步数

总操作数：所有处理器执行的操作总和

工作效率 = 串行操作数 / 并行操作数

算法 并行时间 总操作 工作效率
串行 N N 100%
Kogge-Stone log N N·log N 1/log N
Brent-Kung 2·log N 2N ~50%

何时追求工作效率？

处理器充足时：

  • 操作数主导时间
  • 应该追求工作效率

处理器不足时：

  • 步数主导时间
  • 可以牺牲工作效率换取更少步数

GPU 通常处理器"足够多"，所以工作效率很重要。

前缀和的应用

流压缩（Stream Compaction）

筛选满足条件的元素：

1
2
3
4
5
6
输入:    [3, 1, 7, 0, 4, 1, 6, 3]
条件:    x > 2
标志:    [1, 0, 1, 0, 1, 0, 1, 1]
前缀和:  [0, 1, 1, 2, 2, 3, 3, 4]  // Exclusive
输出位置: 0, -, 1, -, 2, -, 3, 4
输出:    [3, 7, 4, 6, 3]

Exclusive 前缀和给出每个满足条件元素的输出位置！

1
2
3
4
5
6
7
__global__ void compact(int *flags, int *positions, 
                        float *in, float *out, int n) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < n && flags[i]) {
        out[positions[i]] = in[i];
    }
}

基数排序

按位排序，每一位用前缀和确定位置：

1
2
第 0 位 = 0 的元素: 前缀和确定位置
第 0 位 = 1 的元素: 紧随其后

分配工作

根据负载分配任务：

1
2
3
4
5
6
任务大小:  [5, 3, 8, 2, 4]
前缀和:    [0, 5, 8, 16, 18]

线程 7 应该处理哪个任务？
→ 二分查找前缀和数组
→ 任务 1（因为 5 ≤ 7 < 8）

CUB 库

使用 CUB 实现扫描

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
#include <cub/cub.cuh>

void scanWithCub(float *d_in, float *d_out, int n) {
    void *d_temp = nullptr;
    size_t temp_bytes = 0;
    
    // 确定临时存储大小
    cub::DeviceScan::InclusiveSum(d_temp, temp_bytes, d_in, d_out, n);
    
    cudaMalloc(&d_temp, temp_bytes);
    
    // 执行扫描
    cub::DeviceScan::InclusiveSum(d_temp, temp_bytes, d_in, d_out, n);
    
    cudaFree(d_temp);
}

CUB 提供的扫描

函数 功能
DeviceScan::InclusiveSum 包含式求和
DeviceScan::ExclusiveSum 排除式求和
DeviceScan::InclusiveScan 包含式自定义
DeviceScan::ExclusiveScan 排除式自定义

自定义操作

1
2
3
4
5
6
7
8
9
struct MaxOp {
    __device__ float operator()(float a, float b) {
        return (a > b) ? a : b;
    }
};

MaxOp max_op;
cub::DeviceScan::InclusiveScan(d_temp, temp_bytes, 
                                d_in, d_out, max_op, n);

小结

第十一章深入讲解前缀和：

两种变体：Inclusive（含当前）vs Exclusive（不含当前）。Exclusive 的应用更广，特别是确定输出位置。

Kogge-Stone：简单直接，log N 步，但工作量 N·log N，工作效率低。

Brent-Kung：分归约和分发两阶段，工作量 ~2N，接近串行。步数略多但更高效。

分层扫描：大数据需要多 Block 协作。三阶段：Block 内扫描 → 扫描 Section Sums → 加偏移。可递归处理超大数据。

工作效率：并行算法的总操作数不应远超串行。GPU 处理器多，工作效率比步数更重要。

核心应用：流压缩、基数排序、负载分配。前缀和把"输出位置"问题变成并行可解。

前缀和是并行计算的基础原语，理解它对于实现复杂并行算法至关重要。下一章学习合并（Merge）——另一个重要的并行原语。

参考资料：

  • Hwu, W., Kirk, D., & El Hajj, I. (2022). Programming Massively Parallel Processors: A Hands-on Approach (4th Edition). Morgan Kaufmann.
  • Blelloch, G. E. (1990). Prefix Sums and Their Applications. CMU Technical Report.
  • NVIDIA CUB Library - Scan

本文 GitHub 仓库: https://github.com/psmarter/PMPP-Learning