极客工场

Appearance

共 5,379 字预计阅读约 18 分钟

并行前缀和 scan：Hillis-Steele、Blelloch 与 warp shuffle 三种写法 ​

更新时间：2026-09-25。写并行程序时，几乎每个程序员都遇到过那个「看起来没法并行」的循环：第 i 个输出依赖前面所有输入。前缀和（prefix sum，GPU 社区习惯叫 scan）就是这个循环的学名。它是 CUDA 线程层次 之上最经典的协作型 kernel：一个 block 内几百个线程要通过 shared memory 和屏障反复交换中间值。硬件侧为什么 warp 适合干这件事，见 SIMT 与分支分歧；warp 调度器怎样把这些只有几步的 kernel 铺到 SM 上，见 GPU 架构与 warp 调度器。有意思的是，数字电路里早就有同一个思想——超前进位加法器的进位链就是一个前缀网络，见 高级加法器。

本文要回答的问题 ​

  • 前缀和的 inclusive 与 exclusive 两种定义差在哪，为什么 stream compaction 需要的是 exclusive？
  • 串行前缀和明明只有 N-1 次加法，为什么 Hillis-Steele 算法要做 N log N 次加法，步数反而是 log N？
  • Blelloch 算法的「上扫 reduce + 下扫 distribute」各自做什么，为什么总工作量能压到约 2N？
  • 一个 block 只有 shared memory 能协作，block 之间不能同步，上百万元素的全局 scan 怎么拼出来？
  • warp 内 32 个 lane 用 __shfl_up_sync / __shfl_xor_sync 做扫描，为什么只要 5 步、连 shared memory 都不用？

一、前缀和是什么：inclusive 与 exclusive ​

给定数组 x[0..N-1]，inclusive scan（包含自身的前缀和） 定义为：

bash
y[i] = x[0] + x[1] + ... + x[i]        // 即 Σ_{j=0}^{i} x[j]

exclusive scan（不包含自身） 则错开一位：

bash
z[i] = x[0] + x[1] + ... + x[i-1]      // 约定 z[0] = 0

两者只差一次搬运：y[i] = z[i+1]，或者等价地 z[i] = y[i-1]（边界补 0）。任何一种算法算出其中一个，挪动一格就得到另一个，不需要重算。

拿一个 8 元素的小数组走一遍。设 x = [3, 1, 4, 1, 5, 9, 2, 6]，inclusive scan 的结果是 [3, 4, 8, 9, 14, 23, 25, 31]，exclusive scan 则是 [0, 3, 4, 8, 9, 14, 23, 25]。对比两列能看出 exclusive 就是把 inclusive 整体右移一格、开头补 0，末尾的总和 31 不属于任何一格——它通常作为副产品被单独取走。这个例子也能直接演示 stream compaction：假设筛选条件是「大于 2」，谓词数组为 [1, 0, 1, 0, 1, 1, 0, 1]，对谓词做 exclusive scan 得到 [0, 1, 1, 2, 2, 3, 4, 4]；于是五个有效元素 3、4、5、9、6 各自的持有者，按对应位置上的扫描值 0、1、2、3、4 写进输出数组，压紧后恰好是 [3, 4, 5, 9, 6]，没有冲突、不需要原子操作。

串行写法是一个再普通不过的累加循环，N-1 次加法、O(N) 工作量，但依赖链长 N-1：第 i 次加法必须等第 i-1 次完成。在单核上这是最优的，可在几千个线程的 GPU 上，这条链意味着所有线程只能看着一个线程干活。并行前缀和要解决的，就是把这条 N-1 长的链重排成一棵深度 log N 的树，让链上的远距离依赖通过「跨步相加」在几轮内同时完成。

顺带厘清一个术语差异：CPU 并行算法文献里叫 prefix sum 或 prefix scan，NVIDIA 文档里通称 scan，而数据库和编译器文献常把它叫做 prefix aggregation，指的都是同一个东西。scan 的算子默认是加法，但形式化定义里它可以是任意结合算子（associative operator），后面第四节会用到这一点。

