Lecture 8: 数据并行思维（Data-Parallel Thinking）（日期：Oct 16, 2025）

目录 · ← l7 · l9 →

Lecture 8: 数据并行思维（Data-Parallel Thinking）（日期：Oct 16, 2025）

概述：前几讲我们习惯从”worker 做什么、如何把工作分配给 worker”的角度思考并行编程；本讲换一个视角：把算法描述成对数据序列（sequence）的操作——map、filter、fold/reduce、scan/segmented scan、sort、groupBy、gather/scatter 等。核心思想是：这些操作的高性能并行实现已经存在，因此用这些原语写的程序往往能高效跑在并行机器上（前提：不要被内存带宽卡住）。本讲还会深入 scan（前缀和）的并行算法（O(N log N) 朴素算法 vs O(N) 工作高效算法）及其多级实现，并用稀疏矩阵乘法、粒子网格、直方图等实例展示如何把”不规则并行”转化为”规则并行”、把”细粒度同步”转化为”粗粒度同步”。

注意：本讲与 Assignment 3 直接相关——幻灯片明确说 Assignment 3 提供了与本讲 scan_block 类似的 CUDA scan 代码；此外本讲的数据并行思维（map/reduce/scan）也是后续 Assignment 5（最快速 CUDA kernel）的思维基础。

一、核心概念与定义

1. 数据并行模型（Data-parallel model）与序列（Sequence）

• 定义：把计算组织成对元素序列的操作（例如：对序列的所有元素执行同一个函数）。序列（sequence）是有序的元素集合（C++ 的 Sequence<T>、Scala 的 List[T]、Pandas Dataframe、PyTorch/JAX Tensor、Haskell 的 seq T）。关键：与数组不同，程序只能通过特定操作访问序列元素，不能直接按下标访问——这给了实现方重排/并行化的自由。现代著名例子：NumPy 的 C = A + B。
• 现实类比：序列像一条流水线上的”待加工件队列”，工人（worker）只允许通过规定的操作（加工台）触碰工件，不允许伸手到队列中间乱拿——这样流水线才能自由调度。
• 公式/图示：C = A + B（三个等长向量逐元素相加）就是一次 map 级联。

2. Map（映射）

• 定义：高阶函数（以函数为参数的函数）。把无副作用（side-effect free）的一元函数 f :: a -> b 应用到输入序列的每个元素上，产生等长的输出序列。Haskell：map :: (a -> b) -> seq a -> seq b；C++：std::transform；JAX：vmap。
• 现实类比：给一摞文件（输入序列）每份盖章（函数 f），得到一摞盖好章的文件（输出序列）——每份文件互不影响，盖章顺序随便。
• 公式/图示：

  a = [3, 8, 4, 6, 3, 9, 2, 8]      f(x) = x + 10
  b = map f a
  b = [13, 18, 14, 16, 13, 19, 12, 18]   ← 每个元素独立应用 f，顺序可任意
  

3. Filter（过滤）

• 定义：删除序列中不满足谓词（predicate）的元素。输出是输入的子序列（长度 ≤ 输入）。
• 现实类比：安检——只放行”合格”的行李，不合格的留在外面；通过顺序无关紧要。
• 公式/图示：filter f s，其中 f :: a -> Bool。例：过滤掉奇数 → [3,8,4,6,3,9,2,8] → [8,4,6,2,8]。

4. Fold / Reduce（折叠 / 归约）

• 定义：把二元操作 f :: (b,a) -> b 应用到每个元素和一个累加值上，最终把整个序列”折叠”成一个值。fold left：fold :: b -> ((b,a) -> b) -> seq a -> b，初始值（seed）类型为 b。并行 fold 还需要一个额外的二元 combiner 函数 comb :: (b,b) -> b（把子结果合并）；如果 f :: (b,b) -> b 本身就是结合律二元操作，就不需要 combiner。初始值必须是 f 和 comb 的单位元（identity）。
• 现实类比：合唱团报数——每个人把”前一个人报的数 + 自己的数”传下去（串行 fold）；并行版则是每个小组先各自求和，组长再把各组的和加起来（combiner）。
• 公式/图示：

  串行 fold：  fold 10 (+) [3,8,4,6,3,9,2,8] = ((((((10+3)+8)+4)+6)+3)+9)+2)+8 = 53
  并行 fold：  把序列分成若干段，每段并行 fold，再用 comb 合并各段结果：
              sum = comb(comb(seg0, seg1), comb(seg2, seg3))
  

5. Scan / 前缀和（Prefix Sum）

