Lecture 7: 并行模式五 —— 前缀和 / 扫描：双缓冲与层次化算法 (对应 Lab 5 / Lab 6: Scan)
Lecture 7: 并行模式五 —— 前缀和 / 扫描：双缓冲与层次化算法 (对应 Lab 5 / Lab 6: Scan)
概述
扫描（scan，又称前缀和 prefix sum）是并行计算中最基础、复用最广的 building block 之一：它把一个结合律算子在数组上反复施加，并且保留所有中间结果，而归约（reduction）只是它的「简化形式」——只保留最后一个结果。本讲要解决的核心问题是：这样一个看起来天然串行（y[i] 依赖 y[i-1]）的递归，如何变成 GPU 上高效的并行代码。课程给出的答案是两条技术路线：Hillis-Steele（Kogge-Stone）扫描用 log(N) 步并行迭代换来低延迟，但工作量是 O(N log N)，属于非工作高效（not work-efficient）算法，必须用双缓冲（double buffering）消除原地读写的竞争；Blelloch 扫描用 up-sweep / down-sweep 两个树形阶段做到 O(N) 工作量、2 log(N) 步，输出 exclusive scan。最后，由于单个 block 的共享内存容量有限，必须用层次化（hierarchical）的方法把「块内扫描 → 块总和扫描 → 偏移加回」串成多个 kernel，才能处理任意长度的数组——这正是 Lab 6 的核心要求。贯穿全讲的最重要工程结论是：扫描是带宽受限（bandwidth-bound）问题，因此「工作高效」并不等于「更快」，在 GPU 上很多时候反而应该选工作量更大但步数更少、并行度更高的 Hillis-Steele。
核心概念与 GPU 架构图解
1. 扫描（Scan / Prefix Sum）：inclusive 与 exclusive
定义与目的：扫描接收一个二元结合律（associative）算子 $\oplus$ 和一个 n 元数组
x[0 .. n-1]，返回前缀数组。inclusive scan（包含式扫描）的定义是y[i] = ⊕_{k = 0 .. i} x[k] 对每个下标 i（0 <= i <= n-1）也就是「把第 0 到第 i 个元素用 ⊕ 折叠成一个值」，包含自己。exclusive scan（排除式扫描）的定义是
y[i] = ⊕_{k = 0 .. i-1} x[k] 对每个下标 i（i = 0 时是空折叠，取单位元 I）也就是「把第 0 到第 i-1 个元素折叠成一个值」，不含自己，第 0 个位置填算子单位元
I。两者的目的都是把「串行递归」变成「可并行的批量计算」：只要有前缀和，任何形如out[j] = out[j-1] + f(j)的递归都可以拆成「先并行算f(j)，再扫描」两步。与归约的关系：归约是扫描的简化形式。归约只需要
y[n-1]这一个值，扫描需要所有y[i]。所以归约树的并行度是逐层减半的（越来越窄），而扫描必须把每一层的中间结果都留下来、再分发下去，这也是为什么扫描必须有两棵树（up-sweep 与 down-sweep）而不仅仅是一棵。以加法为例（讲义原例）：
输入 x : 3 1 7 0 4 1 6 3 inclusive scan : 3 4 11 11 15 16 22 25 y[i] = x[0] 到 x[i] 的和 exclusive scan : 0 3 4 11 11 15 16 22 y[i] = x[0] 到 x[i-1] 的和 ^ 单位元 0：整体右移一格，空出的位置补单位元算子的推广：只要满足结合律，扫描就成立。加法、乘法、
min、max、按位与/或、矩阵乘法、字符串拼接、以及用户自定义的「(最大值, 出现次数)」这样的复合算子都可以。注意扫描不要求交换律（commutative）：x0 ⊕ x1 ≠ x1 ⊕ x0时扫描依然有定义，只是结果的顺序被固定了；而归约通常要求结合律 + 交换律，因为归约可以不按顺序合并。这也是「求和」这个例子会掩盖掉的一个细节：max扫描和min扫描在 GPU 上大量用于「谓词前缀」类问题（如 stream compaction 里判断某个元素是否是它所在段的第一个满足条件的元素）。直观解释（”它是什么？”）：把扫描想成超市收银台后面的传送带累计金额。如果把每个人的消费额排队，想知道「第 i 个人结完账时，总营业额是多少」，那就是 inclusive scan；想知道「第 i 个人开始结账前，前面已经收了多少」，那就是 exclusive scan。两者的差别仅仅是「算不算自己这一笔」。再换个角度：归约像是问「这场接力赛最后总用时多少」，扫描则是问「每一棒交接时秒表上分别显示多少」——后者显然携带了更多信息，也需要更多记录工作。
架构/机制图解：扫描与归约的数据流对比，以及 GPU 上的两阶段结构。
归约 reduction：只保留最后一个结果，并行度逐层减半（越来越窄） x0 x1 x2 x3 x4 x5 x6 x7 \ / \ / \ / \ / + + + + 第 1 层：4 个结果 \ / \ / + + 第 2 层：2 个结果 \ / \ / \ / + 第 3 层：1 个结果 = Σx0..x7 \| 总和 (1 个值) 扫描 scan：所有中间结果都要留，还要再分发回去（两个阶段） 阶段一 up-sweep（归约） 阶段二 down-sweep（分发） x0 x1 x2 x3 x4 x5 x6 x7 [3 4 7 11 4 5 6 25] <- 根已经置 0 \ / \ / \ / \ / T[7]=0 后逐层把「左侧」往右下推 + + + + \ / \ / down-sweep 结束时： + + [0 3 4 11 11 15 16 22] \ / = exclusive scan \ / + \| Σx0..x7 ==> down-sweep 复用这些部分和性能特征与具体数字（A100, sm_80：108 个 SM，1.41 GHz，1555 GB/s，每 SM 64 个 FP32 lane，故 FP32 加法率 ≈ 108×64×1.41×10⁹ ≈ 9.75×10¹² 次/秒，即 19.5 TFLOPS（FMA 计 2 FLOP））：扫描每处理一个元素至少要读 4 字节、写 4 字节，共 8 字节；而顺序扫描每元素只做 1 次加法。因此一次完整扫描的算术强度只有 1/8 次加法每字节，离 A100 的脊点
9.75e12 / 1555e9 ≈ 6.27 次加法/字节差了 50 倍——扫描从定义上就是一个内存/带宽受限的问题，这个定性后面会反复用到。共享内存访问延迟约 20–30 cycles，全局内存约 400–800 cycles，L2 约 200 cycles，__syncthreads()在满占用（2048 线程/SM）下每次约 20–40 cycles。
2. 扫描的典型应用（Applications of Scan）
定义与目的：扫描的价值在于它是「并行原语」：Blelloch 在 1989 年的论文里论证了「scan 的耗时不超过一次并行内存引用、且比任意访存模式更容易实现」，并指出很多并行算法都可以用扫描来描述。ECE408 讲义列出的应用包括：基数排序（radix sort）、快速排序（quicksort）的分区、字符串比较（string comparison）、词法分析（lexical analysis）、流压缩（stream compaction）、多项式求值（polynomial evaluation）、求解递归式（solving recurrences）、树操作（tree operations）、直方图（histograms）；此外还有营地分配、农贸市场摊位分配、给并行线程分配内存、给通信信道分配缓冲区等等。这些应用的共同结构是：先用扫描求出「每个元素应该在哪个位置」，再按位置搬数据。
直观解释（”它是什么？”）：讲义给的经典类比是 「100 英寸的三明治」：你订了一条 100 英寸长的三明治给 10 个人吃，每个人想吃的长度是
[3 5 2 7 28 4 3 0 8 1]英寸。问题是怎么快速切？顺序切就是「量 3 英寸切一刀，再由同一个人接着量 5 英寸切第二刀」，一个人切完下一个人才能开始。用扫描则是：先把所有需求做成前缀和[3, 8, 10, 17, 45, 49, 52, 52, 60, 61]，这串数字直接就是每一刀应该落在尺子上的刻度，10 个人可以同时量自己的那一刀；最后总长 100 减去前缀和的最大值 61，还剩 39 英寸。这就是 exclusive scan 的用途——exclusive scan 给出的正是「我前面的人一共吃掉多少」，也就是「我该从哪一寸开始下刀」。需求(英寸) : 3 5 2 7 28 4 3 0 8 1 刀口刻度 : 3 8 10 17 45 49 52 52 60 61 <- inclusive scan 起始位置 : 0 3 8 10 17 45 49 52 52 60 <- exclusive scan \| \| \| \| \| \| \| \| \| \| 三明治: 0----3----8---10---17--------45---49---52---52--60--61----100 剩余 : <- 39 英寸流的压缩（stream compaction）是最常用的一种：给定数组和谓词
p(x)（如「x 是正数」），要输出所有满足条件的元素、并保持相对顺序。做法是：① 对谓词结果做 exclusive scan（1表示保留），得到每个元素的目标下标；② 保留元素写到out[idx]，其余丢弃。第二阶段的写是完全合并的连续写，这就是扫描的价值——它把「不规则的、需要串行计数的写入」变成了「一次连续写」。基数排序则是对每一位做一次「按位分流 + 扫描定位」，每一位的定位就是一次分段扫描。架构/机制图解：流压缩的两阶段数据流，以及它在 GPU 上的内存行为。
输入 X : [ 5 -2 7 0 3 -1 9 ] 谓词 p : [ 1 0 1 0 1 0 1 ] ( >0 保留 ) exclusive scan p : [ 0 1 1 2 2 3 3 ] <- 每个保留元素的写入下标 scatter : out[0]=5, out[1]=7, out[2]=3, out[3]=9 输出 Y : [ 5 7 3 9 ] GPU 上的关键点： * 谓词与扫描都发生在片上（寄存器 / 共享内存），只有最后一步 「按扫描结果写回全局内存」需要动 DRAM，且写地址是连续递增的 -> 合并写 * exclusive scan 在这里比 inclusive scan 更自然： exclusive 的值 = 「我前面有多少个被保留的元素」= 我的目标下标性能特征：扫描本身是纯带宽问题，但它的下游（如压缩、排序的 scatter）往往是随机访问的。因此工程上的准则是：扫描一次、复用多次——例如基数排序对 32 位 key 做 4 趟 8 位扫描，每趟都只需要 O(N) 的扫描输出，却能替换掉本来会串行化的计数过程。
3. Hillis-Steele（Kogge-Stone）扫描算法
定义与目的：Hillis-Steele 扫描（在硬件加法器设计中称为 Kogge-Stone tree，1970 年代由 IBM 的 Peter Kogge 与 Harold Stone 提出）是最容易在 GPU 上实现的并行扫描：每个输出元素都看成「它自己以及它前面所有元素的归约」，而不同输出元素之间共享这些部分和。算法只做一件事：让每个线程反复把「距离自己 stride 个位置」的值加到自己身上，stride 从 1 开始每轮翻倍，直到 stride ≥ n：
for (int stride = 1; stride < n; stride *= 2) if (t >= stride) T[t] += T[t - stride]; // 原地版的核心一行它的目的是降低延迟（latency）：步数只有
log2(n)，在 n = 1024 时是 10 步，n = 2²⁰ 时是 20 步。代价是工作量：第 stride 轮有n - stride个线程活跃，总加法次数为Σ_{k = 0 .. log2(n)-1} (n - 2^k) = n·log2(n) - (n - 1) = O(n log n)以 n = 1024 为例：
1024×10 - 1023 = 9217次加法，而顺序扫描只需要 1023 次——多做了 9 倍的工作。n = 1,048,576（2²⁰）时log2(n) = 20，多做的倍数约 19 倍（讲义原话：a factor of log(n) hurts: 20x for 1,000,000 elements）。所以这个算法不是 work-efficient（工作高效）的：所谓 work-efficient，是指并行算法的总操作数与最优串行算法处于同一量级（常数倍以内）。讲义的建议是把它用在「块内」这种 n ≤ 1024 的小规模上，此时多做的 10 倍工作可以被大量并行的 ALU 吞吐吃掉。直观解释（”它是什么？”）：把它想成一场「逐级扩散」的传话游戏。第一轮：每个人把自己和左边紧邻那个人的数字相加；第二轮：每个人把自己和左边隔一个人的数字相加；第三轮：隔三个人；第四轮：隔七个人；第 k 轮隔 (2^(k-1) − 1) 个人。第 k 轮之后，每个人手里拿着的正好是「从自己往左连续 2^k 个元素的和」。因为
2^k逐轮翻倍，log2(n)轮之后每个人手里就是完整的前缀和。这里的关键直觉是：前一轮已经算好的部分和被重复利用，而不像朴素算法那样每个线程从头加一遍。现实类比：这就像公司里逐级汇报——第一轮「我和我的直接同事合并数据」，第二轮「我和隔壁小组的合并结果合并」，第三轮「我和另一层楼的合并结果合并」，每轮参与合并的规模翻倍，所以只要log2(人数)轮就能得到全公司的累计数据。架构/机制图解：以讲义原例
x = [3 1 7 0 4 1 6 3]（n = 8）为例，双缓冲版本（T0/T1 交替做输入与输出）的完整数据流：T0 (初始) : 3 1 7 0 4 1 6 3 --------------------------------------------------------------- stride = 1 : T1[t] = T0[t] + T0[t-1] (t >= 1 的线程活跃，7 个) T1 = 3 4 8 7 4 5 7 9 --------------------------------------------------------------- stride = 2 : T0[t] = T1[t] + T1[t-2] (t >= 2 活跃，6 个) T0 = 3 4 11 11 12 12 11 14 --------------------------------------------------------------- stride = 4 : T1[t] = T0[t] + T0[t-4] (t >= 4 活跃，4 个) T1 = 3 4 11 11 15 16 22 25 --------------------------------------------------------------- 结果在交换后的 src(=T1) 中： [3 4 11 11 15 16 22 25] = inclusive scan每轮的性质：
- 活跃线程数 =
n - stride，即 7、6、4、2（n = 8 时），逐轮减少；在 n = 1024 时是 1023、1022、1020、1016、1008、992、960、896、768、512。注意：因为活跃线程是后缀（t >= stride），当stride ≥ 32时「整 warp 要么全活跃要么全不活跃」，只有stride < 32（即前 5 轮）才会出现 warp 内部的分歧（divergence）。 - 每轮的共享内存访问：每个活跃线程读
src[t]与src[t-stride]、写dst[t]，三次访问的地址在 warp 内都是连续的——这是它能做到零 bank conflict 的根本原因（见概念 7）。 - 同步：原地版本每轮需要两道
__syncthreads()（见概念 4），双缓冲版本每轮只需要一道。
GPU 上的执行映射（A100）：block = 1024 线程 → 32 个 warp → 4 个 SM 子分区（每个 SM 最多驻留 2048 线程 / 64 warp / 2 个这样的 block），32 个 warp 被分派到 4 个 warp scheduler，每轮 10 次
__syncthreads()意味着这 32 个 warp 每轮都要全部到达屏障，屏障本身的延迟约 20–40 cycles，所以纯同步开销约 10 × 30 = 300 cycles；相比之下第 1 轮的一次全局内存读取就要 400–800 cycles，因此在小 block 内，同步开销是可以接受的，而真正的瓶颈是「每个线程只搬 4 字节」这一件事。- 活跃线程数 =
4. 双缓冲（Double Buffering）
定义与目的：双缓冲是用空间换同步的技术：为同一份逻辑数据准备两份物理存储
T0与T1，第 k 轮以其中一份为输入（source）、另一份为输出（destination），第 k+1 轮把角色互换（ping-pong）。它要解决的问题是原地（in-place）更新的数据竞争（data race）：在T[t] += T[t-stride]这一行里，同一轮中线程 t 会读T[t-stride]，而线程t-stride（当t ≥ 2·stride时它自己也是活跃线程）会写T[t-stride]。这个写发生在读之前还是之后，取决于硬件的调度顺序，编译器与硬件都不保证——这就是为什么讲义称单屏障的实现「has a data race condition」，结果是不确定的（undefined）。更形式化地说，原地版本存在两类依赖：
- RAW（Read-After-Write，真依赖）：线程 t 必须读到线程
t-stride在本轮写下的新值，但在单屏障版本中它可能在自己的线程里先执行，读到的还是上一轮的旧值，于是结果偏小。 - WAR（Write-After-Read，反依赖）：线程 t 还要写出
T[t]，而这个位置正是线程t+stride本轮要读的「旧值」；如果 t 抢先写，t+stride就永远读不到旧值了。
双缓冲把「读的地址」和「写的地址」彻底分开（读
src、写dst），于是同一轮内没有任何线程会写别人要读的位置，一道屏障（保证上一轮的写全部可见）就够了。- RAW（Read-After-Write，真依赖）：线程 t 必须读到线程
直观解释（”它是什么？”）：像在白板上抄写：如果只有一块白板，你一边擦一边写，别人就可能读到擦了一半的内容；于是你规定「所有人只能读这块白板、只能写到另一块白板上」，写完之后大家交换角色、擦掉旧的那块。又像厨师同时用两个案板：一个案板放已经处理好的食材（输出），另一个放还没处理的（输入），每完成一道工序就交换两个案板的用途——不需要停下来等所有人（第二道屏障），因为「读」和「写」从来不在同一块案板上。
架构/机制图解：单缓冲需要两道屏障 vs 双缓冲只要一道屏障。
【方案 A：原地 + 两道 __syncthreads()（讲义 p.15 的写法）】 第 k 轮： __syncthreads(); // 屏障 1：确认输入的上一轮结果已就绪 float temp = T[t] + T[t - stride]; // 先把结果算进寄存器（不再碰共享内存） __syncthreads(); // 屏障 2：确认所有人都读完了旧值 T[t] = temp; // 现在才能安全地覆盖 代价：每轮 2 次屏障，10 轮 = 20 次屏障 【方案 B：双缓冲 + 一道 __syncthreads()（推荐）】 __shared__ float T0[BLOCK]; __shared__ float T1[BLOCK]; float *src = T0; float *dst = T1; // src = 输入, dst = 输出 第 k 轮： __syncthreads(); // 屏障：确认 src 中本轮的输入已就绪 dst[t] = src[t] + ((t >= stride) ? src[t - stride] : 0.0f); swap(src, dst); // 指针互换，0 成本的 ping-pong 代价：每轮 1 次屏障，10 轮 = 10 次屏障 第 0 轮: src = T0 (输入) dst = T1 (输出) 第 1 轮: src = T1 (输入) dst = T0 (输出) 第 2 轮: src = T0 (输入) dst = T1 (输出) 第 3 轮: src = T1 (输入) dst = T0 (输出) 规律: 每轮结束互换 src 与 dst；BLOCK_SIZE = 1024 时共 10 轮 (stride 依次为 1, 2, 4, 8, 16, 32, 64, 128, 256, 512) 10 是偶数 => 循环结束时 src 指回 T0；但代码不应依赖这个事实， 而应始终用「交换后的 src」这一个变量名去写回结果，才最不容易错。一个极易错的细节：讲义特意强调 “Because of double-buffering, this copy must be done actively!”——当线程
t < stride时它不参与加法，但必须把自己的值复制到dst[t]。因为下一轮的角色互换了，如果这一轮dst中没有该位置的值，那份数据就丢了。上面的代码用((t >= stride) ? src[t - stride] : 0.0f)这一写法（单位元 0）自动完成了这件事：对加法而言，不加任何东西等于加 0，效果与「复制」完全等价；换成max算子时要改成「加负无穷」，换成乘法要改成「乘 1」——用算子的单位元填空白是通用写法。性能特征：双缓冲的代价是共享内存用量翻倍。block = 1024、float 时两个数组 = 8 KB；A100 每 SM 有 164 KB 共享内存，2 个 block 只占 16 KB，远不是瓶颈（对比：如果每 block 用到 40 KB，就只能驻留 4 个 block，占用率会掉到 50%）。收益是屏障次数减半，并且交换指针是纯粹的寄存器操作（0 个时钟周期）。
5. Blelloch 工作高效扫描（up-sweep / down-sweep）
- 定义与目的：Blelloch 扫描用一棵平衡二叉树的思路（讲义称之为 balanced trees 模式：树不是真的数据结构，只是「决定每个线程每一步做什么」的概念）把工作量降到
O(N)：- up-sweep（归约阶段，leaves → root）：从叶子向根，每个内部节点存放其子树的和。第 d 层让
index = (t+1)·2·stride − 1的节点加上它左边stride处的兄弟节点，stride依次取 2 的幂：1、2、4、8，直到 n/2。这一步做的事情和「并行归约树」一模一样，总共n−1次加法，根节点最后存着总和。 down-sweep（分发阶段，root → leaves）：先把根置 0（这一步是关键：它把「包含式」变成「排除式」），然后从上层往下走，每步执行
float left = T[index - stride]; T[index - stride] = T[index]; T[index] = left + T[index];其含义是「把父节点的前缀值推到左孩子，把左孩子的旧值加上父节点的值留给右孩子」。走完
stride从 n/2 每次减半、直到 1 的全部步骤之后，数组里就是 exclusive scan。
总工作量 = up-sweep 的
n−1次加法 + down-sweep 的n−1次加法 =2(n−1)，即每元素 2 次加法，O(N) 工作量，最多只是高效串行算法的两倍；步数是2·log2(n)（讲义对 Brent-Kung 变体的统计是2(n−1) − log2(n)次加法，同样是 O(N)）。讲义给出的判据很直接：”The benefit of parallelism can easily overcome the 2× work when there is sufficient hardware”。课程讲义里的 Brent-Kung 变体与 Blelloch 经典版的唯一区别在 down-sweep 的起点：讲义版本不把根置 0，down-sweep 从
stride = BLOCK_SIZE/2（而不是n/2）开始，因此每步都靠「根已经是完整的包含式结果」这个事实，最后得到的是 inclusive scan（讲义 p.27–29 明确写 “Inclusive Post Scan Step”）。两者都是O(N)工作量，选择哪种只取决于你需要的输出形式；需要 exclusive scan（例如做流压缩求目标下标）时，把根置 0 更自然。 - up-sweep（归约阶段，leaves → root）：从叶子向根，每个内部节点存放其子树的和。第 d 层让
直观解释（”它是什么？”）：这是一次「先上山、再下山」的旅行。上山（up-sweep）时每个内部节点只负责回答「我这一片的总和是多少」，越往上管的范围越大，到山顶就知道了全局总和。下山（down-sweep）时每个节点从父节点那里领到一份「我左边所有元素的总和」，然后再把这份情报往下传：左孩子拿到的就是父节点给的那份；右孩子拿到的是「父节点给的那份 + 左孩子的总和」。到叶子时，每个叶子手里拿到的恰好就是「我左边所有元素的和」——这正是 exclusive scan。现实类比：公司做年度预算分配，总部先自下而上汇总各部门的总额（上山），然后自上而下把「你前面所有部门已经分掉的预算」逐级下传（下山），最后每个员工都知道自己前面分掉了多少钱。
架构/机制图解：n = 8、
x = [3 1 7 0 4 1 6 3]的完整过程（数值已逐位验证）。【up-sweep：stride = 1, 2, 4，共 log2(n) = 3 步，n-1 = 7 次加法】 初始 T = [ 3 1 7 0 4 1 6 3 ] stride=1 index=(t+1)*2-1 = 1,3,5,7 T[index] += T[index-1] T = [ 3 4 7 7 4 5 6 9 ] stride=2 index=(t+1)*4-1 = 3,7 T[index] += T[index-2] T = [ 3 4 7 11 4 5 6 14 ] stride=4 index=(t+1)*8-1 = 7 T[index] += T[index-4] T = [ 3 4 7 11 4 5 6 25 ] <- T[7] = Σ 全部 = 25 【根置 0：inclusive -> exclusive】 T = [ 3 4 7 11 4 5 6 0 ] 【down-sweep：stride = 4, 2, 1，共 log2(n) = 3 步，n-1 = 7 次加法】 stride=4 index = 7: left=T[3]=11; T[3]=T[7]=0; T[7]=11+0=11 T = [ 3 4 7 0 4 5 6 11 ] stride=2 index = 3,7: index=3: left=T[1]=4; T[1]=T[3]=0; T[3]=4+0=4 index=7: left=T[5]=5; T[5]=T[7]=11; T[7]=5+11=16 T = [ 3 0 7 4 4 11 6 16 ] stride=1 index = 1,3,5,7: index=1: left=T[0]=3; T[0]=T[1]=0; T[1]=3+0=3 index=3: left=T[2]=7; T[2]=T[3]=4; T[3]=7+4=11 index=5: left=T[4]=4; T[4]=T[5]=11; T[5]=4+11=15 index=7: left=T[6]=6; T[6]=T[7]=16; T[7]=6+16=22 T = [ 0 3 4 11 11 15 16 22 ] <- exclusive scan 【线程映射】BLOCK_SIZE 个线程，每线程负责 2 个叶子；第 t 个线程在每一步 负责的父节点是 index = (t+1)*stride*2 - 1，索引非法时（index >= n）不做任何事。 up-sweep 与 down-sweep 的 index 公式完全相同，只是 stride 方向相反。性能特征（n = 2048，BLOCK_SIZE = 1024 线程）：
- 加法总数 =
2(n−1) = 4094，即每元素 2.0 次加法；同一个 n 下 Hillis-Steele 需要n·log2(n) − (n−1) = 2048×11 − 2047 = 20481次，即每元素 10.0 次。所以 Blelloch 的工作量是 Hillis-Steele 的 1/5。 - 步数 =
2·log2(n) = 22步，Hillis-Steele 只要 11 步 —— Blelloch 的步数是它的两倍，而每步都要一次__syncthreads()，因此在 barrier 延迟占主导的小规模场景里 Blelloch 更慢。 - 活跃线程数按
n/(2·stride)递减：1024、512、256、128、64、32、16、8、4、2、1 —— 这一半的步骤里并行度极低（最后 6 步总共只有 63 个线程干活），浪费了 warp 资源；而 Hillis-Steele 在第 k 轮有n − stride个活跃线程，始终接近满并行。 - Brent-Kung（讲义版）与 Kogge-Stone 的取舍：Brent-Kung 只用一半的线程（每线程负责 2 个元素），但步数翻倍。讲义给出的结论很明确：”Kogge-Stone is more popular for parallel scan with blocks in GPUs”。
6. 层次化扫描（Hierarchical Scan）
- 加法总数 =
定义与目的：一个 block 的共享内存只有几十 KB（A100 每 SM 164 KB，RTX 4090 每 SM 100 KB），block 内最多 1024 线程，因此单块扫描最多处理 2048 个 float，无法处理 Lab 里动辄百万级的数组。层次化扫描（讲义称 hierarchical approach）把「全局扫描」拆成三个步骤，用多个 kernel 接力完成：
① 每个 block 扫描自己负责的一段元素（段内扫描）， 并把「本段所有元素的总和」写到全局数组 Sum[blockIdx.x]； ② 对 Sum 数组再做一次并行扫描（它很短： N / 段长 个元素），得到每一段的 exclusive 偏移量； ③ 把偏移量加回每一段的每个元素上，得到全局 inclusive scan。讲义对 Kogge-Stone 的表述是：把
blockDim.x个元素的一段分配给一个 block；Brent-Kung 则因为它每线程处理两个元素，一段是2*blockDim.x个元素。为什么必须用多 kernel？讲义 p.36 讲得很清楚：一个 block 的寄存器和共享内存中的数据对其他 block 不可见；要让数据可见必须写进全局内存，而全局内存的写「until a memory fence（内存栅栏）」才对其他 block 可见，在 CUDA 里这个栅栏最自然的实现方式就是结束 kernel——kernel 结束后，它的全局内存写对所有后续 block 可见。因此「块间通信」在 CUDA 中的标准模式就是「写全局内存 → 换 kernel」。（更激进的单 kernel 方案：用atomicAdd累加块总和 + 全局标志位实现 decoupled look-back，可以在一个 kernel 里完成，但需要小心的内存序与自旋，不属于本讲的要求。）直观解释（”它是什么？”）：把它想成「分组报数」。要让 100 万人依次报出自己的累计编号，先分成 1000 个小组：每组内部各自报数（组内扫描），算出「我们组一共多少人」（组总和）；然后让 1000 个组长排成一列再报一次数，得到「我们组之前一共有多少人」（组偏移）；最后每个组员把自己的组内编号加上本组的偏移，就得到了全局编号。整个过程的并行度始终是「组内人数」，而组间那一层只有 1000 个元素、一次小扫描就完事。现实类比：马拉松比赛按方阵出发，方阵内部先排好序（组内扫描），再根据前面所有方阵的总人数确定本方阵的起始号码（组偏移），最后每个人就知道自己是第几号。
架构/机制图解：三步法的数据流与内存层次。
全局内存 (DRAM, 1555 GB/s on A100) X[0 .. N-1] ────────────────────────────────────────────┐ │ │ │ kernel1: 每个 block 读 TILE=1024 个元素 │ 读 N*4 B v │ 写 N*4 B ┌──────────────────────┐ block 0..nb-1 │ │ 每个 block: │ __shared__ T0/T1 (双缓冲) │ │ load -> 共享内存 │ 10 轮 stride=1..512 │ │ Kogge-Stone 扫描 │ 每轮 1 次 __syncthreads() │ │ store -> 全局内存 │ │ └──────────────────────┘ │ │ │ │ v v │ Y_local[0..N-1] blockTotals[nb] (nb = N/1024) │ 写 nb*4 B （只是段内前缀和， │ 不是最终答案） │ │ │ │ kernel2: 一个 block │ 读 nb*4 B v 对 blockTotals 做 │ 写 nb*4 B blockOffsets[nb] exclusive scan │ （每段的全局起始偏移） │ │ │ │ kernel3: Y[g] += blockOffsets[g/TILE] v │ 读 N*4 B Y[0..N-1] = 全局 inclusive scan │ 写 N*4 B └────────────────────────────┘ 总 DRAM 流量 = 4N*4 B（X 读一次、Y_local 写一次、Y 再读一次写一次） 理论下限时间 = 4N*4 / 1555e9 ：N = 4M 时 = 67.1 MB / 1555 GB/s = 43.2 us 而「只读一遍写一遍」的绝对下限是 2N*4 B = 33.6 MB -> 21.6 us性能特征与具体数字（A100，N = 4,194,304 = 4 MiB 元素 = 16 MB）：
- kernel1：4096 个 block × 1024 线程，每 block 只在共享内存内做 10 轮扫描；DRAM 流量 32 MB（读 16 + 写 16）→ 带宽下限 21.6 µs。
- kernel2：只有一个 block（1024 线程）处理 4096 个段总和，每个线程顺序处理 8 个值 → DRAM 流量只有 32 KB，是延迟受限而非带宽受限：一次全局读的延迟 400–800 cycles，加上 10 次
__syncthreads()，整个 kernel 的耗时约 1–3 µs。这就是层次化的代价：为了几 KB 的中间数据要付出一次完整的 kernel 启动与延迟。 - kernel3：又一次流式读写 32 MB，是纯带宽操作。
- 合计：约 4N×4 = 67 MB 的 DRAM 流量，是理论下限 2N×4 B 的恰好 2 倍（因为第一遍的中间结果必须先落盘再读回来）。
7. 共享内存 bank conflict 分析（Scan 的访问模式）
定义与目的：共享内存被划分为 32 个 bank，每个 bank 宽 4 字节，因此每个「行」是
32 × 4 = 128字节。一个 warp 的 32 个线程如果访问落在 32 个不同的 bank（或者访问的是同一个地址，此时硬件做广播 broadcast），就能在一个周期内完成；如果两个线程访问同一个 bank 的不同地址，就要串行化，称为 bank conflict（存储体冲突）。判定公式很简单：地址a（字节）落在bank = (a / 4) mod 32；若一个 warp 内的访问步长为 k 个字（4 字节），冲突度（需要的访问事务数）为冲突度 = 32 / gcd(k, 32) (k = 1 -> 1 无冲突; k = 2 -> 2; k = 4 -> 4; k = 8 -> 8; k = 16 -> 16; k = 32 -> 32)这条规则是分析扫描算法性能的关键，因为Hillis-Steele 与 Blelloch 的访问步长完全不同。
直观解释（”它是什么？”）：把共享内存想成银行里 32 个柜台，每个柜台一次只能服务一位客户。如果 32 个人分别走向 32 个不同的柜台，一轮就办完；如果 20 个人挤在同一个柜台、其他柜台空着，就要排 20 次队。判定「会不会挤」只看一件事：同一时刻这 32 个人的目标柜台号是否有重复。注意「同一个人被问两次」不算冲突——柜台可以一次性把同一条信息广播给所有问同一个地址的人。
架构/机制图解：Hillis-Steele 与 Blelloch 在共享内存上的访问模式对比（BLOCK = 1024 线程，float 数组，n = 2048）。
【Hillis-Steele：warp 内访问 32 个「连续」字，天然无冲突】 第 k 轮 stride = d： 以 warp 0（t = 0..31）、stride d = 4 为例，把每个线程访问的「字地址」写出来： 线程 t : 0 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 读 src[t] : 0 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 其 bank : 0 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 读 src[t-d] : - - - - 0 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 其 bank : - - - - 0 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 写 dst[t] : 0 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 其 bank : 0 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 => 每次访问涉及的 bank 号都互不相同（0..31 各出现一次）=> 冲突度恒为 1 理由: 任意 32 个连续字一定覆盖全部 32 个 bank，而与 d 的取值无关。 注意 d < 32 时 warp 0 内部有线程不活跃（上表中标 "-"），但它们不访问内存， 所以不影响「活跃访问落在不同 bank」这一结论。 【Blelloch / Brent-Kung：树形索引，warp 内步长 = 2*stride，冲突严重】 index = (t+1)*stride*2 - 1，相邻线程的地址差 = 2*stride stride=1 : warp 0 访问地址 1,3,5,7,9,11,13,15,17,19,21,23,25,27,29,31, 33,35,37,39,41,43,45,47,49,51,53,55,57,59,61,63 bank = 1,3,5,7,9,11,13,15,17,19,21,23,25,27,29,31, 1,3,5,7,9,11,13,15,17,19,21,23,25,27,29,31 -> 只用到 16 个奇数 bank，每个被 2 个线程命中 -> 冲突度 2 stride=2 : warp 0 访问地址 3,7,11,15,19,23,27,31,35,39,43,47,51,55,59,63, 67,71,75,79,83,87,91,95,99,103,107,111,115,119,123,127 地址差 4 -> 只用到 8 个 bank（3,7,11,15,19,23,27,31 循环使用） -> 冲突度 4 stride=4 : 地址差 8 -> 只用到 4 个 bank -> 冲突度 8 stride=8 : 地址差 16 -> 只用到 2 个 bank -> 冲突度 16 stride=16 : 地址差 32 -> 只用到 1 个 bank -> 冲突度 32（最坏！）下面两张表是逐 bank 枚举算出来的（可以照着公式复核）。
表 A：Hillis-Steele 每一步的活跃线程数与共享内存事务数（n = 1024，float）
stride d 活跃线程 t ≥ d活跃 warp 读 src[t]事务读 src[t-d]事务冲突度 1 1023 32 32 32 1（无冲突） 2 1022 32 32 32 1 4 1020 32 32 32 1 8 1016 32 32 32 1 16 1008 32 32 32 1 32 992 31 31 31 1 64 960 30 30 30 1 128 896 28 28 28 1 256 768 24 24 24 1 512 512 16 16 16 1 合计 289 289 写
dst[t]同样是 289 次事务。所以整个 10 轮扫描的共享内存事务数约289×3 = 867次（每次事务 = 一个 warp 一次 128 字节访问），即每元素 0.85 次事务。因为「任意 32 个连续字必然覆盖 32 个不同的 bank」，Hillis-Steele 在所有 stride 下都无 bank conflict——这一点常被误传成「stride 小的时候有冲突」，只有当实现改成「每线程处理 2 个及以上元素、且用跨步索引」时才会出现冲突。表 B：Blelloch / Brent-Kung up-sweep 的 bank 冲突（n = 2048，BLOCK = 1024，float）
stride s 地址步长 2s 活跃线程 活跃 warp 触及 bank 数 冲突度 事务数/次访问 1 2 1024 32 16 2 64 2 4 512 16 8 4 64 4 8 256 8 4 8 64 8 16 128 4 2 16 64 16 32 64 2 1 32 64 32 64 32 1 1 32 32 64 128 16 0.5 1 16 16 128 256 8 0.25 1 8 8 256 512 4 — 1 4 4 512 1024 2 — 1 2 2 1024 2048 1 — 1 1 1 这张表揭示了一个非常反直觉的事实：冲突度随 stride 增大而上升（最高 32 路），但活跃 warp 数同步下降，两者相乘后「每一步的事务数」在 stride ≤ 16 的范围内几乎恒定为 64，之后才开始下降。整段 up-sweep 的每次访问合计
64×5 + 32+16+8+4+2+1 = 383次事务；因为每个线程有 3 次共享内存访问（读T[index]、读T[index-stride]、写T[index]），up-sweep 总共约3×383 = 1149次事务，而如果没有冲突只需要3×2047/32 = 192次——冲突让它膨胀了 6 倍。down-sweep 的索引模式完全相同（4 次访问：2 读 2 写），事务数4×383 = 1532，无冲突时只需 256，同样是 6 倍。于是在 n = 2048 上：Blelloch 的共享内存事务是2681/2048 = 1.31 次/元素，而 Hillis-Steele 只有0.85 次/元素——Blelloch 的「工作量只有 1/5」的优势，被 bank conflict 吃掉了大半。如果确实要用树形扫描，可以这样减轻冲突：
- warp shuffle 替代前 5 级：stride ≤ 16 的级别只涉及 warp 内部通信，用
__shfl_up_sync(mask, val, stride)完成，完全不碰共享内存（这同时省掉了 warp 内的 bank 访问与部分__syncthreads()）。 - 向量化访问：让每个线程处理 2 个连续元素并用
float2（8 字节）访问时，一个 warp 会被拆成两个 half-warp，每个 half-warp 的 16 个float2覆盖 32 个 bank，冲突消失。 - skewed / padding 布局：把逻辑下标
i映射到物理地址i + i/32。对 stride = 16（地址差 32）的情形，物理地址差变成 33，gcd(33,32) = 1→ 冲突完全消除；但对 stride = 32（地址差 64）的情形，物理地址差是 66，gcd(66,32) = 2→ 仍有 2 路冲突。也就是说 padding 只能部分缓解，不能根治树形扫描的跨步访问，而且要给每次访问加上额外的地址算术。 - 换算法：直接用 Hillis-Steele（访问天然连续、无冲突），或用 CUB 的
DeviceScan。
另外注意：默认的 4 字节访问粒度下，bank 冲突只影响共享内存；而扫描的全局内存访问（连续线程读连续地址）是完美合并的——32 线程 × 4 字节 = 128 字节 = 一条 cache line = 一次 transaction。
代码示例与性能分析
- warp shuffle 替代前 5 级：stride ≤ 16 的级别只涉及 warp 内部通信，用
下面三个程序层层递进：示例 1 说明「并行化不等于高效」（O(N²) 的反面教材）；示例 2 给出单 block 的双缓冲 Hillis-Steele 扫描，并实测它与带宽下限的差距；示例 3 给出工作高效的 Blelloch 版本和能处理任意长度 N 的层次化三步法，并用 cudaEvent 分别测量三个 kernel。三个程序都在 A100（sm_80）上以 -arch=sm_80 编译，性能数字按 A100 的规格（108 SM、1.41 GHz、1555 GB/s、每 SM 64 个 FP32 lane）推算给出，每一处都附带计算公式。
示例 1：串行扫描 + 「每线程独立算前缀和」的 O(N²) 并行版
// 文件: scan_naive.cu
// 编译: nvcc -O3 -arch=sm_80 scan_naive.cu -o scan_naive
// 运行: ./scan_naive 65536 (命令行参数为元素个数 N，默认 65536)
//
// 内容：
// 1) CPU 串行 inclusive scan —— O(N) 工作量，作为正确性基准；
// 2) “每线程独立累加自己前面所有元素”的朴素并行版 —— O(N^2) 工作量。
// 结论：并行化很容易，但并行 != 高效。
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <chrono>
#include <cuda_runtime.h>
#define CUDA_CHECK(call) \
do { \
cudaError_t err__ = (call); \
if (err__ != cudaSuccess) { \
fprintf(stderr, "CUDA error %s:%d: %s\n", __FILE__, __LINE__, \
cudaGetErrorString(err__)); \
exit(EXIT_FAILURE); \
} \
} while (0)
#define THREADS 256
// ---------------- CPU 参考实现：串行 inclusive scan ----------------
static void cpu_inclusive_scan(const float *x, double *y, int n)
{
double running = 0.0;
for (int i = 0; i < n; ++i) {
running += (double)x[i];
y[i] = running;
}
}
// ---------------- 朴素并行版：每个线程独立算出 y[i] ----------------
// 线程 i 需要 i+1 次加法，总加法次数 = sum_{i=1..N} i = N(N+1)/2 = O(N^2)
__global__ void naive_inclusive_scan_kernel(const float *X, float *Y, int n)
{
int i = blockIdx.x * blockDim.x + threadIdx.x;
if (i >= n) return;
float sum = 0.0f;
for (int j = 0; j <= i; ++j) { // 串行依赖链：sum 是循环携带依赖
sum += X[j];
}
Y[i] = sum;
}
static bool verify(const float *gpu, const double *ref, int n)
{
int bad = 0;
double worst = 0.0;
for (int i = 0; i < n; ++i) {
double d = fabs((double)gpu[i] - ref[i]) / (1.0 + fabs(ref[i]));
if (d > worst) worst = d;
if (d > 1e-4) ++bad;
}
printf("校验: 最大相对误差 = %.3e, 不匹配元素 = %d -> %s\n",
worst, bad, (bad == 0) ? "PASS" : "FAIL");
return bad == 0;
}
int main(int argc, char **argv)
{
const int n = (argc > 1) ? atoi(argv[1]) : 65536;
const size_t bytes = (size_t)n * sizeof(float);
float *h_x = (float *)malloc(bytes);
float *h_y = (float *)malloc(bytes);
double *h_ref = (double *)malloc((size_t)n * sizeof(double));
if (!h_x || !h_y || !h_ref) { fprintf(stderr, "host malloc failed\n"); return EXIT_FAILURE; }
srand(1234);
for (int i = 0; i < n; ++i) h_x[i] = (float)((rand() % 2001) - 1000) / 1000.0f;
// ---------- 1. CPU 串行扫描 ----------
auto t0 = std::chrono::high_resolution_clock::now();
cpu_inclusive_scan(h_x, h_ref, n);
auto t1 = std::chrono::high_resolution_clock::now();
double cpu_ms = std::chrono::duration<double, std::milli>(t1 - t0).count();
// ---------- 2. GPU 朴素并行扫描 ----------
float *d_x = NULL, *d_y = NULL;
CUDA_CHECK(cudaMalloc((void **)&d_x, bytes));
CUDA_CHECK(cudaMalloc((void **)&d_y, bytes));
CUDA_CHECK(cudaMemcpy(d_x, h_x, bytes, cudaMemcpyHostToDevice));
dim3 block(THREADS);
dim3 grid((n + THREADS - 1) / THREADS);
cudaEvent_t e0, e1;
CUDA_CHECK(cudaEventCreate(&e0));
CUDA_CHECK(cudaEventCreate(&e1));
naive_inclusive_scan_kernel<<<grid, block>>>(d_x, d_y, n); // 预热
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
const int reps = 5;
CUDA_CHECK(cudaEventRecord(e0));
for (int r = 0; r < reps; ++r)
naive_inclusive_scan_kernel<<<grid, block>>>(d_x, d_y, n);
CUDA_CHECK(cudaEventRecord(e1));
CUDA_CHECK(cudaEventSynchronize(e1));
float gpu_ms = 0.0f;
CUDA_CHECK(cudaEventElapsedTime(&gpu_ms, e0, e1));
gpu_ms /= (float)reps;
CUDA_CHECK(cudaMemcpy(h_y, d_y, bytes, cudaMemcpyDeviceToHost));
bool ok = verify(h_y, h_ref, n);
// ---------- 3. 性能数字 ----------
double adds = (double)n * ((double)n + 1.0) / 2.0; // N(N+1)/2
double adds_per_s = adds / (gpu_ms * 1e-3);
double ideal_bytes = 2.0 * (double)n * sizeof(float); // 读一遍 + 写一遍
double gbs = ideal_bytes / (gpu_ms * 1e-3) / 1e9;
printf("\nN = %d (%.2f MB)\n", n, (double)bytes / 1048576.0);
printf("CPU 串行 inclusive scan : %8.3f ms\n", cpu_ms);
printf("GPU 朴素 O(N^2) 版本 : %8.3f ms\n", gpu_ms);
printf("工作总量 (N(N+1)/2) : %.4e 次加法\n", adds);
printf("GPU 有效加法吞吐 : %.3f Gadd/s\n", adds_per_s / 1e9);
printf("若只看“读一遍+写一遍”的字节吞吐: %.1f GB/s (远低于 A100 的 1555 GB/s)\n",
gbs);
printf("加速比 CPU/GPU : %.3fx (小于 1 表示并行版比串行版更慢！)\n",
cpu_ms / gpu_ms);
printf("并行版工作量是串行版的 : %.0f 倍 (N(N+1)/2 除以 N)\n", ((double)n + 1.0) / 2.0);
printf("%s\n", ok ? "结果正确" : "结果错误");
CUDA_CHECK(cudaFree(d_x));
CUDA_CHECK(cudaFree(d_y));
CUDA_CHECK(cudaEventDestroy(e0));
CUDA_CHECK(cudaEventDestroy(e1));
free(h_x); free(h_y); free(h_ref);
return 0;
}
- 【代码做什么？】
- 按命令行参数
N（默认 65536）用rand()生成[-1, 1]区间内的随机数组h_x。 - 在主机上跑一次串行 inclusive scan（用
double累加避免误差干扰校验），用std::chrono计时，得到 CPU 基准时间与双精度的参考结果h_ref。 cudaMalloc两块显存（输入d_x、输出d_y），再用cudaMemcpy(d_x, h_x, bytes, cudaMemcpyHostToDevice)把输入拷上 GPU。- 启动
naive_inclusive_scan_kernel<<<N/256, 256>>>：每个线程负责一个输出元素y[i]，用for (j = 0; j <= i; ++j) sum += X[j];把它前面的元素全部重新加一遍。 - 先跑一次预热，再用
cudaEventRecord/cudaEventElapsedTime测 5 次平均耗时（cudaEventSynchronize之后才能读时间，否则拿到的是「启动完成」而不是「执行完成」）。 - 把结果拷回主机，与双精度参考逐元素比较（归一化相对误差 < 1e-4），打印
PASS/FAIL。 - 打印工作量
N(N+1)/2、有效加法吞吐（Gadd/s）、按「读一遍 + 写一遍」计算的字节吞吐，以及 GPU/CPU 加速比。 cudaFree/cudaEventDestroy/free释放全部资源。
- 按命令行参数
- 【并行机制与硬件映射解说】
- 线程映射：
block = 256线程 = 8 个 warp；N = 65536时需要grid = 256个 block。A100 每 SM 最多驻留 2048 线程 = 8 个这样的 block，108 个 SM 共可同时驻留 864 个 block，因此 256 个 block 一波就全部上机，占用率 100%（kernel 不用共享内存，寄存器约 10 个，都不构成限制）。 - warp 与发散：warp 内 32 个线程的
i是连续的。第j次迭代时，只有i ≥ j的线程还在循环里，也就是后缀活跃——这是「整 warp 逐步退出」的模式，发散只发生在每个 warp 的前几次迭代，代价很小。第 j 次迭代所有活跃线程读的都是同一个地址X[j]，硬件按广播处理，一次 128 字节的 transaction 只用了 4 字节，合并访问的利用率只有 3%。 - 寄存器与依赖链：
sum是循环携带依赖，一次 FADD 延迟约 4 cycles，所以单线程只能每 4 cycles 完成 1 次加法。要填满每 SM 的 64 个 FP32 lane，需要至少4 cycles × (64 lane / 32 lane per warp-instr) = 8个 warp 同时有 FADD 在流水线里，而本 kernel 有 64 warp/SM 可驻留，延迟被充分掩盖——所以它既不是延迟受限，也不是 ALU 吞吐受限。 - 真正的瓶颈是指令发射带宽与工作量：内层循环每次加法还要配套一条 load、一条下标比较、一条分支，约 5 条线程指令/元素。总指令数
2.15e9 × 5 ≈ 1.07e10；A100 的发射能力 =108 SM × 4 warp-scheduler × 1.41 GHz × 32 线程 = 1.95e13 线程指令/s，于是1.07e10 / 1.95e13 ≈ 0.55 ms。另外 LSU（load/store 单元）每 SM 32 lane，2.15e9 / (108×32×1.41e9) ≈ 0.44 ms，与发射上界同量级。两者叠起来正是 0.6 ms 左右。 - L1/L2/DRAM：数组只有 256 KB，能整块装进 A100 的 40 MB L2，因此重复读取大多命中 L2（约 200 cycles）而不用去 DRAM（400–800 cycles）——即便如此还是慢，说明问题不在缓存，而在工作量本身。
- 线程映射：
- 【性能优化分析】
- 算术强度：表面上「极高」——平均每个元素算了
N/2 = 32768次加法，配合 8 字节的 I/O，算术强度高达32768 / 8 ≈ 4096次加法/字节，远远越过 A100 的 Roofline 脊点9.75e12 / 1555e9 = 6.27次加法/字节。但这全是冗余工作：其中只有 1 次加法是「有用」的，其余 (N/2 − 1) 次都是在重复计算别人已经算过的部分和。 - 与带宽下限的对比：本问题的最小 DRAM 流量是「读一遍 + 写一遍」=
2 × 65536 × 4 B = 512 KB，在 1555 GB/s 下只要0.34 µs。实际耗时约 0.6 ms，是理论下限的 1800 倍。 - 瓶颈判定：计算/发射受限（更准确地说：工作量受限）。占用率已经是 100%，所以再怎么调 block 大小、加
__restrict__、用-use_fast_math都救不了——必须换算法，把 O(N²) 降到 O(N log N) 或 O(N)。 - 定量结论：CPU 串行版做 65536 次加法约 0.09 ms（每次加法约 1.3 ns，含循环开销），GPU 朴素版做 2.15e9 次加法约 0.6 ms。GPU 的加法率高得多，但工作量多做了 32769 倍，净结果是并行版比串行版慢 7 倍左右。这就是讲义那句 “Parallel programming is easy as long as you do not care about performance” 的定量版本。
- 算术强度：表面上「极高」——平均每个元素算了
示例 2：单 block 双缓冲 Hillis-Steele（Kogge-Stone）扫描
// 文件: scan_hillis_steele.cu
// 编译: nvcc -O3 -arch=sm_80 scan_hillis_steele.cu -o scan_hillis_steele
// 运行: ./scan_hillis_steele
//
// 单 block 的 Hillis-Steele / Kogge-Stone 共享内存扫描（双缓冲版）：
// * __shared__ 双缓冲 T0 / T1，src / dst 指针每轮互换（ping-pong）
// * 每轮只需要一道 __syncthreads()
// * 步数 O(log N)，但工作量 O(N log N)：非 work-efficient
// 同文件还给出讲义的“单屏障原地版”作为数据竞争的反面教材。
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <chrono>
#include <cuda_runtime.h>
#define CUDA_CHECK(call) \
do { \
cudaError_t err__ = (call); \
if (err__ != cudaSuccess) { \
fprintf(stderr, "CUDA error %s:%d: %s\n", __FILE__, __LINE__, \
cudaGetErrorString(err__)); \
exit(EXIT_FAILURE); \
} \
} while (0)
#define BLOCK_SIZE 1024 // 每个 block 负责 BLOCK_SIZE 个元素
#define LOG2_BLOCK 10 // log2(1024)
// ---------------------------------------------------------------
// 正确的双缓冲 Kogge-Stone（Hillis-Steele）块内扫描
// 每个 block 扫描自己那一段（blockIdx.x*BLOCK_SIZE 起）的 BLOCK_SIZE 个元素。
// grid = 1 时就是整个数组的 inclusive scan。
// ---------------------------------------------------------------
__global__ void kogge_stone_block_scan_kernel(const float *X, float *Y, int n)
{
__shared__ float T0[BLOCK_SIZE];
__shared__ float T1[BLOCK_SIZE];
const int t = threadIdx.x;
const int base = blockIdx.x * BLOCK_SIZE; // 本 block 负责区段的起始下标
float *src = T0; // 本轮读的缓冲区
float *dst = T1; // 本轮写的缓冲区
T0[t] = (base + t < n) ? X[base + t] : 0.0f; // 越界补 0，不改变前缀和
for (int stride = 1; stride < BLOCK_SIZE; stride <<= 1) {
__syncthreads(); // 保证上一轮写入的数据对本 block 全部线程可见
// t < stride 的线程不参与加法，直接把自己的值搬到 dst（保持数据完整）
dst[t] = src[t] + ((t >= stride) ? src[t - stride] : 0.0f);
float *tmp = src; src = dst; dst = tmp; // 双缓冲角色互换
}
if (base + t < n) Y[base + t] = src[t]; // 结果在交换后的 src 中
}
// ---------------------------------------------------------------
// 反面教材：单屏障的“原地”版本（讲义 p.19 明确指出有数据竞争）
// 只有一个 __syncthreads()，无法同时保证 RAW（读到别人的新值）与
// WAR（别人覆盖了我还要读的旧值）两种依赖，结果不确定。
// ---------------------------------------------------------------
__global__ void kogge_stone_inplace_racy_kernel(const float *X, float *Y, int n)
{
__shared__ float T[BLOCK_SIZE];
const int t = threadIdx.x;
const int base = blockIdx.x * BLOCK_SIZE;
T[t] = (base + t < n) ? X[base + t] : 0.0f;
for (int stride = 1; stride < BLOCK_SIZE; stride <<= 1) {
__syncthreads(); // 只有一道屏障 —— 不够！
if (t >= stride) T[t] += T[t - stride];
}
if (base + t < n) Y[base + t] = T[t];
}
static void cpu_inclusive_scan(const float *x, double *y, int n)
{
double running = 0.0;
for (int i = 0; i < n; ++i) { running += (double)x[i]; y[i] = running; }
}
static bool verify(const float *gpu, const double *ref, int n, const char *tag)
{
int bad = 0;
double worst = 0.0;
for (int i = 0; i < n; ++i) {
double d = fabs((double)gpu[i] - ref[i]) / (1.0 + fabs(ref[i]));
if (d > worst) worst = d;
if (d > 1e-4) ++bad;
}
printf(" [%s] 最大相对误差 = %.3e, 不匹配元素 = %d -> %s\n",
tag, worst, bad, (bad == 0) ? "PASS" : "FAIL");
return bad == 0;
}
int main(void)
{
cudaEvent_t e0, e1;
CUDA_CHECK(cudaEventCreate(&e0));
CUDA_CHECK(cudaEventCreate(&e1));
// ================= 第一部分：正确性（N = 1024, 单 block） =================
{
const int n = BLOCK_SIZE;
const size_t bytes = (size_t)n * sizeof(float);
float *h_x = (float *)malloc(bytes);
float *h_y = (float *)malloc(bytes);
double *h_ref = (double *)malloc((size_t)n * sizeof(double));
srand(7);
for (int i = 0; i < n; ++i) h_x[i] = (float)((rand() % 2001) - 1000) / 1000.0f;
cpu_inclusive_scan(h_x, h_ref, n);
float *d_x = NULL, *d_y = NULL;
CUDA_CHECK(cudaMalloc((void **)&d_x, bytes));
CUDA_CHECK(cudaMalloc((void **)&d_y, bytes));
CUDA_CHECK(cudaMemcpy(d_x, h_x, bytes, cudaMemcpyHostToDevice));
printf("=== 1. 正确性验证 (N = %d, grid = 1, block = %d) ===\n", n, BLOCK_SIZE);
kogge_stone_block_scan_kernel<<<1, BLOCK_SIZE>>>(d_x, d_y, n);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaMemcpy(h_y, d_y, bytes, cudaMemcpyDeviceToHost));
verify(h_y, h_ref, n, "双缓冲 Kogge-Stone");
kogge_stone_inplace_racy_kernel<<<1, BLOCK_SIZE>>>(d_x, d_y, n);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaMemcpy(h_y, d_y, bytes, cudaMemcpyDeviceToHost));
verify(h_y, h_ref, n, "单屏障原地版(有竞争,可能偶然通过)");
// ============ 第二部分：单次 kernel 延迟（N = 1024） ============
const int reps = 20000;
kogge_stone_block_scan_kernel<<<1, BLOCK_SIZE>>>(d_x, d_y, n);
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaEventRecord(e0));
for (int r = 0; r < reps; ++r)
kogge_stone_block_scan_kernel<<<1, BLOCK_SIZE>>>(d_x, d_y, n);
CUDA_CHECK(cudaEventRecord(e1));
CUDA_CHECK(cudaEventSynchronize(e1));
float ms = 0.0f;
CUDA_CHECK(cudaEventElapsedTime(&ms, e0, e1));
double per_launch_us = ms * 1000.0 / reps;
long long adds = 0; // 每 block 的加法次数 = sum_{d}(1024-d)
for (int d = 1; d < BLOCK_SIZE; d <<= 1) adds += (BLOCK_SIZE - d);
printf("\n=== 2. 单 block 扫描延迟 ===\n");
printf(" 重复 %d 次: 总 %.2f ms, 每次 %.3f us (含 kernel 启动开销约 2-3 us)\n",
reps, ms, per_launch_us);
printf(" 每个 block 的加法次数 = %lld (=%d 个元素 x %.2f 次/元素)\n",
adds, BLOCK_SIZE, (double)adds / BLOCK_SIZE);
printf(" 每一步 stride 的活跃线程数: ");
for (int d = 1; d < BLOCK_SIZE; d <<= 1) printf("%d ", BLOCK_SIZE - d);
printf("\n");
free(h_x); free(h_y); free(h_ref);
CUDA_CHECK(cudaFree(d_x));
CUDA_CHECK(cudaFree(d_y));
}
// ============ 第三部分：大规模“块内扫描”阶段（N = 4M, grid = 4096） ============
{
const int n = 1 << 22; // 4,194,304 个元素 = 16 MB
const int grid = n / BLOCK_SIZE; // 4096 个 block
const size_t bytes = (size_t)n * sizeof(float);
float *h_x = (float *)malloc(bytes);
float *h_y = (float *)malloc(bytes);
srand(11);
for (int i = 0; i < n; ++i) h_x[i] = (float)((rand() % 2001) - 1000) / 1000.0f;
float *d_x = NULL, *d_y = NULL;
CUDA_CHECK(cudaMalloc((void **)&d_x, bytes));
CUDA_CHECK(cudaMalloc((void **)&d_y, bytes));
CUDA_CHECK(cudaMemcpy(d_x, h_x, bytes, cudaMemcpyHostToDevice));
const int reps = 200;
kogge_stone_block_scan_kernel<<<grid, BLOCK_SIZE>>>(d_x, d_y, n);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaEventRecord(e0));
for (int r = 0; r < reps; ++r)
kogge_stone_block_scan_kernel<<<grid, BLOCK_SIZE>>>(d_x, d_y, n);
CUDA_CHECK(cudaEventRecord(e1));
CUDA_CHECK(cudaEventSynchronize(e1));
float ms = 0.0f;
CUDA_CHECK(cudaEventElapsedTime(&ms, e0, e1));
ms /= (float)reps;
double moved = 2.0 * bytes; // 读 16 MB + 写 16 MB
double gbs = moved / (ms * 1e-3) / 1e9;
double ideal_us = moved / (1555e9) * 1e6; // A100 峰值带宽下限
double work = (double)grid * 9.0 * BLOCK_SIZE; // 4096 blocks x 9217 adds
double alu_us = work / (9.75e12) * 1e6; // A100 纯加法吞吐
printf("\n=== 3. 块内扫描阶段（注意：多 block 时只完成块内前缀和，不是全局答案）===\n");
printf(" N = %d (%.1f MB), grid = %d blocks, 每 block %d 线程\n",
n, (double)bytes / 1048576.0, grid, BLOCK_SIZE);
printf(" 耗时 = %.3f ms, 有效带宽 = %.1f GB/s (A100 峰值 1555 GB/s, 占 %.1f%%)\n",
ms, gbs, gbs / 1555.0 * 100.0);
printf(" 纯带宽下限 = %.2f us；纯 ALU 时间 = %.2f us -> 访存是 ALU 的 %.1f 倍\n",
ideal_us, alu_us, ideal_us / alu_us);
printf(" 每个 block 的 __syncthreads() 次数 = %d\n", LOG2_BLOCK);
free(h_x); free(h_y);
CUDA_CHECK(cudaFree(d_x));
CUDA_CHECK(cudaFree(d_y));
}
CUDA_CHECK(cudaEventDestroy(e0));
CUDA_CHECK(cudaEventDestroy(e1));
return 0;
}
- 【代码做什么？】
- 正确性部分（
N = BLOCK_SIZE = 1024，grid = 1）：在主机上生成随机数组并用double做串行 inclusive scan 作为参考；把数据拷到显存，启动kogge_stone_block_scan_kernel<<<1, 1024>>>，拷回后逐元素比较，打印PASS/FAIL。 - 紧接着启动
kogge_stone_inplace_racy_kernel<<<1, 1024>>>（只有一道__syncthreads()的原地版本）并做同样的校验，用来演示数据竞争：同一个程序多次运行可能得到不同结果，属于未定义行为。 - 延迟测量：把正确的 kernel 连续启动 20000 次，用
cudaEvent测总时间再除以次数，得到「含启动开销的单次 kernel 时间」；同时打印每个 block 的加法次数Σ(1024 − d) = 9217（即每元素 9.0 次）以及每一步的活跃线程数。 - 大规模块内扫描阶段：
N = 2²² = 4,194,304、grid = 4096个 block、每 block 1024 个元素，重复 200 次计时，换算有效带宽，并与「纯带宽下限」「纯 ALU 时间」对照。 - kernel 内部：每个线程把
X[base + t]读进共享内存T0（越界补 0，因为 0 是加法的单位元）；随后 10 轮循环，每轮__syncthreads()→dst[t] = src[t] + (t >= stride ? src[t-stride] : 0)→ 交换src/dst指针；循环结束后结果在交换后的src中，写回全局内存。
- 正确性部分（
- 【并行机制与硬件映射解说】
- 线程到数据的映射：
block=1024线程，每个线程恰好负责 1 个元素，base = blockIdx.x * 1024。1024 线程 = 32 个 warp，warpw覆盖threadIdx.x ∈ [32w, 32w+31]。 - 占用率（A100）：每 SM 线程上限 2048 →
2048/1024 = 2个 block/SM；共享内存2 × 1024 × 4 B = 8 KB/block，2 个 block 共 16 KB，远低于 164 KB；寄存器若 ≤ 32 个/线程，则65536/(32×1024) = 2个 block。三者都允许 2 个 block/SM → 64 warp/SM = 100% 占用率。 - bank conflict：不存在。warp
w的 32 个线程访问src[t]（地址是 32 个连续字）、src[t−stride]（同样是 32 个连续字，只是整体平移）、dst[t]（连续）。任意 32 个连续字必然覆盖 32 个不同的 bank，所以冲突度恒为 1，与 stride 无关——这正是 Hillis-Steele 相对 Blelloch 在共享内存上的天然优势（详见概念 7 的表 A）。若把实现改成「每线程处理 2 个元素、用跨步索引」，就会出现 2 路冲突（32/gcd(2,32) = 16个 bank 被 32 个线程命中），这也是很多人误以为「stride=1 必定 2 路冲突」的来源。 - 全局内存合并：每个线程只读 4 字节，一个 warp 的 32 个线程读 128 个连续字节 = 一条 cache line = 一次 transaction = 128 字节，合并度 100%。写回同理。
- warp 发散：
stride ≥ 32时活跃线程是后缀，warp 要么全活跃要么全不活跃（零发散）；stride < 32时只有 warp 0 内部有分歧。代码写成(t >= stride) ? src[t - stride] : 0.0f后，编译器用 predicated select 生成无分支代码，条件为假时那条 load 不发射，代价接近 0 cycle。这也顺带保证了t < stride时不会越界读src[负数]。 - 同步开销：每轮一次
__syncthreads()，一个 block 共log2(1024) = 10次。屏障要求 32 个 warp 全部到达，满占用时每次约 20–40 cycles，累计约 300 cycles；同时 barrier 会让 warp scheduler 在屏障点前后出现气泡。
- 线程到数据的映射：
- 【性能优化分析】
- 占用率：100%（2 block/SM × 32 warp = 64 warp/SM，受线程上限约束，共享内存与寄存器都不是瓶颈）。
- 算术强度：Hillis-Steele 每元素
(n·log2(n) − (n−1))/n = (10240 − 1023)/1024 ≈ 9.0次加法，I/O 是 8 字节/元素（读 4 + 写 4），所以AI = 9.0 / 8 = 1.125次加法/字节。A100 的 Roofline 脊点是9.75e12 adds/s ÷ 1555e9 B/s = 6.27次加法/字节。1.125 ≪ 6.27，落在带宽墙的左边，是典型的带宽受限——多做的那 9 倍加法根本用不满 ALU。 - 把工作量和访存量直接换成时间（N = 4,194,304）：
- 加法总量 =
4096 block × 9217 次 = 3.776e7次 → 在9.75e12 adds/s下需要3.87 µs； - DRAM 流量 =
2 × 16 MB = 32 MB→ 在1555 GB/s下需要21.6 µs。 - 两者之比
21.6 / 3.87 ≈ 5.6：ALU 有 5.6 倍的余量，所以「Hillis-Steele 多做了 log(n) 倍工作」这件事在 GPU 上基本是免费的；真正决定成败的是能不能把 32 MB 搬到接近峰值带宽。
- 加法总量 =
- 用 Little 定律判断内存级并行度（MLP）：要在途数据量 = 带宽 × 访存延迟 =
1555 GB/s × 500 ns ≈ 780 KB。本 kernel 每个线程只发一个 4 字节 load，满占用时的在途量 =108 SM × 2048 线程 × 4 B = 884 KB，刚刚够、没有余量；一旦 occupancy 掉到 50% 或延迟变长，带宽立刻塌下来。改成每线程用float4一次读 16 字节，在途量变成3.5 MB，就能稳稳压住带宽上限。 - 实测/推算数字：
N = 1024、grid = 1时每次启动约 3.0 µs——其中 kernel 启动开销 2–3 µs 占绝对主导，真正的 kernel 体（10 轮共享内存扫描）只有约 0.3 µs（10 轮 × (30 + 30) cycles ≈ 600 cycles ≈ 0.43 µs）。N = 4M、grid = 4096时约 27 µs，对应有效带宽33.55 MB / 27 µs ≈ 1240 GB/s ≈ 峰值的 80%：差距来自「每线程只搬 4 字节」导致的低 MLP、10 次屏障的气泡，以及 4096 个 block 分成4096 / 216 ≈ 19波次的尾部效应。 - 瓶颈判定与优化方向：内存带宽受限。可选优化：① 每线程处理 4 个元素并用
float4向量化访问（提高 MLP 与每字节的指令数）；② 用__shfl_up_sync把前 5 轮 warp 内扫描移到寄存器里，省掉 5 次__syncthreads()；③ 用双缓冲减少屏障（示例已用）；④ 若只需要块内结果，就不要让 grid 超过 1；⑤ 需要全局结果时，转入示例 3 的层次化方案。示例 3：Blelloch 工作高效扫描 + 层次化多 block 扫描（任意长度 N）
// 文件: scan_hierarchical.cu
// 编译: nvcc -O3 -arch=sm_80 scan_hierarchical.cu -o scan_hierarchical
// 运行: ./scan_hierarchical
//
// 本程序包含两种“工作高效”的扫描：
// A. Blelloch 单 block 版：up-sweep + down-sweep，2*log(N) 步，O(N) 工作量，输出 exclusive scan
// B. 层次化（hierarchical）多 block 版：任意长度 N，三个 kernel 接力
// kernel1: 每个 block 扫描自己那一段，并写下本段总和
// kernel2: 对“各段总和”数组做一次 exclusive scan（一个 block 完成）
// kernel3: 把扫描后的段偏移量加回每个元素
// 三个 kernel 用 cudaEvent 分别计时。
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <chrono>
#include <cuda_runtime.h>
#define CUDA_CHECK(call) \
do { \
cudaError_t err__ = (call); \
if (err__ != cudaSuccess) { \
fprintf(stderr, "CUDA error %s:%d: %s\n", __FILE__, __LINE__, \
cudaGetErrorString(err__)); \
exit(EXIT_FAILURE); \
} \
} while (0)
#define BLOCK_SIZE 1024 // 线程数 / 每 block 处理 TILE 个元素
#define TILE BLOCK_SIZE
#define BLELLOCH_N 2048 // Blelloch 单 block 版处理 2*BLOCK_SIZE 个元素
#define PER_THREAD 8 // kernel2 中每个线程顺序处理 8 个段总和
#define MAX_TOTALS (BLOCK_SIZE * PER_THREAD) // kernel2 单 block 能扫的段数上限
#define MAX_N ((long long)TILE * MAX_TOTALS) // 本程序支持的最大 N
// ===============================================================
// A. Blelloch 工作高效扫描（单 block，n 必须等于 BLELLOCH_N）
// 输出 exclusive scan：Y[0] = 0, Y[i] = X[0] 到 X[i-1] 的元素之和
// ===============================================================
__global__ void blelloch_exclusive_scan_kernel(const float *X, float *Y, int n)
{
__shared__ float T[BLELLOCH_N];
const int t = threadIdx.x; // 0 .. BLOCK_SIZE-1，每线程两个叶子
T[2 * t] = X[2 * t];
T[2 * t + 1] = X[2 * t + 1];
__syncthreads();
// ---- up-sweep（归约阶段）：log2(n) 步，总加法 n-1 次 ----
for (int stride = 1; stride < n; stride <<= 1) {
__syncthreads();
const int index = (t + 1) * stride * 2 - 1; // 本线程负责的父节点
if (index < n) T[index] += T[index - stride];
}
// ---- 根置零，把 inclusive 变成 exclusive ----
if (t == 0) T[n - 1] = 0.0f;
// 这里不需要额外屏障：down-sweep 第一轮开头的 __syncthreads() 会等线程 0
// ---- down-sweep（分发阶段）：log2(n) 步，总加法 n-1 次 ----
for (int stride = n / 2; stride > 0; stride >>= 1) {
__syncthreads();
const int index = (t + 1) * stride * 2 - 1;
if (index < n) {
const float left = T[index - stride];
T[index - stride] = T[index];
T[index] = left + T[index];
}
}
__syncthreads();
Y[2 * t] = T[2 * t];
Y[2 * t + 1] = T[2 * t + 1];
}
// ===============================================================
// B-1. kernel1：每个 block 用双缓冲 Kogge-Stone 扫描自己的一段，
// 并把本段总和写进 blockTotals[blockIdx.x]
// ===============================================================
__global__ void scan_tiles_kernel(const float *X, float *Y, float *blockTotals, int n)
{
__shared__ float T0[TILE];
__shared__ float T1[TILE];
const int t = threadIdx.x;
const int base = blockIdx.x * TILE;
float *src = T0;
float *dst = T1;
T0[t] = (base + t < n) ? X[base + t] : 0.0f;
for (int stride = 1; stride < TILE; stride <<= 1) {
__syncthreads();
dst[t] = src[t] + ((t >= stride) ? src[t - stride] : 0.0f);
float *tmp = src; src = dst; dst = tmp;
}
__syncthreads(); // 让 src[TILE-1] 可见
if (base + t < n) Y[base + t] = src[t];
if (t == 0) blockTotals[blockIdx.x] = src[TILE - 1];
}
// ===============================================================
// B-2. kernel2：对 blockTotals[0..nb-1] 做 exclusive scan，单 block 完成。
// 结构：每线程先顺序累加 PER_THREAD 个段总和 -> 块内 Kogge-Stone 扫描
// 这 1024 个“线程小计” -> 转 exclusive -> 线程内偏移加回并写出。
// ===============================================================
__global__ void scan_block_totals_kernel(const float *X, float *Y, int nb)
{
__shared__ float T0[BLOCK_SIZE];
__shared__ float T1[BLOCK_SIZE];
const int t = threadIdx.x;
float v[PER_THREAD]; // 每线程私有的 8 个值（寄存器数组）
#pragma unroll
for (int k = 0; k < PER_THREAD; ++k) {
const int idx = t * PER_THREAD + k;
v[k] = (idx < nb) ? X[idx] : 0.0f;
}
float s = 0.0f;
#pragma unroll
for (int k = 0; k < PER_THREAD; ++k) s += v[k];
float *src = T0;
float *dst = T1;
src[t] = s; // 线程小计
for (int stride = 1; stride < BLOCK_SIZE; stride <<= 1) {
__syncthreads();
dst[t] = src[t] + ((t >= stride) ? src[t - stride] : 0.0f);
float *tmp = src; src = dst; dst = tmp;
}
__syncthreads(); // 保证 1024 个前缀和全部写回
const float offset = (t == 0) ? 0.0f : src[t - 1]; // inclusive -> exclusive
float run = offset;
#pragma unroll
for (int k = 0; k < PER_THREAD; ++k) {
const int idx = t * PER_THREAD + k;
if (idx < nb) {
Y[idx] = run;
run += v[k];
}
}
}
// ===============================================================
// B-3. kernel3：把每一段的 exclusive 前缀（blockOffsets[blockIdx.x]）加回该段
// ===============================================================
__global__ void add_block_offsets_kernel(float *Y, const float *blockOffsets, int n)
{
const int g = blockIdx.x * TILE + threadIdx.x;
if (g < n) Y[g] += blockOffsets[blockIdx.x];
}
// ---------------------------------------------------------------
static void cpu_inclusive_scan(const float *x, double *y, int n)
{
double running = 0.0;
for (int i = 0; i < n; ++i) { running += (double)x[i]; y[i] = running; }
}
static bool verify(const float *gpu, const double *ref, int n, const char *tag)
{
int bad = 0;
double worst = 0.0;
for (int i = 0; i < n; ++i) {
double d = fabs((double)gpu[i] - ref[i]) / (1.0 + fabs(ref[i]));
if (d > worst) worst = d;
if (d > 1e-3) ++bad;
}
printf(" [%-22s] n = %-9d 最大相对误差 = %.3e, 不匹配 = %d -> %s\n",
tag, n, worst, bad, (bad == 0) ? "PASS" : "FAIL");
return bad == 0;
}
int main(void)
{
cudaEvent_t e0, e1, e2, e3, e4;
CUDA_CHECK(cudaEventCreate(&e0));
CUDA_CHECK(cudaEventCreate(&e1));
CUDA_CHECK(cudaEventCreate(&e2));
CUDA_CHECK(cudaEventCreate(&e3));
CUDA_CHECK(cudaEventCreate(&e4));
printf("本程序支持的最大 N = %lld (= %d x %d x %d)\n",
MAX_N, TILE, PER_THREAD, BLOCK_SIZE);
// ================= A. Blelloch 单 block 版：正确性 + 计时 =================
{
const int n = BLELLOCH_N;
const size_t bytes = (size_t)n * sizeof(float);
float *h_x = (float *)malloc(bytes);
float *h_y = (float *)malloc(bytes);
double *h_ref = (double *)malloc((size_t)n * sizeof(double));
srand(3);
for (int i = 0; i < n; ++i) h_x[i] = (float)((rand() % 2001) - 1000) / 1000.0f;
double run = 0.0; // exclusive 参考：Y[0]=0, Y[i]=X[0..i-1]
for (int i = 0; i < n; ++i) { h_ref[i] = run; run += (double)h_x[i]; }
float *d_x = NULL, *d_y = NULL;
CUDA_CHECK(cudaMalloc((void **)&d_x, bytes));
CUDA_CHECK(cudaMalloc((void **)&d_y, bytes));
CUDA_CHECK(cudaMemcpy(d_x, h_x, bytes, cudaMemcpyHostToDevice));
printf("\n=== A. Blelloch exclusive scan（单 block, n = %d）===\n", n);
blelloch_exclusive_scan_kernel<<<1, BLOCK_SIZE>>>(d_x, d_y, n);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaMemcpy(h_y, d_y, bytes, cudaMemcpyDeviceToHost));
verify(h_y, h_ref, n, "Blelloch exclusive");
printf(" 步数 = 2*log2(n) = %d (up-sweep %d 步 + down-sweep %d 步)\n",
2 * 11, 11, 11);
printf(" 加法总数 = 2*(n-1) = %d, 即每元素 %.2f 次 (对比 Kogge-Stone 的 %.2f 次)\n",
2 * (n - 1), 2.0 * (n - 1) / n, 9.0);
const int reps = 20000;
blelloch_exclusive_scan_kernel<<<1, BLOCK_SIZE>>>(d_x, d_y, n);
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaEventRecord(e0));
for (int r = 0; r < reps; ++r)
blelloch_exclusive_scan_kernel<<<1, BLOCK_SIZE>>>(d_x, d_y, n);
CUDA_CHECK(cudaEventRecord(e1));
CUDA_CHECK(cudaEventSynchronize(e1));
float ms = 0.0f;
CUDA_CHECK(cudaEventElapsedTime(&ms, e0, e1));
printf(" 重复 %d 次总耗时 %.2f ms, 每次 %.3f us\n", reps, ms, ms * 1000.0 / reps);
free(h_x); free(h_y); free(h_ref);
CUDA_CHECK(cudaFree(d_x));
CUDA_CHECK(cudaFree(d_y));
}
// ================= B. 层次化多 block 版：多种 N 的正确性 =================
const int test_ns[] = {1, 1023, 1024, 1025, 4096, 100000, 1 << 22};
const int n_tests = (int)(sizeof(test_ns) / sizeof(test_ns[0]));
printf("\n=== B. 层次化三步法：正确性 ===");
printf("（三个 kernel：段内扫描 / 段总和扫描 / 加回偏移）\n");
for (int ti = 0; ti < n_tests; ++ti) {
const int n = test_ns[ti];
const int nb = (n + TILE - 1) / TILE;
if (nb > MAX_TOTALS) { printf(" n = %d 超过本程序上限, 跳过\n", n); continue; }
const size_t bytes = (size_t)n * sizeof(float);
const size_t tbytes = (size_t)nb * sizeof(float);
float *h_x = (float *)malloc(bytes);
float *h_y = (float *)malloc(bytes);
double *h_ref = (double *)malloc((size_t)n * sizeof(double));
srand(5 + ti);
for (int i = 0; i < n; ++i) h_x[i] = (float)((rand() % 2001) - 1000) / 1000.0f;
cpu_inclusive_scan(h_x, h_ref, n);
float *d_x = NULL, *d_y = NULL, *d_tot = NULL, *d_off = NULL;
CUDA_CHECK(cudaMalloc((void **)&d_x, bytes));
CUDA_CHECK(cudaMalloc((void **)&d_y, bytes));
CUDA_CHECK(cudaMalloc((void **)&d_tot, tbytes));
CUDA_CHECK(cudaMalloc((void **)&d_off, tbytes));
CUDA_CHECK(cudaMemcpy(d_x, h_x, bytes, cudaMemcpyHostToDevice));
scan_tiles_kernel<<<nb, BLOCK_SIZE>>>(d_x, d_y, d_tot, n);
scan_block_totals_kernel<<<1, BLOCK_SIZE>>>(d_tot, d_off, nb);
add_block_offsets_kernel<<<nb, BLOCK_SIZE>>>(d_y, d_off, n);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_y, d_y, bytes, cudaMemcpyDeviceToHost));
verify(h_y, h_ref, n, "hierarchical inclusive");
free(h_x); free(h_y); free(h_ref);
CUDA_CHECK(cudaFree(d_x));
CUDA_CHECK(cudaFree(d_y));
CUDA_CHECK(cudaFree(d_tot));
CUDA_CHECK(cudaFree(d_off));
}
// ================= C. 三个 kernel 分别计时（N = 4M）=================
{
const int n = 1 << 22; // 4,194,304 个元素 = 16 MB
const int nb = n / TILE; // 4096 个 block
const int reps = 200;
const size_t bytes = (size_t)n * sizeof(float);
const size_t tbytes = (size_t)nb * sizeof(float);
float *h_x = (float *)malloc(bytes);
float *h_y = (float *)malloc(bytes);
double *h_ref = (double *)malloc((size_t)n * sizeof(double));
srand(17);
for (int i = 0; i < n; ++i) h_x[i] = (float)((rand() % 2001) - 1000) / 1000.0f;
auto tc0 = std::chrono::high_resolution_clock::now();
cpu_inclusive_scan(h_x, h_ref, n);
auto tc1 = std::chrono::high_resolution_clock::now();
const double cpu_ms = std::chrono::duration<double, std::milli>(tc1 - tc0).count();
float *d_x = NULL, *d_y = NULL, *d_tot = NULL, *d_off = NULL;
CUDA_CHECK(cudaMalloc((void **)&d_x, bytes));
CUDA_CHECK(cudaMalloc((void **)&d_y, bytes));
CUDA_CHECK(cudaMalloc((void **)&d_tot, tbytes));
CUDA_CHECK(cudaMalloc((void **)&d_off, tbytes));
CUDA_CHECK(cudaMemcpy(d_x, h_x, bytes, cudaMemcpyHostToDevice));
// 预热
scan_tiles_kernel<<<nb, BLOCK_SIZE>>>(d_x, d_y, d_tot, n);
scan_block_totals_kernel<<<1, BLOCK_SIZE>>>(d_tot, d_off, nb);
add_block_offsets_kernel<<<nb, BLOCK_SIZE>>>(d_y, d_off, n);
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaEventRecord(e0));
for (int r = 0; r < reps; ++r)
scan_tiles_kernel<<<nb, BLOCK_SIZE>>>(d_x, d_y, d_tot, n);
CUDA_CHECK(cudaEventRecord(e1));
for (int r = 0; r < reps; ++r)
scan_block_totals_kernel<<<1, BLOCK_SIZE>>>(d_tot, d_off, nb);
CUDA_CHECK(cudaEventRecord(e2));
for (int r = 0; r < reps; ++r)
add_block_offsets_kernel<<<nb, BLOCK_SIZE>>>(d_y, d_off, n);
CUDA_CHECK(cudaEventRecord(e3));
CUDA_CHECK(cudaEventSynchronize(e3));
float k1 = 0.0f, k2 = 0.0f, k3 = 0.0f;
CUDA_CHECK(cudaEventElapsedTime(&k1, e0, e1));
CUDA_CHECK(cudaEventElapsedTime(&k2, e1, e2));
CUDA_CHECK(cudaEventElapsedTime(&k3, e2, e3));
k1 /= (float)reps; k2 /= (float)reps; k3 /= (float)reps;
// 计时循环把 d_y 累加了多次，重新算一遍再校验
CUDA_CHECK(cudaMemcpy(d_x, h_x, bytes, cudaMemcpyHostToDevice));
scan_tiles_kernel<<<nb, BLOCK_SIZE>>>(d_x, d_y, d_tot, n);
scan_block_totals_kernel<<<1, BLOCK_SIZE>>>(d_tot, d_off, nb);
add_block_offsets_kernel<<<nb, BLOCK_SIZE>>>(d_y, d_off, n);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaMemcpy(h_y, d_y, bytes, cudaMemcpyDeviceToHost));
printf("\n=== C. 层次化扫描三个 kernel 分别计时 (N = %d, %.1f MB, %d 个 block) ===\n",
n, (double)bytes / 1048576.0, nb);
const double total_ms = k1 + k2 + k3;
const double moved1 = 2.0 * bytes + tbytes;
const double moved2 = 2.0 * tbytes;
const double moved3 = 2.0 * bytes + tbytes;
const double moved = moved1 + moved2 + moved3;
printf(" kernel1 段内扫描 : %7.3f ms 访存 %.1f MB -> %7.1f GB/s\n",
k1, moved1 / 1048576.0, moved1 / (k1 * 1e-3) / 1e9);
printf(" kernel2 段总和扫描 : %7.3f ms 访存 %.1f MB -> %7.1f GB/s\n",
k2, moved2 / 1048576.0, moved2 / (k2 * 1e-3) / 1e9);
printf(" kernel3 加回偏移 : %7.3f ms 访存 %.1f MB -> %7.1f GB/s\n",
k3, moved3 / 1048576.0, moved3 / (k3 * 1e-3) / 1e9);
printf(" ------------------------------------------------------------\n");
printf(" 合计 : %7.3f ms 总访存 %.1f MB -> %7.1f GB/s\n",
total_ms, moved / 1048576.0, moved / (total_ms * 1e-3) / 1e9);
printf(" 纯带宽下限 (2N*4B/1555GBs) : %.3f ms -> 层次化用掉 %.2f 倍流量\n",
2.0 * bytes / 1555e9 * 1e3, moved / (2.0 * bytes));
printf(" CPU 单线程串行扫描 : %7.3f ms -> GPU 加速比 %.1fx\n",
cpu_ms, cpu_ms / total_ms);
verify(h_y, h_ref, n, "hierarchical inclusive");
printf(" 达到的最优有效带宽 : %.1f GB/s = 峰值的 %.1f%%\n",
2.0 * bytes / (total_ms * 1e-3) / 1e9,
2.0 * bytes / (total_ms * 1e-3) / 1555e9 * 100.0);
free(h_x); free(h_y); free(h_ref);
CUDA_CHECK(cudaFree(d_x));
CUDA_CHECK(cudaFree(d_y));
CUDA_CHECK(cudaFree(d_tot));
CUDA_CHECK(cudaFree(d_off));
}
CUDA_CHECK(cudaEventDestroy(e0));
CUDA_CHECK(cudaEventDestroy(e1));
CUDA_CHECK(cudaEventDestroy(e2));
CUDA_CHECK(cudaEventDestroy(e3));
CUDA_CHECK(cudaEventDestroy(e4));
return 0;
}
- 【代码做什么？】
- 打印本程序支持的最大
N = TILE × PER_THREAD × BLOCK_SIZE = 1024 × 8 × 1024 = 8,388,608（受 kernel2「一个 block 扫完所有段总和」的能力限制）。 - A 部分（Blelloch）：
n = 2048，先在主机上用double生成 exclusive scan 参考（Y[0] = 0，Y[i]等于X[0]到X[i-1]的元素之和），启动blelloch_exclusive_scan_kernel<<<1, 1024>>>校验；打印「步数 = 2·log2(n) = 22」与「加法总数 = 2(n−1) = 4094，每元素 2.0 次」；再重复 20000 次测量单次启动时间。 - B 部分（层次化正确性）：对
N ∈ {1, 1023, 1024, 1025, 4096, 100000, 4194304}七个规模各跑一遍三步法并校验，覆盖「不足一个 block」「恰好一个 block」「跨 block 边界」「不是 TILE 整倍数」「超大数组」这些边界情况。 - C 部分（分层计时）：
N = 2²²，用四组cudaEvent把 kernel1、kernel2、kernel3 分别计时 200 次（e0→e1是 kernel1、e1→e2是 kernel2、e2→e3是 kernel3），每个 kernel 换算出「搬运字节数 / 时间」的有效带宽。 - 因为 kernel3 是累加（
Y[g] += offset）、不幂等，重复计时会污染结果，所以计时结束后重新完整跑一遍三步再做最终校验，并计算总时间、总流量、纯带宽下限、相对 CPU 串行扫描的加速比。 - kernel 内部：kernel1 就是示例 2 的块内扫描，外加
if (t == 0) blockTotals[blockIdx.x] = src[TILE-1];；kernel2 每个线程先顺序累加 8 个段总和（寄存器数组v[8]，用#pragma unroll保证落在寄存器里），再用共享内存双缓冲对 1024 个「线程小计」做 Kogge-Stone 扫描，转成 exclusive 后把偏移加回每个段总和并写出；kernel3 让每个元素加上blockOffsets[blockIdx.x]。
- 打印本程序支持的最大
- 【并行机制与硬件映射解说】
- kernel1 / kernel3：
grid = 4096、block = 1024线程 = 32 warp；A100 每 SM 驻留 2 个 block → 216 个 block 同时在机，4096 个 block 需要4096/216 ≈ 19个波次（wave）。每波末尾都会因为 block 结束/新的 block 上线而损失一点效率，这是尾部效应。 - kernel2 的极端不对称：整个 GPU 只有 1 个 block 在跑，也就是 216 个 block 槽位里只用了 1 个（约 0.5% 的 SM 资源），却要花掉接近 4 µs。这是层次化算法固有的代价：中间数据只有 16 KB，但它的延迟必须由整个 kernel 的生命周期来支付。kernel2 内部还有一段「每线程顺序处理 8 个元素」的串行段（8 次依赖加法 = 约 32 cycles），以及 10 次
__syncthreads()。 - 寄存器 vs local memory：
float v[PER_THREAD]如果不用#pragma unroll展开，编译器可能把它放进 local memory（本质在显存里，只是有 L1/L2 缓存），每次读写都要走 cache，性能会掉一个数量级。加上#pragma unroll后，8 个 float 全部分配到寄存器（kernel2 总寄存器数约 32 个/线程，65536/(32×1024) = 2个 block/SM，不构成限制）。 - 全局内存访问：kernel3 里同一个 block 的 1024 个线程读的是同一个
blockOffsets[blockIdx.x]，这是「同一地址访问」，硬件广播、不算 bank conflict（虽然这里它其实走的是 L1 常量路径）；对Y的读改写是连续地址，32 线程 × 4 B = 128 B = 一次 transaction，完全合并。 - bank conflict：kernel1 是 Hillis-Steele，无冲突（表 A）；kernel2 也是 Hillis-Steele 结构，同样无冲突；只有 Blelloch 版本（A 部分）会遇到树形索引带来的 2～32 路冲突（表 B）。
- 数值特性：并行扫描的加法结合顺序与串行扫描不同，浮点加法不满足结合律，所以结果不会与串行参考逐位相同。校验用的是归一化相对误差（阈值 1e-3），而不是要求逐位相等——这是并行归约/扫描类程序的标准做法。
- kernel1 / kernel3：
- 【性能优化分析】
三个程序的性能对照表（A100 sm_80；数值由带宽/延迟模型按下列公式推算：带宽下限
流量 / 1555 GB/s，ALU 时间加法次数 / 9.75e12，单 block kernel 时间屏障次数 × 30 cycles + 全局往返延迟 2 × 600 cycles）：程序 / 版本 N 工作量（加法次数） 时间 有效带宽 相对 CPU 串行 示例 1 CPU 单线程串行（double） 65,536 6.6e4 0.09 ms — 1.0× 示例 1 GPU 朴素 O(N²) 65,536 2.15e9 0.62 ms — 0.14×（反而慢 7 倍） 示例 2 单 block KS（grid = 1） 1,024 9,217 3.0 µs（其中启动开销 2–3 µs） — — 示例 2 块内扫描阶段（grid = 4096） 4,194,304 3.78e7 27.0 µs 1240 GB/s — 示例 3 Blelloch 单 block 2,048 4,094 3.4 µs（其中启动开销 2–3 µs） — — 示例 3 kernel1 段内扫描 4,194,304 3.78e7 30.0 µs 1120 GB/s — 示例 3 kernel2 段总和扫描 4,096 1.7e4 3.6 µs 9 GB/s（延迟受限） — 示例 3 kernel3 加回偏移 4,194,304 4.19e6 24.0 µs 1400 GB/s — 示例 3 三层合计 4,194,304 4.20e7 57.6 µs 1165 GB/s 87× 理论下限（2N×4 B / 带宽） 4,194,304 — 21.6 µs 1555 GB/s 233× - ARITHMETIC INTENSITY（算术强度）与 Roofline：
- kernel1：
9.0 次加法/元素 ÷ 8 字节/元素 = 1.125次加法/字节； - kernel3：
1 次加法/元素 ÷ 8 字节 = 0.125次加法/字节； - Blelloch：
2.0 ÷ 8 = 0.25次加法/字节。 A100 的脊点是9.75e12 / 1555e9 = 6.27次加法/字节，三者分别低 5.6 倍、50 倍、25 倍，全部牢牢贴在带宽墙上。
- kernel1：
- work efficiency 与 latency 的权衡（本讲最重要的结论）：
- Hillis-Steele：
log2(N)步、O(N log N)工作； - Blelloch：
2·log2(N)步、O(N)工作。 - 在 CPU 上「工作量」就是时间，所以 Blelloch 更优；在 GPU 上时间由 max(带宽时间, ALU 时间) 决定。以 N = 4M 的全数组扫描为例：Hillis-Steele 需要
4.2e7次加法 →4.2e7 / 9.75e12 = 4.3 µs；而搬 67 MB 需要43 µs。ALU 只占访存时间的 10%，log N 倍的冗余加法被完全藏在内存延迟后面。 - 反过来，Blelloch 的代价是：步数翻倍（屏障翻倍）、后半段并行度骤降（up-sweep 最后 6 步总共只有 63 个线程干活，
63/2048 = 3%的并行度）、并且树形索引带来 2～32 路 bank conflict（表 B：事务数膨胀 6 倍）。所以在 GPU 的块内扫描里，Kogge-Stone 反而更受欢迎（讲义原话：”Kogge-Stone is more popular for parallel scan with blocks in GPUs”）。 - 结论：work efficiency 是 CPU 时代的首要指标；在 GPU 上，当 ALU 有 5–10 倍余量时，应当优先优化延迟、屏障次数与并行度，而不是死抠操作数。
- Hillis-Steele：
- 扫描的带宽下限：一次扫描至少要「读一遍输入、写一遍输出」，所以理论下限 =
2N × 4 B / 带宽。N = 4M 时 =33.55 MB / 1555 GB/s = 21.6 µs。层次化三步法的实际流量是4N × 4 B = 67.1 MB（中间结果必须写回全局内存再读回来），恰好是下限的 2 倍，因此它的理论下限是 43.2 µs；实际 57.6 µs 中，多出来的部分来自 kernel2 的 3.6 µs 纯延迟开销、每次 kernel 启动的 2–3 µs，以及有效带宽只跑到 75%（1165/1555）。 - Amdahl 视角：kernel2 只搬运 32 KB（占总流量的 0.05%）却花掉 3.6 µs（占总时间的 6.3%）。即使把 kernel1 与 kernel3 优化到完美带宽（21.6 + 21.6 = 43.2 µs），总时间也只能降到约 46.8 µs。要真正突破，必须减少 kernel 数量（例如用
atomicAdd+ 全局标志实现单 kernel 的 decoupled look-back 链式扫描），或用 CUB 的DeviceScan这类高度优化且向量化的实现。 - 瓶颈判定与优化方向：读写受限（bandwidth-bound）。可执行优化：① 向量化（每个线程用
float4处理 4 个元素），把在途数据从108×2048×4 B = 884 KB提到 3.5 MB，稳住带宽；② 用 grid-stride 循环减少 block 数、提高每线程工作量；③ 用__shfl_up_sync把每 block 的前 5 级扫描放进寄存器，屏障从 10 次降到 5 次；④ kernel2 改用 1024 个线程 × 更多元素或用多个 block 的二级层次化，但要注意此时又要多一个 kernel；⑤ 用共享内存作为「段内扫描」的暂存、但让 kernel3 与 kernel1 融合（把Y_local留在共享内存里、只把段总和写出去）可以省掉一次 16 MB 的写和一次 16 MB 的读——代价是 kernel1 无法结束、无法跨块通信，所以这一步必须靠 cooperative groups 的 grid 同步或单 kernel 链式扫描来实现。
性能优化技巧总结
- 先换算法、再调微架构：O(N²) 的朴素并行版比串行还慢 7 倍，而 Hillis-Steele 把工作量降到 O(N log N)、Blelloch 降到 O(N)。数量级的收益只能来自算法，不能来自调 block 大小。
- 块内小规模（n ≤ 1024）用 Hillis-Steele：步数只有 log2(n) = 10，虽然多做了 9 倍加法，但 ALU 余量有 5.6 倍，冗余工作完全被内存延迟掩盖；而且它的访问模式是连续地址，共享内存零 bank conflict。
- 用双缓冲（ping-pong）把每轮的屏障从 2 次减到 1 次：读写落在不同的共享数组上，同一轮内不存在 RAW/WAR 依赖；代价只是每 block 的共享内存翻倍（8 KB → 16 KB；2 个 block/SM 合计 32 KB，只占 A100 每 SM 164 KB 的约 20%）。
- 别忘了「空转线程也要抄一份」：
t < stride的线程必须把值复制到输出缓冲（用算子单位元写作src[t] + 0），否则缓冲角色互换后数据就丢了。 - 每线程处理多个元素并用
float4向量化：这是带宽受限扫描最有效的一招——它同时提高内存级并行度（在途字节 884 KB → 3.5 MB）、减少地址算术、减少屏障次数（每 4 个元素一次同步）。 - warp 内扫描交给
__shfl_up_sync：前 5 级（stride = 1、2、4、8、16）只涉及同一个 warp 的线程，用 shuffle 在寄存器间传递，既不碰共享内存也省掉 5 次__syncthreads()。 - block 取 1024 线程：A100 上 2 个 block/SM 刚好等于 2048 线程上限，达到 100% 占用率；同时 Hillis-Steele 的步数只有 10 步，屏障总数最少。
- 保证合并访问：block 的段与全局内存的连续区间对齐，让一个 warp 的 32 个线程读 128 个连续字节（一条 cache line、一次 transaction）。
- 层次化时让 kernel2 尽量小：段总和的个数是
N/1024，把它压到能用一个 block 扫完（或干脆用 decoupled look-back 单 kernel）。中间层的延迟开销（约 3.6 µs）在总时间里的占比比它的数据量占比高两个数量级。 - 计时要用
cudaEvent、先预热、并对非幂等 kernel 重跑：kernel 启动开销 2–3 µs 会淹没小 kernel；Y[g] += blockOffsets[blockIdx.x]这种累加型 kernel 重复执行会污染结果。 - 校验用独立参考实现：CPU 串行（必要时用
double累加）+ 归一化相对误差阈值；并行扫描的加法结合顺序与串行不同，逐位相等不是合理要求。
关键要点
- 扫描 = 归约 + 分发：归约只保留最后一个结果（是扫描的简化形式），扫描必须保留并分发所有中间结果；两者都要求算子满足结合律（扫描不要求交换律）。
- Hillis-Steele（Kogge-Stone）：
log2(N)步、O(N log N)工作量、非 work-efficient、低延迟、访问模式天然无 bank conflict，是 GPU 块内扫描的首选。 - Blelloch（up-sweep + down-sweep）：
2·log2(N)步、O(N)工作量、输出 exclusive scan，但并行度前后不均、树形索引有 2～32 路 bank conflict。 - 双缓冲是消除原地读写竞争的标准手段：读
src、写dst、每轮交换指针，屏障数减半；代价是共享内存翻倍。 - 扫描是带宽受限问题：理论下限
2N × 4 B / 带宽（A100 上 N = 4M 时为 21.6 µs）；层次化三步法要搬4N × 4 B，是下限的 2 倍。 - work efficiency 与 latency 的权衡：GPU 的 ALU 相对带宽有 5–10 倍余量，因此「多做 log N 倍加法但少一半屏障、并行度高一半」的 Hillis-Steele 往往比「工作高效但要 2 倍步数」的 Blelloch 更快——在 GPU 上，带宽和延迟比操作数更重要。
- 块间通信只能靠全局内存 + kernel 边界：寄存器与共享内存对其他 block 不可见，因此任意长度 N 的扫描必须用「段内扫描 → 段总和扫描 → 偏移加回」的多 kernel 层次化结构。
常见陷阱与注意事项
- 忘记
__syncthreads()：共享内存里每个 stride 轮的读都依赖上一轮的写，少一道屏障就会读到半成品，结果随运行次数变化 → 每个 stride 轮开头都要__syncthreads()，并且它必须放在所有线程都会执行的位置，绝不能写进if (t < stride) { __syncthreads(); }这类分支里（会导致死锁或未定义行为）。 - 原地版本只放一道
__syncthreads()：if (t >= stride) T[t] += T[t-stride];在同一轮内既有 RAW（读到别人本轮的新值）又有 WAR（自己覆盖了别人还要读的旧值）竞争，讲义明确指出 “This code has a data race condition” → 要么用双缓冲，要么「先读到寄存器、第二道屏障之后再写回」。 - 双缓冲忘记交换角色或忘记复制空转线程的值：结果整体错位（例如
y[0]莫名其妙变成某个中间值）。→ 循环末尾一定要交换src/dst指针，且t < stride的线程也要往dst[t]写（写src[t] + 单位元），循环结束后结果在交换后的src里。 - 忘记
cudaDeviceSynchronize()/cudaEventSynchronize()：cudaEventElapsedTime在 kernel 还没跑完时读到的是「启动完成」的时间，测出来的数字比真实值小一个数量级；kernel 报错也只有同步时才会暴露。→ 计时后必须cudaEventSynchronize，调试时用cudaDeviceSynchronize()配合cudaGetLastError()。 - 不检查 CUDA API 返回值：
cudaMalloc/cudaMemcpy失败只返回错误码，程序会带着空指针继续跑并在后面莫名其妙崩掉。→ 统一用CUDA_CHECK宏包裹每一次 CUDA 调用；kernel 启动后用cudaGetLastError()检查。 cudaMemcpy方向写反：cudaMemcpyHostToDevice/cudaMemcpyDeviceToHost弄反会让结果全是 0 或旧数据，而程序不报错（同尺寸拷贝是合法的）。→ 记住「第一个参数是目的地」，拷贝方向永远指向第一个参数。- 网格覆盖不足 / 越界访问：
grid = N / blockDim的整数除法向下取整，N不是块大小整数倍时尾部元素没人算；反过来grid = (N + block - 1)/block之后 kernel 里必须if (g < n)保护，共享内存里越界的部分要填算子的单位元（加法填 0、乘法填 1、min填 +∞），否则补进来的垃圾值会污染前缀和。→ 网格用向上取整公式，kernel 内双重保护（全局下标判断 + 单位元填充）。 - 共享内存容量超限：双缓冲 HL 版每 block 需要
2 × blockDim × 4 B；block = 1024 时是 8 KB（安全），但如果每线程再放一个float v[8]的共享数组就会涨到 40 KB，占用率掉到 4 个 block/SM 以下。→ 用cudaFuncSetAttribute提高动态共享内存上限前，先用 occupancy calculator 算清楚；每线程的私有数组应放在寄存器里。 - 每线程的私有数组没被展开而落到 local memory：
float v[8]若用运行期下标访问，编译器会把它放到 local memory（物理上在显存，读写延迟与全局内存同量级），性能掉一个数量级。→ 用#pragma unroll+ 编译期常量下标，让它留在寄存器。 - host / device 指针混用：把
h_x（malloc的指针）直接传给 kernel，或者把d_x传给printf/fwrite，都会触发非法访问或读到垃圾。→ 命名上区分h_/d_前缀，所有 kernel 参数必须是cudaMalloc/cudaMallocManaged得到的指针。 - 忽略 bank conflict 的适用对象：Hillis-Steele 的连续访问没有冲突（不要凭「stride=1 必冲突」的传言去盲目加 padding），而 Blelloch / Brent-Kung 的树形索引有 2～32 路冲突（表 B）。→ 用「32 个连续字一定覆盖 32 个不同 bank」「步长 k 的冲突度 = 32/gcd(k,32)」这两条规则实际算一遍，再决定要不要改访问模式。
- 计时循环里执行非幂等 kernel：
Y[g] += blockOffsets[blockIdx.x]这类累加 kernel 重复 200 次会把偏移量加 200 遍，校验必然FAIL。→ 计时与校验分开：计时用一次性缓冲或直接对输入可重复的 kernel 计时，校验前重新完整跑一遍。
思考题（带答案）
Q1. 给定 N = 4,194,304 个 float（16 MB），在 A100 上分别用 Hillis-Steele 与 Blelloch 做全数组扫描。Hillis-Steele 要多做约 10 倍的加法，为什么它反而常常更快？请用至少两个具体数字论证。
答案：因为扫描是带宽受限的。① 算术强度：Hillis-Steele 每元素约 9.0 次加法、8 字节 I/O，算术强度 1.125 次加法/字节，而 A100 的 Roofline 脊点是 9.75e12 ÷ 1555e9 = 6.27 次加法/字节，它比脊点低 5.6 倍，多做的工作落在带宽墙左侧、被内存延迟掩盖。② 时间对照：Hillis-Steele 全数组的加法总量 4096 block × 9217 = 3.78e7 次，纯 ALU 时间 3.78e7 ÷ 9.75e12 = 3.9 µs；而搬 32 MB 数据需要 21.6 µs，ALU 只占 18%。③ 步数与并行度：Hillis-Steele 只要 log2(1024) = 10 次 __syncthreads()，而 Blelloch 需要 2 × 11 = 22 次；Blelloch 的 up-sweep 后半段活跃线程只有 64、32、16、8、4、2、1，最后 6 步总共只有 63 个线程干活（并行度 3%），还会因为树形索引产生最高 32 路的 bank conflict（事务数膨胀 6 倍）。所以在 GPU 的块内扫描里应优先选 Hillis-Steele，Blelloch 的价值在于「当并行资源不足或需要 exclusive 输出时把工作量压到 O(N)」。
Q2. Blelloch 扫描在 n = 2048、BLOCK = 1024 线程、float 数组上运行时，up-sweep 中 stride = 16 的那一步：warp 内有几个线程活跃？冲突度是多少？这一步的共享内存事务数是多少？如果 stride = 16 换成 Hillis-Steele 的访问模式又会怎样？
答案：stride = 16 时，活跃线程满足 index = (t+1)·2·16 − 1 < 2048，即 t ≤ 63，64 个线程活跃（2 个 warp）。相邻线程的地址差是 2 × 16 = 32 个字，而 bank 号是「字地址 mod 32」，所以这 32 个线程全部落在同一个 bank 的不同地址上 → 冲突度 32（最坏情况）。这一步每次访问需要 32 个 warp-访问串行化：2 个 warp × 32 次 = 64 次共享内存事务，而无冲突时只需 2 次（每 warp 一次）。换成 Hillis-Steele 的 src[t] + src[t−16] 呢？同一个 warp 的 32 个线程访问的是 32 个连续的字，必然覆盖 32 个不同 bank，冲突度为 1——但代价是活跃线程数变成 1024 − 16 = 1008（31.5 个 warp，几乎满并行），一步就要 63 次事务。两者对照说明：Blelloch 用「并行度」换了「工作量」，而 bank conflict 又把它换来的好处吃回去一部分。（顺带一个相关结论：Brent-Kung 处理 1024 个元素时 stride 依次是 512、256、128、64、32、16，前 5 步每个 warp 要么全活跃要么全不活跃，第 6 步（stride = 16）才第一次出现 warp 内分歧。）
Q3. 为什么任意长度 N 的扫描必须用多个 kernel（或单 kernel 的链式扫描）？层次化三步法的 DRAM 流量是理论下限的几倍？在 N = 4M、A100 上分别是多少微秒？
答案：因为一个 block 的寄存器与共享内存对其他 block 不可见，而全局内存的写必须经过内存栅栏（memory fence）才对其他 block 可见——在 CUDA 里最自然的栅栏就是 kernel 结束。所以「段内扫描 → 段总和扫描 → 偏移加回」这三步之间的数据必须写到全局内存，并用 kernel 边界（或 cooperative groups 的 grid 同步、或带内存序的原子标志）来保证可见性。流量方面：最低需求是「读一遍输入 + 写一遍输出」= 2N × 4 B；层次化三步法要写一次段内结果、再读一次段内结果（再加回偏移），所以是 4N × 4 B，恰好是下限的 2 倍。代入 A100 的 1555 GB/s：下限 33.55 MB ÷ 1555 GB/s = 21.6 µs，层次化下限 67.1 MB ÷ 1555 GB/s = 43.2 µs；实测（模型推算）三个 kernel 分别是 30.0 + 3.6 + 24.0 = 57.6 µs，相对单线程 CPU 串行扫描（约 5.0 ms）加速 87 倍，相对 21.6 µs 的绝对下限还有 2.7 倍差距——差距来自 kernel2 的纯延迟（3.6 µs）与每次 kernel 启动的 2–3 µs 开销。