为什么非要跟这个循环较劲？因为它是一大批并行算法的内层构件。stream compaction（流压缩） 要把数组里满足条件的元素压紧到输出前段：先对谓词（满足为 1）做 exclusive scan，每个有效元素直接拿到自己的输出下标，再写一次即可。radix sort（基数排序） 每一位分桶时，先算出当前位为 0 的谓词，对它 scan 得到 0 组的目标位置；1 组的位置再加上 0 组总数，于是两类元素各走各的下标、一次写齐，不需要比较、也不需要锁。每一轮处理一位或几位，32 位整数排序通常十几轮内完成，每轮的核心都是 scan。求括号是否匹配则把左括号记 +1、右括号记 -1，scan 的过程中若出现负值说明右括号多余，末尾不为零说明左右不配对，一次扫描同时回答两个问题。此外多项式求值（Horner 法则的并行版）、稀疏矩阵 CSR 格式的行偏移数组构建、并行调度器给任务包分配输出区间，内核都是同一个 scan。可以说凡是「让每个元素知道前面有多少东西、从而知道自己去哪」的场合，scan 都在背后。

二、Hillis-Steele：跨步邻居相加，log N 步 ​

最直观的并行化由 Hillis 和 Steele 在 1986 年给出。第 k 步（k 从 0 开始），每个元素把自己和左边距离 2^k 的元素相加：

bash
for k = 0, 1, ..., log2(N)-1:
    并行地对所有 i：
        若 i >= 2^k，则 x[i] = x[i] + x[i - 2^k]

以 N=8 为例，传播过程如下图，没有画出箭头的格子在该步保持原值：

chip-design/intermediate/parallel-prefix-scan.md D2 示意图

图里只画了几个有代表性的格子，关键在箭头跨度的变化：步 1 每个结果汇总连续 2 个输入，步 2 汇总连续 4 个，步 3 汇总连续 8 个。右下角在步 3 就拿到了全部 8 个元素的总和——而串行链走到这里需要 7 拍。为什么不会漏？因为跨步是倍增的：步 k 结束时，位置 i 的值恰好覆盖长度 2^(k+1) 的窗口，下一步把两个这样的窗口拼起来。一共 log2 N 步，对 N=8 是 3 步，对 N=1024 是 10 步。

还可以用前面的 8 元素数组亲手验证一步：跨 1 之后数组变为 [3, 4, 5, 5, 6, 14, 11, 8]（每格是相邻两数之和），跨 2 之后再变为 [3, 4, 8, 9, 11, 19, 16, 19]，跨 4 之后得到 [3, 4, 8, 9, 14, 23, 25, 31]，与串行结果完全一致。注意第一步、第二步的中间结果看起来「不像前缀和」，比如中间出现过 5、14、11 这些数，但每个数始终是某个连续窗口的和，窗口每步扩大一倍，三步后窗口覆盖全数组。

这种「跨步倍增、每格都是一个窗口和」的连接方式，在数字电路里有个更响亮的名字——Kogge-Stone 前缀网络。用 4 条 lane、两级跨步画出它的层级连接结构：

chip-design/intermediate/parallel-prefix-scan.md D2 示意图

读这张图抓住两条主线：一是每个节点最多扇入两条边——无论 N 多大，每拍只做一次两输入加法，这就是 Kogge-Stone 网络在电路里好落地的原因；二是跨步逐层翻倍——第 1 级只连邻居（跨 1），第 2 级连距离 2 的格子，若有第 3 级就连距离 4。lane 数扩到 8、16 时图案完全平移复制，层数只按 log2 N 增长，这与 高级加法器 里超前进位链用的是同一张网。

代价在工作量上。每步所有 N 个线程都可以发起一次加法（边界附近的空转忽略），总加法数约为 N log2 N。N=1024 时约 1 万次加法，而串行只需要 1023 次——多做了整整 10 倍；N 到一百万时是 20 倍。GPU 上这些多出来的加法本身很便宜（整数加法只占发射槽，不占访存带宽），真正贵的是每步之间两次屏障的等待。用并行算法的术语说，它的步深（step complexity） 是 O(log N)，但工作量（work） 是 O(N log N)，相对串行下界 O(N) 不是 work-efficient（工作量最优）的。它的好处是简单：block 内把数组放进 shared memory，每步先读齐邻居、过一次屏障、再写回即可：

cuda
__shared__ float s[1024];
s[i] = x[i];
__syncthreads();
for (int k = 1; k < n; k <<= 1) {
    float v = s[i];
    float u = (i >= k) ? s[i - k] : 0.0f;
    __syncthreads();   // 保证所有人读完旧值
    s[i] = v + u;      // 再统一写新值
    __syncthreads();
}