• 定义：给定结合律二元操作 ⊕，inclusive scan 输出 [a0, a0⊕a1, a0⊕a1⊕a2, ...]，即每个输出元素是”从开头到当前位置（含）的所有元素的累计结果”；exclusive scan 输出 [I, a0, a0⊕a1, ...]（不含当前元素，第一个元素是单位元 I）。当 ⊕ = + 时称为 prefix sum（前缀和）。
• 现实类比：电影院排队——每个人想知道”我前面有多少人”。如果每个人只问前面那个人”你前面有多少人”，队伍排多长就得等多久（串行 scan）；并行 scan 让一小群人各自数完自己的小组后，再用”前面所有小组的总人数”一次性修正每个人的号（工作高效并行 scan）。
• 公式/图示：

  in   = [3, 8, 4, 6, 3, 9, 2, 8]
  scan_inclusive(+) = [3, 11, 15, 21, 24, 33, 35, 43]
  scan_exclusive(+) = [0, 3, 11, 15, 21, 24, 33, 35]
  

6. 并行 scan 的两个算法：Work 与 Span

• 定义：Work = 算法执行的总操作数；Span = 最长串行依赖链长度（关键路径）。① 朴素并行 scan（Hillis-Steele 风格，每步跨距翻倍）：Work = O(N log N)（比串行算法还低效！），Span = O(log N)。② 工作高效（work-efficient）scan（Blelloch 算法）：up-sweep（自底向上建树）+ down-sweep（自顶向下修正） 两阶段，Work = O(N)，Span = O(log N)。幻灯片特别提醒：要注意常数因子（”but what is the constant?”）。
• 现实类比：Work 像”总工时”，Span 像”最短工期”。一个算法可以工时很高但工期很短（全员加班并行），也可以工时最优但工期较长。GPU 上有时故意选”低效”算法，因为 SIMD 利用率更高（见下）。
• 公式/图示（Blelloch up-sweep + down-sweep 骨架，⊕ = +）：

  Up-sweep（自底向上，d 从 0 到 log2n-1）：
    forall k 步长为 2^(d+1)：
      a[k + 2^(d+1) - 1] = a[k + 2^d - 1] + a[k + 2^(d+1) - 1]
  Down-sweep（自顶向下，d 从 log2n-1 到 0）：
    x[n-1] = 0
    forall k 步长为 2^(d+1)：
      tmp = a[k + 2^d - 1]
      a[k + 2^d - 1] = a[k + 2^(d+1) - 1]
      a[k + 2^(d+1) - 1] = tmp + a[k + 2^(d+1) - 1]
  

• 补充：两处理器（共享内存）实现——”分而治之 + 加基数”：幻灯片第 18-19 页展示了一个更贴近真实机器的做法：P1 串行 scan 前半段 [a0..a7]，P2 串行 scan 后半段 [a8..a15]（两段完全独立，可并行）；然后 P2 把前半段的总和 base = a0-7 加到后半段的每个元素上。分析：

  P1: 串行 scan a0..a7          （8 次加法）
  P2: 串行 scan a8..a15         （8 次加法，与 P1 并行）
  P2: 后半段每个元素加 base     （8 次加法，只需 base 一个通信值）
  Work = O(N)，常数只有 1.5（2N 次加法中的一半与另一半并行 + N/2 次修正）
  数据访问：空间局部性极好（连续内存）；P2 读 base 的开销在大规模
            NUMA 系统上可能更贵，但小型多核系统上几乎可忽略
  

  这个例子说明：scan 的”最优”实现强烈依赖目标机器的规模与内存层级——两核机器用”串行一半 + 修正一半”；GPU 用”warp 内 SIMD scan + 块内协作”；大规模并行机器才值得用完整的两阶段树形算法。

7. Segmented Scan（分段扫描）

• 定义：scan 的推广：对输入序列的连续分区同时各自做 scan（例如对 [[1,2],[6],[1,2,3,4]] 做分段 exclusive scan 得 [[0,1],[0],[0,1,3,6]]）。常用 start-flag 表示法：一个 flag 序列标记每个分段的起点，与数据序列平行存放。
• 现实类比：一排人分成几组，每组内部各自报数——”组内前面的人总数”，而不是所有人混在一起报数。
• 公式/图示：

  嵌套序列 A = [[1,2,3],[4,5,6,7,8]]
  flag: 1 0 0 \| 1 0 0 0 0
  data: 1 2 3 \| 4 5 6 7 8
  

8. Gather / Scatter（聚集 / 散开）

• 定义：gather(index, input, output)：output[i] = input[index[i]]（按索引序列取数据，可以并行、无冲突）；scatter(index, input, output)：output[index[i]] = input[i]（把数据按索引序列放到目标位置，可能冲突）。硬件支持：AVX2（2013）支持 SIMD gather 但不支持 scatter；AVX-512 有 scatter 指令；GPU 硬件支持 gather/scatter——但都比连续向量的 load/store 贵得多。
• 现实类比：gather 像”按购物清单（index）从货架（input）取货”，每件货独立可取；scatter 像”把货按收货地址（index）投递”，两个货可能抢同一个地址（冲突）。
• 公式/图示：

  gather:  output = [data[3], data[12], data[4], data[9], ...]
  scatter: output[index[i]] = input[i]     ← 多个 i 可能写同一位置 → 需要原子/排序
  