注意屏障的位置：必须先让全部线程把 s[i] 和 s[i-k] 读进寄存器，再写回。漏掉中间这个屏障，写得快的线程就会覆盖邻居还没读走的旧值，结果随机出错。另一种等价写法是双缓冲：准备 s0、s1 两块 shared 数组，每轮从一块读、写另一块，轮次间交换指针。双缓冲省掉了「读后屏障」这一轮等待，代价是多占一倍 shared memory——对 scan 这种数组长度常取 1024 的 kernel，float 下就是多 4KB，多数 block 配置还承受得起；寄存器或 shared 吃紧、occupancy 已经掉到边缘时，则退回单数组版本。选哪种本质上是在 shared 容量和屏障次数之间做交换。

三、Blelloch：上扫 reduce，下扫 distribute，2N 工 ​

想把工作量从 N log N 压回 2N，需要换一棵树形算法。Blelloch scan（也叫 work-efficient scan）分两遍走同一棵二叉树，先上后下，下面以 N=8 为例。

上扫（up-sweep，reduce 阶段） 从叶子向根走：每一级把左孩子的段和加进右孩子，根节点最终得到全体总和。

chip-design/intermediate/parallel-prefix-scan.md D2 示意图

上扫的加法数是 N/2 + N/4 + ... + 1 = N-1，和一次普通树形归约完全相同。

下扫（down-sweep，distribute 阶段） 是算法的关键。先把根节点清空为 0（这一步决定了算出的是 exclusive scan），然后从根向叶子走：每一级上，父节点把自己的值传给左孩子，把「左孩子原值 + 自己的值」传给右孩子。

chip-design/intermediate/parallel-prefix-scan.md D2 示意图

下扫的分发规则乍看古怪，值得把含义点透：树上下传到某个内部节点的值，恰好是「该节点所管区间左边全部元素」的和。那么左孩子区间的元素，左边还是这些东西，直接照传；右孩子区间的元素，除了这些还要加上左孩子整段的和，也就是左孩子在上扫阶段记下的那个段和。「左传父值、右传父值加左段和」就完成了段和向两半边的交接。仍用 [3,1,4,1,5,9,2,6] 跟踪：根节点先置 0；下扫一级后，管理前四个元素的节点拿到 0、管理后四个的节点拿到前半段和 9；再下一级，第六、七两个元素所在节点拿到 9+9=18 这类偏移；叶子层最终就是前面列出的 exclusive 结果 [0,3,4,8,9,14,23,25]。

下扫同样是 N/2 + N/4 + ... + 1 = N-1 次加法（每个内部节点做一次求和）。两遍合计 2N - 2 次加法，与串行的 N-1 只差一个约 2 的常数因子，所以叫 work-efficient。步深是两遍树高，共 2 log2 N 步，比 Hillis-Steele 多一倍。这里有个重要的权衡：在 GPU 上「步」意味着一次屏障、一次全 block 的等待，Blelloch 虽然总加法少，但每步只有少数线程在干活（上扫第一步一半线程活跃，越往上越稀），负载不均；N 不大时，简单饱满的 Hillis-Steele 在实测中常常不输甚至更快。另一个细节是前缀和并不只属于加法：任何满足结合律的算子都能扫，比如最大值、异或、矩阵连乘；但算子没有交换律时（矩阵乘法），必须严格保持左右次序，跨步相加的方向不能反。

写成代码，数组长度需要是 2 的幂，不足则补零。上扫：

cuda
for (int d = 1; d < n; d <<= 1) {
    if ((i + 1) % (2 * d) == 0) s[i] += s[i - d];
    __syncthreads();
}
if (i == n - 1) s[i] = 0;      // 置零，得到 exclusive
__syncthreads();

下扫把步长反过来，每个命中节点先保下左孩子旧值再分发：

cuda
for (int d = n / 2; d >= 1; d >>= 1) {
    if ((i + 1) % (2 * d) == 0) {
        float t = s[i - d];
        s[i - d] = s[i];       // 左孩子拿父值
        s[i] += t;             // 右孩子拿父值+左孩子旧值
    }
    __syncthreads();
}

想要 inclusive 版本，要么不做根节点置零、按 inclusive 规则下扫，要么更简单——exclusive 结果每格再加回原始输入。

四、超出一个 block：中间数组与递归 scan ​

一个 block 的 shared memory 和屏障只能覆盖本 block 的线程，block 之间在单次 kernel 启动内没有屏障（CUDA 线程层次 解释了这条硬件约束）。全局 scan 因此拆成三段：

  1. 块内扫描：每个 block 对自己分到的一段（典型 1024 或 2048 个元素）用 Blelloch 做 exclusive scan，同时把本块的总和写进一个全局中间数组 totals[blockIdx.x]。
  2. 扫描块总和：对 totals 数组再做一次 scan，得到每个 block 的前缀块和——也就是该块的统一偏移量。block 数量不多时这仍是一个 scan kernel。
  3. 加偏移：再启动一个 kernel，每个 block 把自己的偏移量加到本块全部元素上。

如果数据大到中间数组本身也装不进一个 block（例如 N=1600 万、每块 1024，中间数组约 1.6 万元素），第 2 步就对 totals 递归套用同一套三段流程，形成多级；实践中两三级就收敛。用刚才的 8 元素例子类比多级流程：把它想成两个各含 4 元素的 block，第一块总和是 9、第二块总和是 22，对 [9, 22] 扫描得偏移 [0, 9]，于是第二块的局部 exclusive 结果 [0, 5, 14, 16] 整体加上 9，得到全局的 [9, 14, 23, 25]，与一次扫描的结果逐格相同。块内结果和全局结果之间只差一个统一偏移，这正是三段式正确性的核心。

CUB 库的 DeviceScan 内部就是这种「block scan + 单块扫描中间数组 + 偏移修正」的三 kernel 结构，还会用单 block 占用尽量多元素的策略减少级数。拆成三次启动要付 launch 开销（每次几微秒到十几微秒），换来的是无需全局屏障、任意规模都能跑；曾经有人尝试用全局原子标志让 block 间在一个 kernel 内自旋等待，实测延迟和可移植性都更差，主流实现都选择了多级启动。shared memory 的 bank conflict 在这套流程里影响明显：Blelloch 跨步访问容易让多个线程落到同一 bank，常用填充（padding）或调整索引规避，内存层级的代价分析见 GPU 内存层次与带宽。

五、warp shuffle：32 元素只要 5 步 ​

最小规模的 scan 发生在一个 warp（32 个 lane）内部。此时连 shared memory 都可以不用——warp 内 lane 天然锁步执行，直接用寄存器级的 shuffle 指令交换数据。Hillis-Steele 的跨步相加对应五条 __shfl_up_sync：

cuda
int v = x[lane];
v += __shfl_up_sync(0xffffffff, v, 1);   // 跨 1
v += __shfl_up_sync(0xffffffff, v, 2);   // 跨 2
v += __shfl_up_sync(0xffffffff, v, 4);   // 跨 4
v += __shfl_up_sync(0xffffffff, v, 8);   // 跨 8
v += __shfl_up_sync(0xffffffff, v, 16);  // 跨 16
// v 即 32 元素 inclusive scan；exclusive 再向上挪一位、lane0 置 0

5 条指令、5 拍，数据始终在寄存器里流动，没有 shared memory 的读写、没有 bank conflict、也不需要屏障（warp 内天然锁步，执行到此处时全部 lane 同步推进）。等价地也可以用 __shfl_xor_sync 按 1/2/4/8/16 异或掩码传播，写法对称，在归约类操作里更常见。

__shfl_xor_sync 的异或版本里，lane 与 lane lane ^ delta 交换并累加，delta 依次取 1、2、4、8、16；五次之后每格也持有全 warp 的聚合，写法在「每个人都要拿到同一个总和」的归约场景更顺手，而 scan 一般用 up 版本直接得到递增序列。要从 inclusive 变 exclusive，只需把结果再向上挪一位：e = (lane == 0) ? 0 : __shfl_up_sync(mask, v, 1)。