9. GroupByKey（按键分组）与 Sort

• 定义：groupByKey：Seq (key, T) -> Seq (key, Seq T)，把相同 key 的元素聚成”序列的序列”。Sort：按 key 排序整个序列。两者是构建直方图、稀疏矩阵、粒子网格等结构的基础。
• 现实类比：把一摞混着各种颜色的纸牌按花色分堆（groupBy），或按点数排成一条龙（sort）。
• 公式/图示：

  (1,3),(2,8),(2,4),(1,6),(3,3),(1,9),(1,2),(2,8)
      groupByKey ──► (1,[3,6,9,2]), (2,[8,4,8]), (3,[3])
  

10. Work / Span（工作量 / 跨度）与并行度

• 定义：见第 6 点。本讲反复用 Work/Span 分析算法：scan 的”理论并行度与元素数成线性关系”，但实践上”只使用填满机器执行资源所需的并行度即可”（避免过度并行带来的额外通信与同步开销）。
• 现实类比：Work = 总工作量，Span = 关键路径。想快，两个都要小；但 GPU 上有时宁可 Work 大一点也要让 SIMD 通道全部忙起来。
• 公式/图示：并行度上界 ≈ Work / Span（对于理想机器）。

11. 带宽受限（Bandwidth-bound）

• 定义：如果程序的总执行时间由”需要搬多少数据 / 机器带宽”决定，而不是由计算量决定，则称带宽受限。数据并行方案通常需要多趟遍历数据（map → sort → scan → 取末元素……），每一趟都是全量带宽消耗——”if you can avoid being bandwidth bound”是本讲原语的潜台词。
• 现实类比：数据就是”料”，带宽就是”运料的卡车”。算法再精巧，卡车不够，工地就得干等。
• 公式/图示：T ≥ max(总操作数 / 峰值算力, 总字节数 / 带宽)（roofline 思想，第 9 讲展开）。

12. 数据并行思维的核心方法论

• 定义：① 把不规则并行转化为规则并行（如用 sort+scan 替代锁）；② 把细粒度同步转化为粗粒度同步（如每线程一个原子操作 → 块内归约 + 少量原子）；③ 代价：多次遍历数据、额外带宽与存储。
• 现实类比：与其让几千人同时抢一支笔签到（细粒度锁），不如先按姓名拼音排序，再数每个首字母有多少人（sort + scan）——大家都只读不抢，最后合并一次。
• 公式/图示：现代大数据/并行系统的基石：CUDA Thrust、Pandas Dataframe 操作、JAX、Apache Spark / Hadoop 都以这些原语为核心。

案例：把 100 万粒子放进 16 个格子的五种解法（幻灯片第 40-46 页）

问题：按 2D 位置把 1M 个粒子放入 16 格均匀网格（构建”二维数组的列表”），供 N-body 近邻查询（只需查周围格子）。五条路线的权衡光谱，就是本讲方法论的全部：

解法 5 的三步（对应概念 12 的”sort 聚拢 + 差分定位”）：

Step 1 (map):    每粒子算所在格子号 grid_cell[i]
Step 2 (sort):   按格子号排序（粒子索引数组随排序置换）
Step 3 (差分):   每线程比较 grid_cell[i] 与 grid_cell[i-1]，
                 不同即写入 cell_starts[this_cell] 与 cell_ends[prev_cell]
                 （首尾元素特判），最终每个格子得到 [start, end) 区间

结论：解法 5 保持大并行度、完全消除细粒度同步，用一次 sort 和额外遍历（带宽）换来了 GPU 上可行的实现——”This solution maintains a large amount of parallelism and removes the need for fine-grained synchronization… at the cost of a sort and extra passes over the data (extra BW)!”

二、代码示例与详细解说

示例 1：Map——CUDA 内核与 C++ std::transform 的对照

// map.cu —— 编译：nvcc -o map map.cu
// 功能：b[i] = f(a[i])，其中 f(x) = x + 10（map 原语的一次实现）
#include <cuda_runtime.h>

__global__ void mapAdd10(int N, const int* a, int* b)
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < N)                       // 边界保护：N 不必是 block 大小的整数倍
        b[i] = a[i] + 10;            // 无副作用：只读 a[i]，只写 b[i]
}

// host 端：mapAdd10<<<ceil(N/256), 256>>>(N, devA, devB);

// C++ 标准库中的 map（幻灯片原例）：
#include <algorithm>
int f(int x) { return x + 10; }
int a[] = {3, 8, 4, 6, 3, 9, 2, 8};
int b[8];
std::transform(a, a + 8, b, f);      // 输入起止迭代器 + 输出起始迭代器 + 一元函数

【代码做了什么？】

• CUDA 版：每个线程负责一个下标 i，读 a[i]、写 b[i]，f 是无副作用的纯函数。if (i < N) 是”线程数显式、数据集合大小不决定线程数”的典型护栏（对应第 7 讲：launch 的线程数不由数据规模决定）。
• C++ 版：std::transform(first1, last1, d_first, unary_op) 就是标准库的 map——输入两个迭代器界定序列、输出起始迭代器、一元操作符。
• 两者都在表达同一件事：b = map f a。

【并行机制解说】

• 为什么 map 可任意并行：f 无副作用，所以对每个元素的处理互不干扰、顺序任意。map 的并行化策略（幻灯片）：

  map f s =
      把序列 s 分成 P 个子序列
      对每个子序列 s_i（并行地）：
          out_i = map f s_i
      out = 拼接所有 out_i
  

• CUDA 的实现：block/grid 就是”划分子序列”，每个线程独立应用 f——零通信、零同步，GPU 上天然高效。
• 对应概念：map、数据并行模型、序列操作。

示例 2：树形 Reduce——CUDA 树归约（Work O(N)，Span O(log N)）

// tree_reduce.cu —— 编译：nvcc -o tree_reduce tree_reduce.cu
#include <cuda_runtime.h>

#define THREADS_PER_BLK 256

// 每个 block 把 256 个元素的子树归约为一个部分和
__global__ void treeReduce(const float* input, float* partial, int N)
{
    __shared__ float sdata[THREADS_PER_BLK];
    int tid = threadIdx.x;
    int i   = blockIdx.x * blockDim.x + tid;

    sdata[tid] = (i < N) ? input[i] : 0.0f;   // 越界补 0（0 是 + 的单位元）
    __syncthreads();

    // 树形归约：跨度每次减半，8 轮后 sdata[0] = 本 block 的和
    for (int s = THREADS_PER_BLK / 2; s > 0; s >>= 1) {
        if (tid < s)
            sdata[tid] += sdata[tid + s];     // 两两相加，跨度减半
        __syncthreads();
    }
    if (tid == 0)
        partial[blockIdx.x] = sdata[0];
}
// host 端再对 partial[0..numBlocks-1] 做第二次归约（可复用同一 kernel 或拷回 CPU）

【代码做了什么？】

• 每个线程装载一个元素到 shared memory，然后执行 8 轮”跨度减半”的树形相加：第 1 轮 128 个线程把 256 个元素两两相加 → 128 个部分和；第 2 轮 64 个线程 → 64 个部分和……第 8 轮 1 个线程 → 1 个和。整个过程像一个倒置的二叉树。
• 每轮之后 __syncthreads() 保证写方完成、读方再读。

【并行机制解说】

• 树形结构与 Work/Span：

  level 0:  256 个元素
  level 1:  128 个部分和（128 线程并行）
  level 2:   64 个部分和（64 线程并行）
  ...
  level 8:    1 个总和（1 线程）
  Work = O(N)（每元素参与一次加法），Span = O(log N)（8 轮串行依赖链）
  

• 这是 fold/reduce 原语的并行实现：每个 block 是一个”子归约”，block 之间靠第二次 kernel 合并——对应”并行 fold = 分段 fold + combiner 合并”。
• 注意每轮闲置一半线程（tid >= s 的直接跳过）——朴素实现的固有开销，幻灯片强调”work-efficient scan 的常数因子”时与此同理。
• 对应概念：reduce/fold、树形归约、Work/Span、shared memory 协作。

示例 3：scan_warp——32 元素 SIMD 并行 scan（Hillis-Steele 风格，幻灯片原码）

// 在 32 个 CUDA 线程（一个 warp）上执行 exclusive scan。
// 调用约定：由 32 个线程共同调用，ptr 指向 shared memory 中的 32 个元素；
// 每个线程返回自己下标对应的 exclusive scan 结果
// （完成后 ptr[] 里保存的是 inclusive scan 结果）。
__device__ int scan_warp(int *ptr, const unsigned int idx)
{
    const unsigned int lane = idx % 32;   // 线程在 warp 内的编号（0..31）
    __syncwarp();                         // warp 内同步（比 __syncthreads 便宜）

    for (int i = 0; i < 5; i++) {         // 5 步，因为 2^5 = 32
        int shift = 1 << i;               // 跨距 1, 2, 4, 8, 16
        if (lane >= shift) {
            int tmp1 = ptr[idx - shift];
            int tmp2 = ptr[idx];
            __syncwarp();                 // 读旧值之前先同步，防止读到被覆盖的新值
            ptr[idx] = tmp1 + tmp2;
            __syncwarp();                 // 写完后同步，防止下一轮读到半更新状态
        }
    }
    return (lane > 0) ? ptr[idx - 1] : 0; // exclusive：返回前一个位置的 inclusive 值
}

【代码做了什么？】