整个 block（比如 256 线程 = 8 个 warp）做 scan 时，常用两级拼装：每个 warp 先用 shuffle 扫出自己的 32 个结果和 warp 总和；8 个 warp 总和只有 8 个数，由 lane 0 们写进 shared memory，用一个 warp 再扫一次（甚至直接串行累加，8 次加法只有 8 拍）；过一次屏障后，各 warp 把自己的偏移加回 32 个元素。算笔账：每元素 5 次 shuffle 加法加 1 次偏移加法，warp 总和的扫描摊到 256 个元素上可忽略；全程只过一两次屏障，shared memory 只存 8 个总和，流量小到可以忽略。这就是 CUB 的 BlockScan 在小 tile 上的基本套路，也是手写 scan 最值得背下来的模板。

三种写法放在一起对比，选择依据是规模和层级：

算法步深（屏障轮数）总工作量（加法数）并行饱满度适用规模与层级
Hillis-Steele（shared 版）log2 N约 N log2 N每步近全部线程工作，饱满block 内小 N（数百到一两千），代码简单时首选
Blelloch（work-efficient）2 log2 N2N − 2上扫下扫逐步变稀，不均大 N、工作量敏感处；递归 scan 的块内内核
warp shufflelog2 32 = 5每 lane 5 次，共 160 次warp 全 lane 锁步，无 shared 流量单 warp 32 元素；block 级 scan 的第一级

以 N=1024 记几个具体数字：串行 1023 次加法、依赖链长 1023；Hillis-Steele 10 步、约 1 万次加法；Blelloch 20 步、2046 次加法；warp shuffle 对 32 个元素 5 条指令出结果。步深和工作量是两个独立坐标轴，选算法就是按数据规模在这两轴上取舍。实践中常见的混合搭配是：warp 内一律 shuffle，block 内用 Hillis-Steele 或两级 shuffle，跨 block 才请出 Blelloch 递归——让每个层级用自己最划算的那一种，而不是一个算法包打天下。

常见坑 ​

  • inclusive / exclusive 混用：stream compaction 的写入下标必须是 exclusive，用了 inclusive 会整体错位一格，还会丢掉第一个元素。
  • 原地覆盖：Hillis-Steele 每轮必须先读齐两个旧值、过屏障再写；直接 s[i] += s[i-k] 会读到被邻居提前改坏的值。
  • N 不是 2 的幂：Blelloch 的索引假设数组长度为 2 的幂，要么补零到下一个幂、要么在边界条件上单独处理，否则上扫阶段就会越界。
  • 浮点不结合：浮点加法不满足结合律，并行重排后的结果与串行累加末位可能不同；要求可复现结果时考虑定点、Kahan 补偿或接受并行口径。
  • 中间数组未归零：递归 scan 时 inactive 的 padding 元素和 totals 边界要显式置 0，否则第二级扫描会把垃圾值传播成全局偏移。
  • 块总和取错格子：exclusive scan 不会把块总和留在任何输出格里，必须在置零根节点之前单独保存它；忘了保存，第二级扫描拿到的全是 0，偏移修正形同虚设。
  • 屏障数错配：循环内的条件分支只让部分线程写 shared，但 __syncthreads() 必须所有线程都到达，不能把屏障放进 if 里，否则部分线程绕开屏障会导致整个 block 挂死。

这一篇在知识体系里的位置 ​

这一篇属于「集成电路大类 · 芯片设计」进阶层，是理解 GPU 协作型 kernel 的枢纽：它把 SIMT 与分支分歧 的锁步执行模型用成了算法构件，也为更大规模的数据并行算法铺路。线程怎样被组织成 grid 和 block，见 CUDA 线程层次；shared memory、bank conflict 所在的存储金字塔，见 GPU 内存层次与带宽；SM 怎样调度这些只有几步的 kernel、warp 切换怎样藏住延迟，见 GPU 架构与 warp 调度器。在电路侧，同一个前缀网络早已出现在 高级加法器 的超前进位逻辑里——算法与硬件在这里共享同一张树。

总结 ​

并行前缀和把 N-1 长的串行依赖链重排成 log N 深的传播树：Hillis-Steele 跨步邻居相加，log N 步、N log N 工、简单饱满，适合 block 内小规模；Blelloch 上扫归约加下扫分发，2 log N 步、约 2N 工，是大规模递归 scan 的内核；warp 内用 __shfl_up_sync 五步在寄存器里扫完 32 个元素，block 再大就靠中间数组把各块总和递归扫描、回加偏移——stream compaction、radix sort 这些「并行去依赖」的算法，内核都是它。