• 目标：32 个元素的前缀和（exclusive）。每步跨距翻倍：第 1 步每个线程把 ptr[idx] += ptr[idx-1]（跨距 1），第 2 步 ptr[idx] += ptr[idx-2]（跨距 2）……第 5 步跨距 16。5 步之后，ptr[idx] 里就是从 0 到 idx 的 inclusive scan。
• 每个线程最后返回 ptr[idx-1]（lane>0）或 0（lane=0）——即 exclusive scan 结果。
• 注意 __syncwarp() 的使用：在”读旧值”和”写新值”之间各放一次，保证同一 warp 内 32 个线程步调一致（warp 内同步比块内 __syncthreads() 便宜，因为只有一个 warp）。

【并行机制解说】

• 这是朴素并行 scan（每轮所有线程可并行，跨距翻倍）的教科书实现。它的复杂度：Work = N log N = 32×5 = 160 次加法——比串行 scan 的 31 次加法多得多！
• 为什么幻灯片说”work-efficient 的 scan 在这里反而不划算”：work-efficient（Blelloch）算法需要 up-sweep + down-sweep 两趟，指令数比这个实现多 2 倍以上，而且在 SIMD 上利用率低（up-sweep 时每轮线程数减半，down-sweep 时同样有半空通道）。在一个 warp 内（32 个 lane 全忙）跑 N log N 的 Hillis-Steele 反而 SIMD 利用率最高——”并行度只要够填满机器就行”的活例子。
• 对应概念：scan、inclusive/exclusive、Work vs Span、SIMD 利用率。

示例 4：scan_block——多 warp 协作的块级 scan（幻灯片原码，Assignment 3 同款）

// 在 CUDA thread block 内执行 scan。假设 ptr 指向 shared memory，
// 数组长度 == block 内线程数（Assignment 3 中提供了类似代码）。
__device__ void scan_block(int* ptr, const unsigned int idx)
{
    const unsigned int lane    = idx % 32;   // 线程在 warp 内的编号
    const unsigned int warp_id = idx >> 5;   // 线程所在 warp 在 block 内的编号

    // Step 1: 每个 warp 先对各自的 32 个元素做 scan_warp（部分 scan）
    int val = scan_warp(ptr, idx);
    // （所有线程都参与：同 warp 线程通过 shared 缓冲 ptr 通信）

    // Step 2: 每个 warp 的 31 号线程把本 warp 的扫描结果（最后一个元素）
    //         拷进 block 级 shared 的紧凑区域 ptr[0..numWarps-1]
    if (lane == 31) ptr[warp_id] = ptr[idx];
    __syncthreads();

    // Step 3: 只有 warp 0 对这 numWarps 个"基数"做一次 scan_warp，
    //         得到每个 warp 的偏移基数（inclusive）
    if (warp_id == 0) scan_warp(ptr, idx);
    __syncthreads();

    // Step 4: 所有线程把本 warp 的基数加到自己的部分 scan 结果上
    if (warp_id > 0)
        val = val + ptr[warp_id - 1];        // ptr[warp_id-1] 是前序 warp 的累计基数
    __syncthreads();

    ptr[idx] = val;                          // 写回最终 exclusive scan 结果
}

【代码做了什么？】

• Step 1：每个 warp 内先做 32 元素的 scan_warp（示例 3）——得到”warp 内局部 scan”。
• Step 2：每个 warp 的 31 号线程把该 warp 的局部总和（= 局部 inclusive scan 的最后一个值）写入 ptr[warp_id]（紧凑区域，只占 numWarps 个槽）。
• Step 3：warp 0 对 numWarps 个基数再做一次 scan_warp，得到”每个 warp 前面所有 warp 的总和”（warp 间的累计基数，inclusive 存放在 ptr[]）。
• Step 4：每个非 0 号 warp 的线程把自己的局部 scan 结果加上”前面 warp 的累计基数”，得到全局 exclusive scan。最后 __syncthreads() 保证写完才让下一个使用 ptr 的阶段开始。

【并行机制解说】

• 这是多级（heterogeneous）scan 策略的典范：warp 内用一种算法（SIMD 友好的 Hillis-Steele），warp 之间用另一种策略（先收缩成 numWarps 个基数、再扫描基数、再广播）。幻灯片称之为”算法在不同层级采用不同策略”——这是 scan 实现的关键洞见。
• 同步点：两次 __syncthreads()（Step 2→3、3→4）是块级屏障；scan_warp 内部的 __syncwarp() 是 warp 级屏障。块内通信完全走 shared memory。
• 对应概念：scan、分段扫描思想、多级实现、块内协作。

从块级 scan 到百万级 scan：三 kernel 流水（幻灯片第 24 页）

超过一个 block 能容纳的元素（如 100 万元素、每 block 1024 元素）时，”扫描基数”这一步本身也超过一个 block 的能力，需要把”计算 → 汇总 → 修正”拆成三个 kernel launch（每次 launch 之间有一次隐式全局屏障）：

Kernel Launch 1（逐块局部 scan）：
  Block 0  Scan ──► 局部 scan 结果 + 基数 base[0]
  Block 1  Scan ──► 局部 scan 结果 + 基数 base[1]
  ...      ...
  Block N-1 Scan ──► 局部 scan 结果 + 基数 base[N-1]
                （每个 block 只依赖自己的数据——并行）

Kernel Launch 2（扫描基数，规模小，一个 block 足够）：
  base[0..N-1] 做一次 scan ──► 每个 block 的累计偏移

Kernel Launch 3（逐块修正）：
  Block 0  Add base[0]  ──► 全局 scan 结果
  Block 1  Add base[1]  ──► 全局 scan 结果
  ...      ...
  Block N-1 Add base[N-1] ──► 全局 scan 结果
                （每个 block 又只依赖自己的数据——并行）

要点：每次 kernel launch 返回时隐式同步了所有线程，所以”基数必须先于修正”的依赖由 launch 边界天然满足——这就是 CUDA 中”用粗粒度同步（kernel 边界）替代细粒度同步（跨 block 通信）”的标准手法；其代价是数据被多遍历一遍（带宽）。

示例 5：用 gather + map + segmented scan 实现稀疏矩阵乘法（幻灯片实例）

目标：y = M·x，其中 M 是稀疏矩阵（大部分元素为 0），x 是稠密向量。
CSR（compressed sparse row）表示：
  values     = [3, 1, 2, 4, 2, 6, 8]        // 所有非零元素（按行展开）
  cols       = [0, 2, 1, 2, 1, 2, 3]        // 每个非零元素的列号
  row_starts = [0, 2, 3, 4]                 // 每行在 values 中的起始下标

5 步数据并行算法（每步都是一种原语）：
Step 1 (gather):  gathered[i] = x[cols[i]]
      gathered = [x0, x2, x1, x2, x1, x2, x3]
Step 2 (map):    products[i] = values[i] * gathered[i]
      products = [3x0, x2, 2x1, 4x2, 2x1, 6x2, 8x3]
Step 3 (构造 flag): 由 row_starts 生成 start-flag：每行起点为 1
      flags = [1, 0, 1, 1, 1, 0, 0]
Step 4 (segmented scan): 对 (products, flags) 做 inclusive segmented scan（+）
      [3x0, 3x0+x2, 2x1, 4x2, 2x1, 2x1+6x2, 2x1+6x2+8x3]
Step 5 (取每段末元素): 取每个 flag=1 段落的最后一个元素
      y = [3x0+x2, 2x1, 4x2, 2x1+6x2+8x3]

【代码做了什么？】

• 把”每行一个点积”（各行的非零数不同 → 不规则并行）改写为四条整齐的序列操作：先 gather 出每列对应的 x 值，再 map 乘上 values，再按行分段做 segmented scan 求和，最后取每段末元素作为该行结果。
• 复杂度：所有步骤都是 O(非零元素数) 的规则数据并行操作——没有锁、没有每行独立的线程（那会造成行间负载不均，破坏 SIMD）。

【并行机制解说】

• 为什么这是”数据并行思维”的胜利：直接按行并行，不同行非零数不同 → SIMD 利用率灾难（warp 内各 lane 工作量不同）；改写为”按非零元素并行 + segmented scan”后，每个 lane 的工作量完全相同——把不规则并行变成了规则并行（对应概念：segmented scan、gather、map）。
• 代价：多次遍历数据（gather → map → scan → 取末），带宽开销上升——这正是第 50 页总结里”带宽 hungry”的注脚。
• 对应概念：segmented scan、gather、map、数据并行思维。

示例 6：数据并行直方图——map + sort + scan 的组合拳（幻灯片第 47-49 页）

// histogram.cu（示意）—— 只用 map、sort、scan 类原语构造大规模并行直方图
// 目标：统计 input[0..N-1] 落入 NUM_BINS 个 bin 的数量（f 把值映射为 bin id）

// 阶段 1（map）：每个线程算自己元素的 bin id
__global__ void compute_bin(float* input, int* bin_ids)
{
    int thread_index = blockIdx.x * blockDim.x + threadIdx.x;
    bin_ids[thread_index] = f(input[thread_index]);
}

// 阶段 3（map + 差分检测）：排序后，扫描相邻元素，定位每个 bin 的起点
__global__ void find_starts(int* bin_ids, int* bin_starts)
{
    int thread_index = blockIdx.x * blockDim.x + threadIdx.x;
    if (thread_index == 0 || bin_ids[thread_index] != bin_ids[thread_index - 1])
        bin_starts[bin_ids[thread_index]] = thread_index;   // 该 bin 首次出现的下标
}

// 阶段 4（每 bin 一线程）：由 bin 起点算出 bin 大小（要跳过空 bin）
__global__ void bin_sizes(int* bin_starts, int* histogram_bins, int num_items, int num_bins)
{
    int thread_index = blockIdx.x * blockDim.x + threadIdx.x;   // 每线程一个 bin
    if (bin_starts[thread_index] == -1) {
        histogram_bins[thread_index] = 0;                       // 空 bin
    } else {
        // 找下一个非空 bin 的起点，两者之差即本 bin 大小
        int next_idx = thread_index + 1;
        while (next_idx < num_bins && bin_starts[next_idx] == -1)
            next_idx++;
        histogram_bins[thread_index] = (next_idx < num_bins)
            ? bin_starts[next_idx] - bin_starts[thread_index]
            : num_items - bin_starts[thread_index];
    }
}

// host 端流程：
//   bin_ids[N]  ← 由 compute_bin 填充；bin_starts[NUM_BINS] 初始化为 -1
//   sort(N, bin_ids, sorted_bin_ids);          // 阶段 2：按 bin id 排序（元素聚拢）
//   launch<<<N>>>     find_starts(sorted_bin_ids, bin_starts);   // 阶段 3
//   launch<<<NUM_BINS>>> bin_sizes(bin_starts, histogram_bins, N, NUM_BINS); // 阶段 4

【代码做了什么？】

• 阶段 1（map）：compute_bin 每线程算一个元素的 bin id——纯并行。
• 阶段 2（sort）：按 bin id 排序，让同一 bin 的元素聚拢成连续段（排序时元素数组本身也要随之置换，但直方图只关心计数）。
• 阶段 3（差分检测）：排序后”相邻两个元素的 bin id 不同”就说明这里是某 bin 的起点——find_starts 每线程检查自己与前一个元素，把起点下标写进 bin_starts[bin_id]。thread_index == 0 与最后一个元素是两个特例。
• 阶段 4（每 bin 一线程）：bin_sizes 由”本 bin 起点 − 下一个非空 bin 起点”得大小；用 while 跳过中间的空 bin（bin_starts == -1）。

【并行机制解说】

• 为什么不用”每个元素原子地给 bin 加 1”：原子方案不是错误（第 7 讲论证过原子互斥合法），但在 GPU 上争用是灾难——数据分布集中时几千个线程抢同一 bin 的原子单元，串行化到吞吐趋零。数据并行方案零细粒度同步：全程只有 map/sort/scan 式的并行遍历与只读比较，全部线程从不互斥写同一位置。
• 三种原语各司其职：map 负责”逐元素计算”，sort 负责”把相同 key 聚拢”（这是 groupByKey 的核心机制），scan/差分负责”把连续段边界变成可并行的计数”。
• 代价：sort 本身是超线性工作 + 多次遍历数组的带宽开销（”at the cost of a sort and extra passes over the data (extra BW)!”）。
• 对应概念：map、sort、groupByKey（排序聚拢）、数据并行直方图、带宽代价。

三、关键要点

1. 用”对序列的操作”思考算法：map、filter、reduce、scan、segmented scan、sort、groupBy、gather/scatter 的高性能并行实现已经存在（Thrust、Pandas、JAX、Spark/Hadoop 全是它们撑起来的），把程序改写成这些原语的组合，就自动获得了并行能力——前提是别被带宽卡住。
2. 理解依赖是关键：x = a + b; y = b*7; z = (x-y)*(x+y)——没有依赖的操作可以并行（a+b 与 b*7 并行），有依赖的必须等待（z 依赖 x、y）。scan 的难点全在”每个输出都依赖前面所有输入”这条长依赖链上。
3. 并行 scan 是 Work/Span 权衡的教科书：朴素算法 Work O(N log N)、Span O(log N)；工作高效算法（up-sweep+down-sweep）Work O(N)、Span O(log N)。但 GPU 上常选”低效”算法：warp 内 32 个 lane 全忙的 N log N 版，比”高效但指令多 2 倍以上、通道利用率低”的 Blelloch 版更快。
4. 多级实现（hierarchy）：scan 在不同层级用不同策略——warp 内 SIMD scan、块内 warp 间协作（收缩基数 + 再扫描 + 广播）、跨块用多 kernel launch。目标永远是：减少 work、减少通信/同步、匹配内存层级。
5. 数据并行思维的三大功效与代价：把不规则并行→规则并行；把细粒度同步→粗粒度同步（甚至零同步）；代价是多趟遍历 → 带宽压力大。粒子网格的五种解法完美演示了这条光谱：全局锁（细粒度、高争用）→ 每格锁 → 按格并行（并行度不足）→ 部分结果+合并 → sort+scan 的纯数据并行方案（无锁、大并行度，代价是一次 sort 和额外带宽）。

四、常见陷阱与注意事项

1. 以为 scan 的输出可以”独立并行计算”：map 可以任意乱序，但 scan 的每个输出都依赖前缀，不能直接并行——必须用跨距翻倍或 up/down-sweep 等专门算法；把 inclusive/exclusive 搞混（exclusive 的第一个元素是单位元 I，不是 in[0]）也是高频错误。
2. 忽略结合律要求：并行 scan/reduce 要求 ⊕ 是结合律操作，且（并行 fold 时）初始值必须是单位元。浮点加法”近似结合”但不严格结合——并行归约的顺序不同，结果可能与串行版本有微小数值差异（这是语义上可接受的，但要知道）。
3. __syncwarp() / __syncthreads() 放错位置：scan_warp 中”读旧值→同步→写新值→同步”的顺序一旦错乱（比如漏掉写后同步），warp 内线程会读到半更新的数据。块级 scan 中 Step 2→3、3→4 的 __syncthreads() 一个都不能少。
4. 带宽盲区：数据并行方案（如稀疏矩阵乘法、粒子网格的 sort 方案）需要多趟全量遍历，每趟都是 O(N) 带宽。若算术强度低，最终瓶颈是带宽而不是并行度——”efficient implementations only leverage as much parallelism as required”。
5. 把 block 当线程池 / 线程数当元素数：GPU 需要海量并行度（V100 可并发 163,840 个 CUDA 线程）；只有 16 个格子的”按格并行”方案并行度严重不足，而”每粒子一原子操作”又争用爆炸——要在并行度、同步开销、带宽之间取平衡（粒子网格 5 种解法就是这份权衡的完整案例）。
6. 对不规则数据强行 SIMD：行长度不同的稀疏矩阵如果”一行一线程”，warp 内 lane 工作量差异巨大（SIMD 利用率灾难）；正确姿势是先压平成规则序列 + segmented scan（示例 5）。

五、思考题（带答案）

Q1：为什么幻灯片说朴素并行 scan（Hillis-Steele，Work = N log N）”比串行算法还低效”，却在 CUDA 的 scan_warp 中仍然使用它？在什么条件下你会选择 work-efficient（Blelloch）算法？

A1：从 Work（总操作数）看，N log N > N，朴素算法确实更”费算力”；但 GPU 的性能模型不是只看 Work，而是看执行时间 ≈ max(Work/并行度, Span) × 常数。在单个 warp（32 lane）里，Hillis-Steele 每轮 32 个 lane 全忙，5 轮搞定（Span = 5），时间 ≈ 160/32×每指令时间；而 Blelloch 需要 up-sweep（每轮线程减半）+ down-sweep（同样有半空轮），虽然总加法数约 2N，但指令数多 2 倍以上且 SIMD 通道利用率低，实际更慢。选择原则：当可用并行度（lane/核数）≥ 数据规模时（如 warp 内 32 元素），用高并行度低 span 的朴素算法；当数据规模远大于可用并行度、且每元素工作较重时，用 work-efficient 算法把总工作量压下来。这也是”efficient implementations only leverage as much parallelism as required”的注脚。

Q2：在”粒子网格”问题（把 100 万个粒子按 2D 位置放进 16 个格子）的五种解法中，为什么”数据并行（sort+scan）”方案最终胜出？它牺牲了什么？

A2：前四种方案要么争用严重（全局锁：数千线程抢一把锁）、要么并行度不足（16 个格子只有 16 个任务，远小于 GPU 需要的数万）、要么浪费算力（每格扫描全部粒子 = 16 倍粒子-格子判定）。sort+scan 方案：Step 1 map 算每个粒子的格子号（并行于粒子）；Step 2 按格子号 sort（粒子索引同步置换）；Step 3 并行扫描排序后的数组，用”相邻元素格子号不同”定位每格 start/end。它的优点：并行度始终是 O(粒子数)，且完全不需要锁（排序天然把同类聚在一起，start/end 由位置关系确定）。牺牲：一次全局 sort 和多次额外遍历带来的带宽开销（”at the cost of a sort and extra passes over the data (extra BW)!”）。这正体现了本讲总结：用规则并行替代不规则并行、用粗粒度（甚至零）同步替代细粒度同步，代价是带宽。

Q3：如何只用 map、sort、segmented scan 构造一个大规模并行直方图？为什么不能用”每个元素原子地给对应 bin 加 1”替代？

A3：数据并行直方图三阶段：① map：compute_bin——每个线程算 bin_ids[i] = f(input[i])；② sort：按 bin_id 排序得到 sorted_bin_ids（同类聚拢）；③ find_starts + bin_sizes：并行扫描排序数组，bin_starts[bin_id] 记录每个 bin 首次出现的下标；再按”下一个非空 bin 的 start 减当前 start”算出每个 bin 的大小（bin_sizes kernel 用 while 跳过空 bin）。至于”原子加 1”方案：它并非不正确（第 7 讲已论证原子互斥不破坏块间调度），但在 GPU 上争用是灾难——如果数据分布集中，成千上万个线程同时对少数几个 bin 做 atomicAdd，原子单元会串行化到接近零吞吐；而 sort+scan 方案完全没有细粒度同步，全部线程只读。代价依然是 sort 和额外遍历的带宽。