Lecture 26: MPI, OpenMP, Cilk Implementation

目录 · ← l24 · l26 →

Lecture 26: MPI, OpenMP, Cilk Implementation

1. 章节标题与概述

Lecture 26: MPI, OpenMP, Cilk Implementation

• 本讲核心问题：前两周我们从”抽象层”讨论了三大并行编程模型（shared address space（共享地址空间）、message passing（消息传递）、data parallel（数据并行）），本讲把镜头拉到”实现层”：MPI、OpenMP、Cilk 这三种具体编程系统究竟由谁、用什么机制把程序员书写的并行抽象映射到硬件上？具体要回答三个问题：(1) 一个 cilk_spawn / #pragma omp parallel for / MPI_Send 调用在运行时到底做了什么？(2) 三种系统在 decomposition（分解）、assignment（分配）、orchestration（编排）、mapping（映射）四个阶段分别把责任交给了程序员、编译器、运行时还是硬件？(3) 这些实现机制各自的性能上界由什么决定——是计算吞吐、内存带宽、同步延迟，还是关键路径（span）？
• 涉及的主要硬件/软件机制：
  • 硬件侧：集群互连网络（InfiniBand/NIC、fat-tree 拓扑、bisection bandwidth（对分带宽）、oversubscription（超额订阅））、单节点多核 + 共享 LLC 与 NUMA、cache line（64 字节）与 cache coherence（缓存一致性）协议、原子 read-modify-write（RMW）指令（lock xadd / cmpxchg）、消息延迟 α 与带宽 β。
  • 软件侧：MPI 库实现中的 eager（急切）与 rendezvous（会合）协议、unexpected-message queue（非预期消息队列）、非阻塞 Isend/Irecv 与 Waitall、集合通信的树形算法；OpenMP 的编译器 outlining（轮廓化）变换 + 线程池运行时（libgomp/libomp）+ worksharing 循环切块 + barrier（屏障）实现（集中式 sense-reversing、combining tree、dissemination）；Cilk 的 cilk_spawn/cilk_sync → continuation stealing（续体窃取） + 每 worker 一个 lock-free dequeue + randomized work stealing（随机工作窃取）+ reducers（hyperobject）。
• 在并行计算知识体系中的角色：本讲是”抽象 → 实现”的闭合环节。前面课程已经建立了三大模型（models）、依赖与同步（synchronization）、缓存一致性与内存序（coherence/consistency）、互连网络（interconnects）、work-span 与 Amdahl 定律等工具；本讲说明这些工具如何共同落成真实代码：message passing 落地为 MPI 的通信协议与集合通信算法，shared address space 落地为 OpenMP 的线程池 + 屏障 + 调度器，fork-join/divide-and-conquer 落地为 Cilk 的 work-stealing 运行时。它也是后续”性能调优 / 项目实现”的直接方法论来源：先用最简单的实现跑通，再依据瓶颈（带宽、延迟、span、同步次数）选择机制。
• 配套材料：
  • Fall 2026 课程主页与日程表：https://www.cs.cmu.edu/~418/、https://www.cs.cmu.edu/~418/schedule.html —— 已公开（无需登录）。日程表列出 Lecture 26 主题 “MPI, OpenMP, Cilk implementation”（Nov 3），并注明 Project proposal due。
  • 本讲的 Fall 2026 讲义 PDF：在 https://www.cs.cmu.edu/~418/lectures/ 下尚未发布（Fall 2026 已公开 PDF 只覆盖 Why Parallelism、ILP、Modern Multi-Core、Programming Models、CUDA、Parallel Programming Basics、Performance Optimization I/II、Interconnects、Coherence、Synchronization、Lock-Free、Heterogeneity、Virtual Memory、Consistency、Parallel Deep Learning 等）。日程表中 Lecture 26 一行引用的是历史学期（Fall 2022）的 “Part A slides”（26a-msgpassing.pdf）与 “Part B slides”（26b-omp.pdf），它们位于 /afs/cs/academic/class/15418-f22/public/lectures/ 之下，访问需要 CMU 登录 —— 属未公开。
  • 抽取尝试（本地）：extracted/f22_26-msgpassing.txt、extracted/f22_26a-msgpassing.txt、extracted/s23_26b-omp.txt 三个文件都只含登录页 HTML（”Your Browser does not support javascript. Please use this link.”，88 字节），无讲义正文，不能作为事实来源，已跳过。
  • 录像：Fall 2026 日程表中 video 链接均被 HTML 注释隐藏（”video from a previous offering; uncomment when posted for Fall 2026”），属未发布；Ed 讨论区、Autolab、Canvas/Gradescope 需登录，属需登录未公开。
  • 本笔记的事实基础：Fall 2026 已公开讲义 extracted/06_progbasics.txt（Lecture 7: Parallel Programming Basics，页眉标注 Fall 2026）与 extracted/07_progperf1.txt（页眉标注 Fall 2025，属讲义沿用的历史学期字样，内容为 Performance Optimization Part 1: Work Distribution and Scheduling）。前者提供：三大模型对比、分解/分配/编排/映射四阶段、Amdahl 定律与其算例、锁与屏障、message passing 网格求解器（ghost cell（影子单元）、同步 send/recv 死锁、非阻塞 send/recv）；后者提供：静态/半静态/动态分配、任务粒度权衡、共享工作队列、每 worker 分布式队列与 work stealing、Cilk Plus 的 cilk_spawn/cilk_sync 语义、continuation stealing vs child stealing、每 worker dequeue（本地压栈尾部、窃取头部）、parallel slack（并行松弛 ≈ 8）。本讲其余 MPI/OpenMP/Cilk 实现细节依据公开成熟知识（MPI 标准与实现文档、OpenMP 规范与 libgomp/libomp 实现、Cilk-5 / OpenCilk 相关论文与技术资料）撰写，并已与上述讲义中的术语和结论对齐。
  • Fall 2026 授课教师：Brian Railing 与 Dimitrios Skarlatos；课程由 Kayvon Fatahalian 创建。

2. 核心概念与硬件/软件架构图解

2.1 谁负责哪一步：从抽象到实现的责任分层

• 定义与目的：并行程序从”问题”到”执行”要经过四个阶段——decomposition（分解，把问题切成任务）→ assignment（分配，把任务交给 worker）→ orchestration（编排，组织通信/同步/数据布局/调度）→ mapping（映射，把 worker 落到硬件执行单元）。三大编程系统的差别，本质上就是这四个阶段的责任归属不同：MPI 把四件事几乎全交给程序员；OpenMP 把 assignment 交给运行时（调度子句可调）、orchestration 部分交给编译器 + 运行时；Cilk 把 assignment 和 mapping 完全交给运行时（work-stealing 调度器）。
• 直观解释（”它是什么？”）：把它想象成装修一栋楼。decomposition = 决定有哪些活（水电、刷墙、铺地板）；assignment = 决定哪个工头带哪批活；orchestration = 规定工头之间怎么交接（要不要开会、材料放哪）；mapping = 实际派几个工人、派到哪几个房间。MPI 模式下你是总包 + 每个工头 + 每个工人，所有交接都必须你手写电话（显式消息）；OpenMP 模式下你只写”这面墙刷两遍”，公司（运行时）自动分工头、自动开会（隐式屏障）；Cilk 模式下你只说”这活能并行”，调度中心（work-stealing runtime）看到谁闲就把活塞给谁。
• 架构/机制图解（软件执行模型）：

                     decomposition   assignment      orchestration        mapping
                     (谁切任务)       (谁派任务)       (谁管通信/同步)       (谁落硬件)
   ┌──────────────┬───────────────┬───────────────┬───────────────────┬──────────────────┐
   │  裸 pthread  │   程序员       │   程序员       │   程序员           │  OS 调度器 → 核   │
   │  共享内存    │               │               │  (锁/屏障/原子)     │                  │
   ├──────────────┼───────────────┼───────────────┼───────────────────┼──────────────────┤
   │  MPI         │   程序员       │   程序员       │   程序员           │  MPI 库 → NIC/网络│
   │  消息传递    │  (域分解)      │  (rank 编号)   │  (send/recv 显式)  │  (进程↔节点映射)  │
   ├──────────────┼───────────────┼───────────────┼───────────────────┼──────────────────┤
   │  OpenMP      │   程序员       │ 运行时(调度子句)│ 编译器插桩+运行时  │  线程池 → core/HT │
   │  共享地址空间│  (#pragma 标注)│ static/dynamic │  (隐式屏障/atomic) │  (OMP_PLACES)    │
   ├──────────────┼───────────────┼───────────────┼───────────────────┼──────────────────┤
   │  Cilk        │   程序员       │ 运行时         │ 运行时            │  运行时 worker 池  │
   │  fork-join   │  (spawn/sync)  │ (work stealing)│ (隐式 sync/串行化) │  → 执行上下文     │
   ├──────────────┼───────────────┼───────────────┼───────────────────┼──────────────────┤
   │  ISPC / CUDA │   程序员       │ 系统/硬件      │ 系统 + 内置原语    │  编译器 → SIMD lane│
   │  (数据并行)  │  (foreach/grid)│  (静态/块调度) │  (reduceAdd 等)    │  硬件 → SM/block  │
   └──────────────┴───────────────┴───────────────┴───────────────────┴──────────────────┘
       越往左 = 程序员责任越重、控制力越强、第一条正确程序越难写
       越往右 = 系统责任越重、开发越快，但性能上限受运行时机制约束

• 关键操作与性能特征：责任越是交给系统，运行时开销就越”隐形但真实”——OpenMP 每次进入 parallel for 是一次 fork-join（线程唤醒 + 屏障，量级 1–10 µs）；Cilk 每次真的发生窃取（steal）要访问别的 worker 的 dequeue（量级 0.1–5 µs）；MPI 每次通信要付 α 启动延迟（节点内 ~0.5–2 µs，跨节点 ~1–5 µs）。判断标准不是”谁更优雅”，而是”每次迭代的计算量能否摊薄这些开销的固定成本”。

表 1：三种编程系统在四个阶段与失败模式上的对比

2.2 MPI 实现：独立地址空间上的显式消息

• 定义与目的：MPI（Message Passing Interface）是一个库标准（不是语言，也不是编译器特性），它把”通信”结构化为一对一对的消息（message）：发送方调用 MPI_Send，接收方调用 MPI_Recv，库负责把数据从发送方地址空间搬到接收方地址空间。它解决的核心问题是：在节点之间没有共享内存（甚至没有统一地址空间）时，如何表达通信与同步——正如讲义所说，message passing 模型下”每个线程拥有自己的地址空间，没有共享地址空间抽象，线程通过收发消息来通信与同步”。
• 直观解释（”它是什么？”）：把集群想象成两家没有共享仓库的公司，只能靠快递传货。MPI_Send = 把货打包交给快递（数据从你的缓冲区复制到网络缓冲区）；MPI_Recv = 从前台取货（数据从网络缓冲区复制到你的缓冲区）。关键点：你不是”共享内存 + 网络硬件”，而是”复制 + 复制 + 复制”——每次通信至少三次数据搬运（用户缓冲 → 库缓冲 → 网络 → 库缓冲 → 用户缓冲），所以 MPI 的优化重点永远是少发、发大块（bulk transfer），正如讲义强调的”communicate entire rows at a time (not individual elements)”。
• 架构/机制图解（硬件结构 + 消息路径）：

   Rank 0 (CPU)                                    Rank 1 (CPU)
 ┌────────────────────────┐                     ┌────────────────────────┐
 │ 用户缓冲 sendbuf       │                     │ 用户缓冲 recvbuf       │
 │  (C 数组, 页锁定可选)  │                     │                        │
 └──────────┬─────────────┘                     └───────────▲────────────┘
            │ ① memcpy (库内部, ~10 GB/s)                   │ ⑤ memcpy
 ┌──────────▼─────────────┐                     ┌───────────┴────────────┐
 │ MPI 库缓冲 /           │                     │ MPI 库: 匹配队列       │
 │ unexpected queue       │                     │  posted recv /         │
 │ 快路径: 命中 posted    │                     │  unexpected msg        │
 └──────────┬─────────────┘                     └───────────▲────────────┘
            │ ② doorbell (MMIO 写)                          │ ④ 中断/轮询
 ┌──────────▼─────────────┐   fat-tree / 交换芯片  ┌────────▼───────────┐
 │ NIC (HCA)              │ ③ 包序列化 + 链路传输  │ NIC (HCA)          │
 │  ─ 发送队列 (WQE)      │◄══════════════════════►│                    │
 │  ─ 完成队列 (CQE)      │   InfiniBand HDR/NDR   │                    │
 └────────────────────────┘                        └────────────────────┘
      时间轴(延迟构成, 小消息 ~1–5 µs):
      ├─ 软件路径(库+驱动) ~0.5–2 µs ─┼─ NIC 序列化&仲裁 ~0.1 µs ─┼─
      ├─ 链路传播(每跳 ~5 ns/m 级) ───┼─ 交换机排队(拥塞) ────────┼─
      ├─ 接收侧软件路径 ~0.5–2 µs ────┼─ 匹配 + memcpy ──────────┤
      大消息时: 时间 ≈ α + n/BW, 传播与排队被带宽项淹没

• 关键操作与性能特征：
  • 延迟（latency）：小消息成本 ≈ α（启动延迟），与消息大小几乎无关。节点内共享内存 MPI 约 0.3–1 µs，InfiniBand 跨节点约 1–5 µs。α 是”每条消息”的固定税，因此把 1000 个 8 字节的消息合成 1 个 8 KB 的消息通常能提速一个数量级以上。
  • 带宽（bandwidth）：大消息成本 ≈ n/BW。单端口 HDR 200 Gb/s ≈ 25 GB/s（理论），实际可达约 10–20 GB/s；带宽项由链路速率、NIC DMA 引擎、PCIe 通道与内存带宽共同决定。
  • 消息速率（message rate）：很多实现能到 5–20 M msg/s/NIC；这是”小消息”场景的真实上界，与带宽无关。
  • eager vs rendezvous：库在消息大小超过阈值（典型 8–64 KB）时从 eager 切换到 rendezvous，以避免占用不可控的接收缓冲。

表 2：MPI 发送/接收协议对比（实现层视角）

2.3 MPI 域分解与 ghost cell：表面-体积比

• 定义与目的：把 N×N 网格按行切成 size 块，每块带两行 ghost cell（影子单元）——它们不是本 rank 拥有的数据，而是从邻居地址空间复制来的副本。ghost cell 存在的唯一目的：让 stencil（模板）计算在本地地址空间里就能拿到所有邻居值，而不必每次访问都发消息。讲义明确说：”data replication is now required to correctly execute the program”，并且”information in ghost cells is ‘owned’ by other threads”。
• 直观解释（”它是什么？”）：像拼图时把邻块边缘那一条先拓印到自己的板子上。你不需要知道邻居怎么拼，只需要边缘那一条。
• 架构/机制图解（软件执行模型：环状 halo 交换）：

        全局 N x N 网格按行分给 4 个 rank（blocked 分配）
   ┌───────────────────────────────────────────────┐
   │ 全局边界 (固定值 1.0)                          │  ← 只有边界 rank 需要
   ╞═══════════════════════════════════════════════╡
   │  rank 0: rows 0 .. r-1       [上行 ghost=边界] │
   │                              下行 ghost ◄──────┼── rank1 发来的第 1 行
   ╞═══════════════════════════════════════════════╡
   │  rank 1: rows r .. 2r-1     上行 ghost ◄───────┼── rank0 发来的最后 1 行
   │                              下行 ghost ◄──────┼── rank2 发来的第 1 行
   ╞═══════════════════════════════════════════════╡
   │  rank 2: rows 2r .. 3r-1    上行 ghost ◄───────┼── rank1 发来的最后 1 行
   │                              下行 ghost ◄──────┼── rank3 发来的第 1 行
   ╞═══════════════════════════════════════════════╡
   │  rank 3: rows 3r .. N-1     上行 ghost ◄───────┼── rank2 发来的最后 1 行
   │                              [下行 ghost=边界] │
   └───────────────────────────────────────────────┘

  一次迭代的数据流(每个 rank):
     ① Isend 上行 1 行 ──►  ② Isend 下行 1 行   (先发起, 不等待)
     ③ Irecv 上行 ghost ◄── ④ Irecv 下行 ghost
     ⑤ Waitall(4 个请求)  ◄── 通信期间 CPU 可做别的活(如同步旧数据)
     ⑥ 计算本地 rows x (N-2) 个内点
     ⑦ Allreduce 收敛判据(所有 rank 得到同一结论, 避免有人提前退出)

  表面-体积比:
     通信字节 / 计算字节 = (2 * N * 4B) / (rows * N * 20B) = 0.4 / rows
     rows 越大 → 通信占比越低 (这就是"weak scaling 友好"的根源)

• 关键操作与性能特征：通信量与”表面积”成正比，计算量与”体积”成正比，因此每个 rank 负责的行数 rows 越大，通信/计算比越低。这也是讲义中”blocked assignment requires less data to be communicated between processors”的定量版本——相比 interleaved（交错）分配，blocked 分配的边界长度从 O(N·P) 降到 O(N)。

2.4 MPI 集合通信：用树把 O(P) 变成 O(log P)

• 定义与目的：MPI_Reduce、MPI_Allreduce、MPI_Bcast、MPI_Scatter/Gather 等集合操作（collective）如果用 P 个 rank 的朴素”全都发给 rank 0”或”rank 0 发给所有人”实现，延迟随 P 线性增长，且 rank 0 的 NIC 成为热点。实现层的解法：树形算法（binomial tree / recursive doubling / halving），把步数从 O(P) 降到 O(log P)。
• 直观解释（”它是什么？”）：像公司逐级汇报：不是 64 个员工都跑去找 CEO，而是两两合并、16 人合 8 人、8 人合 4 人……6 轮就到 CEO。广播则相反：CEO 只讲给 2 个副总，副总各讲给 2 个经理，逐层扩散。
• 架构/机制图解（集合通信的树形实现）：

  MPI_Bcast(root=0) 的 binomial tree (P = 8, 3 轮完成, O(log P) 深度)

   轮 1          轮 2                轮 3
   ┌───┐        ┌───┐              ┌───┐
   │ 0 │───────►│ 1 │─────────────►│ 3 │
   └─┬─┘        └───┘              └───┘
     │
     │          ┌───┐              ┌───┐
     └─────────►│ 2 │─────────────►│ 5 │
                └─┬─┘              └───┘
                  │
                  │                ┌───┐
                  └───────────────►│ 6 │
                                   └───┘
   同时 rank1 在第 2 轮可以把数据给 rank 4/5 的分支...（每轮活跃 rank 数翻倍）
   总步数 = ceil(log2 P)；每步成本 ≈ α + n/BW
   T_bcast ≈ log2(P) * (α + n/BW)          ← 延迟优先(小消息)
   若 n 很大, 更优的是"分段流水线/递归倍增"变体, 使 T ≈ α*logP + n/BW

  MPI_Allreduce(P=8) = reduce(树形向上, log P 步) + bcast(树形向下, log P 步)
     ≈ 2 * log2(P) * α + 2 * (P-1)/P * n/BW
  朴素实现 = 7 个 rank 各发一次给 0, 再等 7 次广播: 14 步 + rank0 NIC 热点 ← 差

• 关键操作与性能特征：集合通信的实现在”延迟最优”（小消息：log P 步）与”带宽最优”（大消息：分段 + 环形/递归倍增，使每条链路流量为 (P−1)/P·n）之间切换，切换依据同样是消息大小。对使用者的直接结论：每轮迭代只做必要的集合操作；收敛判据（如 Jacobi 的 diff）必须用 Allreduce 而不是 Reduce + 广播两次，且能用 Allreduce 的场合不要用”大家都 Reduce 到 rank 0”。

2.5 OpenMP 实现：编译器 outlining + 线程池 + 屏障 + 调度器

• 定义与目的：OpenMP 是指令制导（directive-based）的共享地址空间编程接口。编译器把 #pragma omp parallel for 的循环体轮廓化（outlining）成一个独立函数，生成”启动并行区 + 按线程切分迭代空间 + 结束时屏障”的代码，运行时库（GCC 的 libgomp，Clang/ICC 的 libomp）负责线程池、屏障与调度。它解决的问题是：让共享内存并行化的样板代码由编译器生成，程序员只描述”哪段循环可并行”。
• 直观解释（”它是什么？”）：像公司里的”会议制”：老板（主线程）喊一声”开会”，已经在休息室待命的员工（线程池）进会议室；会上把活按规则一分（worksharing），干完必须都到齐（barrier）才能散会。关键实现细节是”员工不是每次招新人，而是常驻休息室”——线程创建开销（数十微秒）在第一次并行区付一次，之后每次 fork-join 只需唤醒 + 屏障（微秒级）。
• 架构/机制图解（软件执行模型 + 硬件映射）：

  OpenMP fork-join 执行模型与线程池 (OMP_NUM_THREADS=4, 8 个硬件线程上下文)

  主线程                                          worker 池(懒创建: 首次 parallel 时)
  ┌─────┐   fork(唤醒)   ┌──────┐┌──────┐┌──────┐┌──────┐
  │ T0  │───────────────►│ T0   ││ T1   ││ T2   ││ T3   │
  └──┬──┘                └──┬───┘└──┬───┘└──┬───┘└──┬───┘
     │  parallel for         │       │       │       │
     │  schedule(static,4)   ▼       ▼       ▼       ▼
     │                    ┌───────────────────────────────┐
     │                    │ 迭代空间 0..15 按块切分        │
     │                    │ T0: 0-3  T1: 4-7  T2: 8-11 T3:12-15 │
     │                    └───────────────┬───────────────┘
     │  implicit barrier                  │ 全部到齐
     ◄────────────────────────────────────┘
     │  join (worker 回到自旋/休眠)
     ▼
  继续串行代码

  映射(4 线程 -> 8 硬件上下文):
   ┌──────────────────────────────────────────────────────────┐
   │ Socket 0 (NUMA 0)              │ Socket 1 (NUMA 1)        │
   │ ┌────────┐┌────────┐┌────────┐ │ ┌────────┐┌────────┐ ... │
   │ │Core0+T0││Core1+T1││Core2+T2│ │ │Core4+T3││Core5   │     │
   │ │L1  L2  ││L1  L2  ││L1  L2  │ │ │L1  L2  ││L1  L2  │     │
   │ └───┬────┘└───┬────┘└───┬────┘ │ └───┬────┘└───┬────┘     │
   │     └─────────┴─────────┘      │     └─────────┘          │
   │          共享 L3 (LLC)         │      共享 L3 (LLC)       │
   └───────────────┬────────────────┴───────────┬─────────────┘
                   └────── 片间互连 / 内存 ───────┘
   OMP_PROC_BIND=close 让线程贴紧核心、不迁移; OMP_PLACES=cores 绑定到物理核
   → 绑定失效会带来跨 NUMA 远端内存访问, 延迟翻倍、带宽减半

• 关键操作与性能特征：
  • fork/join 成本：并行区进入 + 退出（含屏障）典型 1–10 µs，取决于线程是否在自旋（OMP_WAIT_POLICY=ACTIVE 延迟低但抢 CPU；PASSIVE 省电但唤醒延迟高）。
  • 屏障是”最保守的依赖表达”：讲义指出 barrier “divide computation into phases”，且 “all computations by all threads before the barrier complete before any computation in any thread after the barrier begins”。它的代价是所有线程的最长者决定相位长度（负载不均被放大），且每次都是一次全局同步。
  • 调度开销与并行度：schedule(dynamic, c) 的计数器是串行的，因此 chunk 越小，理论可达到的并行度上限越低（见第 4 节公式）。schedule(guided) 用递减的 chunk（约 剩余/线程数）在两者之间折中。
  • 数据环境：reduction(+:sum) 的实现是”每线程私有副本 + 结束时合并”，因此不要在并行区内加锁更新 sum（讲义中的求解器示例正是用 myDiff 局部累加、循环外只加一次锁，把锁次数从 O(N²) 降到 O(线程数)）。

表 3：OpenMP 调度策略对比

2.6 屏障的实现：从集中式到树形

• 定义与目的：屏障（barrier）是”相位同步”的硬件/软件原语。实现层的核心问题：P 个线程如何在没有专用硬件的情况下高效达成”全部到齐”。朴素实现用一个计数器 + 锁 → 所有线程争抢同一条 cache line → 串行化；工程实现用感知反转（sense-reversing）避免”屏障重用”的正确性问题，用树形/传播式（combining tree / dissemination）把 O(P) 的关键路径降到 O(log P)。
• 直观解释（”它是什么？”）：集中式屏障像只有一支签到笔的签到台，64 个人排队签；树形屏障像分组签到：4 人一组先在小表上签到，组代表再去上一层签，只有 log₂64 = 6 层。
• 架构/机制图解（状态机 + 数据流）：

  集中式 sense-reversing barrier —— 状态迁移图

        arrive(i)                 last thread sees
        count++                   count == P
   ┌──────────────┐   count<P   ┌──────────────────┐
   │  ARRIVING    │────────────►│  LAST_ARRIVAL    │
   │ (线程自旋在   │             │ (最后一个线程:    │
   │  local_sense)│             │  count=0;        │
   └──────┬───────┘             │  sense = !sense) │
          │ local_sense          └────────┬─────────┘
          │ == global sense               │ 翻转 global sense
          │ (被唤醒)                       ▼
   ┌──────▼───────┐   epoch++   ┌──────────────────┐
   │  RELEASED    │◄────────────│  BROADCAST       │
   │ (可进入下一相位)│            │ (所有自旋线程看到  │
   └──────────────┘             │  sense 变化)      │
                                └──────────────────┘

   伪代码(每个线程):
     local_sense = !local_sense;            // 每相位翻转 ⇒ 不会被上一相位误唤醒
     if (fetch_add(&count, 1) == P-1) {     // 我是最后一个
         count = 0;                          // 必须先清零再翻 sense
         global_sense = local_sense;         // 释放
     } else {
         while (global_sense != local_sense) ;   // 自旋(可加 pause/退避)
     }

  树形 barrier 的数据流 (P = 16, 4 叉树): ARRIVE 自底向上计数, RELEASE 自顶向下唤醒
      叶子(线程)        中间层(每 4 个叶子到齐→代表上行)      根
     ┌──┬──┬──┬──┐        ┌──────────┐                ┌──────────┐
     │L0│L1│L2│L3│──┐     │ 内部节点 │  上行计数       │   根     │
     └──┴──┴──┴──┘  ├────►│ count==4 │───────────────►│ count==4 │
     ┌──┬──┬──┬──┐  │     └────┬─────┘                └────┬─────┘
     │L4│L5│L6│L7│──┘          │ RELEASE: 统一翻 sense      │
     └──┴──┴──┴──┘◄────────────┘（自顶向下逐层唤醒）◄────────┘
   关键路径 O(log_k P) 层 × 2(上行+下行) 次 cache line 传输; 每次 ≈ 30–80 ns (同 socket)
   P=64: 集中式 ≈ 64 次串行争抢 × ~50 ns ≈ 3.2 µs
         4 叉树  ≈ 3 层 × 2 × ~50 ns ≈ 0.3 µs   ← 快一个数量级

• 关键操作与性能特征：屏障成本 ≈（关键路径上的 cache line 传输次数）×（一次跨核传输延迟 30–100 ns）+ 唤醒延迟。集中式屏障在 P ≤ 16 时通常够快（一次原子操作 + 广播）；P 增大到几十上百时树形/传播式明显更优。对使用者的结论：屏障频率比屏障实现更重要——如果每次迭代的计算只有几微秒，即使 0.3 µs 的树形屏障也会吃掉 10% 以上的时间。

2.7 Cilk 实现：continuation stealing + 每 worker dequeue

• 定义与目的：Cilk 用两个关键字表达 fork-join：cilk_spawn foo(args);（”调用 foo，但调用者可以继续异步执行”，即 fork）与 cilk_sync;（”等本函数内所有 spawn 完成”，即 join）；并且”每个含 cilk_spawn 的函数末尾都有一个隐式 cilk_sync“，因此”当一个 Cilk 函数返回时，与它相关的所有工作都已完成”。实现层的核心是：把逻辑上的 fork-join 变成”每 worker 一个队列 + 空闲者偷活”的调度器，使得串行化（删掉 spawn/sync）后的语义与串行程序完全一致（serial elision，串行消除）。
• 直观解释（”它是什么？”）：Cilk 运行时像一个共享办公区：每个人手边有一摞”待办条”（自己的 dequeue）。谁有活就干自己摞在最上面的那张（后进先出，缓存友好）；一旦自己的摞空了，就去别人摞的最下面抽一张（那张是”最早搁置、包含最多后续工作”的大活）。抽活（steal）要打电话协调，比较贵，所以只在真的没活时才做。
• 架构/机制图解（软件执行模型：continuation stealing + dequeue）：

   cilk_spawn foo(); bar(); cilk_sync;   的实际执行 (work-first: 先跑 child)

   Thread 0 的调用栈                      Thread 0 的 dequeue (deque)
   ┌──────────────────┐                  头(head, 被窃取端)
   │  main 帧          │                  ┌───────────────┐
   │  ...             │                  │ cont: bar()   │ ◄── 刚 push 的续体
   ├──────────────────┤  ── push ──►     └───────────────┘
   │  foo() 帧        │                  尾(tail, 本地 push/pop)
   │  (立即执行 child) │                  ↑ 本地线程在尾部操作: 后进先出 → 深优 → 局部性好
   └──────────────────┘                  ↑ 窃取者在头部操作: 先进先出 → 偷到大块工作
                                         ↓ 头尾分离 ⇒ 两端在不同 cache line ⇒ 冲突小

   ① T0 执行 spawn: 把续体 bar() 压入自己 dequeue 尾部, 立即开始跑 foo()  (continuation stealing)
   ② T1 空闲 (自己队列空) → 随机选一个受害者(如 T0) → 从 T0 队列头部偷
   ③ T1 偷到 "bar()", 压入自己队列尾, 弹出执行
   ④ 双方各自继续; 若 T1 又空, 再偷

   ┌─────────────┐  steal(头部)  ┌─────────────┐
   │ T0 侧 dequeue│◄──────────────│ T1 侧 dequeue│
   │ tail: 本地   │──────────────►│ tail: 本地   │
   └─────────────┘  抢到就搬走    └─────────────┘
       所有 worker: while (有活) { 弹尾部执行; 空则随机窃取 }

   内存布局(伪共享与竞争的关键):
   ┌──────────────┬──────────────┬───────────────────────┐
   │ head (64B)   │ tail (64B)   │ buffers[] (环形数组)   │
   │ 窃取者写/读   │ 本地线程写    │ 本地写尾/窃取读头      │
   └──────────────┴──────────────┴───────────────────────┘
     head 与 tail 必须分处不同 cache line, 否则本地 push 会被窃取者的读
     拖成 cache line 乒乓 (每次 ~50–100 ns, 本该 ~1 ns)

• 关键操作与性能特征：
  • continuation stealing（续体窃取） vs child stealing（子任务窃取）：讲义给出两条抉择——”Run child first: record continuation for later execution”（continuation stealing）与”Run continuation first: record child for later execution”（child stealing）。Cilk 选择先跑 child，理由是讲义明确的两条：(a) 对”循环里逐次 spawn”的代码，调用者线程只需创建 1 个可被偷的条目（代表剩余所有迭代的续体），而 child stealing 会”在开始执行之前就为所有迭代创建了工作”，占 O(N) 空间；(b) 无窃取时，continuation stealing 的执行顺序与”删掉 spawn”的程序完全一致（深度优先调用图遍历）；并且可以证明”work queue 的存储量不超过单线程栈存储量的 T 倍”（T = 线程数）。
  • 窃取粒度的智慧：从 dequeue 的头部窃取 ⇒ 偷到的是调用树中”更靠上、更大”的工作，从而摊薄窃取成本、并与本地线程错开访问区域（减少竞争）。
  • locality：run-child-first + 本地尾操作 ⇒ 本地线程始终工作在自己调用树的局部；讲义总结为”maximizes locality”。
  • T1/T∞ 与 greedy 调度：work-span 模型给出 T_P ≥ max(T1/P, T∞)，而贪心/工作窃取调度器满足 T_P ≤ T1/P + T∞；并行度 = T1/T∞。讲义强调 “parallel slack = ratio of independent work to machine’s parallel execution capability（实践中 ~8 是好比例）”——即可独立工作量最好是机器并行能力的约 8 倍，既够填满机器、又不至于粒度过细。

2.8 Work-span 模型与三种实现机制的统一视角

• 定义与目的：Work（T1）= 全部操作数（单处理器时间）；Span（T∞）= 关键路径长度（无限处理器时间）；Parallelism = T1/T∞。这三个量是判断”该用哪种机制”的统一语言：MPI 关心的是每条 rank 上的 T1 与通信量；OpenMP 关心的是屏障之间的并行度与串行段；Cilk 关心的是 T∞（span）与窃取次数。
• 直观解释（”它是什么？”）：把计算想成做饭：T1 = 所有切菜炒菜的总工时；T∞ = 从”烧水”到”装盘”这条最长依赖链上必须逐分钟等待的时间；并行度 = 总工时 / 关键链长度 = “理论上最多几个厨师有意义”。
• 架构/机制图解（计算 DAG 与关键路径）：

   计算 DAG (每个节点=一个单位工作, 箭头=依赖)
        ┌──────────────────────────────────────────────┐
   T∞ → │ A ──► B ──► C ──────────────► J ──► K ──► L  │ ← 关键路径(span)
        │              \                 ▲           │
        │               \                │           │
        │  D ──► E ──► F ──► G ──► H ─► I            │ ← 旁支工作(可被窃取)
        └──────────────────────────────────────────────┘
   T1 = A..L 的总节点数 = 12
   T∞ = A-B-C-J-K-L      = 6
   并行度 = 12/6 = 2      ← 即使有 64 个核, 也最多 2 倍加速!

   T_P 的界 (greedy / work-stealing 调度器):
     下界:  T_P ≥ max( T1/P , T∞ )
     上界:  T_P ≤ T1/P + T∞
   ⇒ 加速比 S_P = T1/T_P ≤ T1/(T1/P) = P   且   S_P ≤ T1/T∞ = 并行度
   ⇒ 想扩展, 必须先缩短 span (把串行段并行化), 而不是"多 spawn 几下"

   与 Amdahl 定律的关系:
     Amdahl: S ≤ 1/S_serial          (S_serial = 固有串行比例)
     work-span: S ≤ T1/T∞            (span 就是"无法并行的那条链")
     两者是同一件事的两种写法: span 长 ⇒ 串行比例大 ⇒ 扩展性差

• 关键操作与性能特征：对三种实现机制的统一判据是：机制引入的开销必须远小于它带来的并行收益。定量地，用 P 个 worker 一定要满足 T1/P ≫ 单次机制开销（MPI 的一次消息、OpenMP 的一次屏障、Cilk 的一次窃取）。这条判据把第 3 节的三个代码示例串成一条主线。

3. 代码示例与性能分析

3.1 示例一：MPI 五点 Jacobi + halo 交换（非阻塞发送/接收）

/* mpi_jacobi.c —— MPI 2D 五点 Jacobi 迭代 + ghost cell(halo) 交换
 *
 * 编译(必须用 MPI 编译器包装器, 由它决定 -O3 与链接库):
 *     mpicc -O3 -march=native mpi_jacobi.c -o mpi_jacobi -lm
 * 运行:
 *     mpirun -np 8 --bind-to core ./mpi_jacobi 4096 200
 *     (OpenMPI 用 --bind-to core; MPICH 用 -binding core 或 srun --cpu-bind=cores)
 * 说明: -O3 让内层循环向量化; --bind-to core 避免 rank 在核间漂移破坏 NUMA 局部性
 */
#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

int main(int argc, char **argv)
{
    MPI_Init(&argc, &argv);

    int rank, size;
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);   /* 本进程编号 0..P-1 */
    MPI_Comm_size(MPI_COMM_WORLD, &size);   /* 进程总数 P */

    const int    N     = (argc > 1) ? atoi(argv[1]) : 4096;  /* 全局网格 N x N */
    const int    ITERS = (argc > 2) ? atoi(argv[2]) : 200;   /* 迭代上限 */
    const double TOL   = 1e-6;                                /* 收敛阈值 */

    if (N % size != 0) {
        if (rank == 0) fprintf(stderr, "错误: N 必须能被进程数整除\n");
        MPI_Abort(MPI_COMM_WORLD, 1);
    }
    const int    rows  = N / size;                      /* 本地负责的行数 */
    const size_t plane = (size_t)(rows + 2) * N;        /* 含上下两行 ghost */

    double *A = (double *)malloc(plane * sizeof(double));
    double *B = (double *)malloc(plane * sizeof(double));
    if (!A || !B) { fprintf(stderr, "rank %d: 内存不足\n", rank); MPI_Abort(MPI_COMM_WORLD, 2); }

    double *cur = A, *nxt = B;      /* 每个缓冲区都是"2 行 ghost + rows 行有效 + 2 行 ghost" */

    /* ---------- 初始化: 全局边界 = 1.0, 内部 = 0.0 ---------- */
    for (int i = 0; i < rows + 2; ++i) {
        const int gi = rank * rows + (i - 1);        /* 全局行号; ghost 行为 -1 或 N */
        for (int j = 0; j < N; ++j) {
            const int on_bnd = (gi <= 0 || gi >= N - 1 || j == 0 || j == N - 1);
            cur[(size_t)i * N + j] = on_bnd ? 1.0 : 0.0;
        }
    }

    /* 邻居: 用 MPI_PROC_NULL 表示"不存在", 对它的收发自动变成空操作 */
    const int up   = (rank == 0)        ? MPI_PROC_NULL : rank - 1;
    const int down = (rank == size - 1) ? MPI_PROC_NULL : rank + 1;

    MPI_Request req[4];
    const double t0 = MPI_Wtime();
    int iter = 0;

    for (; iter < ITERS; ++iter) {
        double *ghost_top = cur;                              /* 上行 ghost */
        double *my_first  = cur + (size_t)N;                  /* 本地第 1 行 */
        double *my_last   = cur + (size_t)rows * N;           /* 本地最后 1 行 */
        double *ghost_bot = cur + (size_t)(rows + 1) * N;     /* 下行 ghost */

        /* ---------- 步骤 1: 先发起 4 个非阻塞通信, 再等 ---------- */
        int k = 0;
        MPI_Irecv(ghost_top, N, MPI_DOUBLE, up,   0, MPI_COMM_WORLD, &req[k++]);
        MPI_Irecv(ghost_bot, N, MPI_DOUBLE, down, 1, MPI_COMM_WORLD, &req[k++]);
        MPI_Isend(my_first, N, MPI_DOUBLE, up,   1, MPI_COMM_WORLD, &req[k++]);
        MPI_Isend(my_last,  N, MPI_DOUBLE, down, 0, MPI_COMM_WORLD, &req[k++]);
        MPI_Waitall(k, req, MPI_STATUSES_IGNORE);

        /* ---------- 步骤 2: 计算本地内点(整行 bulk 通信已完成) ---------- */
        double mydiff = 0.0;
        for (int i = 0; i < rows; ++i) {
            const double *up_r   = cur + (size_t)i * N;
            const double *mid_r  = cur + (size_t)(i + 1) * N;
            const double *down_r = cur + (size_t)(i + 2) * N;
            double       *out    = nxt + (size_t)(i + 1) * N;

            out[0] = mid_r[0];              /* 全局左右边界固定 */
            out[N - 1] = mid_r[N - 1];
            double d = 0.0;
            for (int j = 1; j < N - 1; ++j) {
                const double v = 0.25 * (mid_r[j - 1] + mid_r[j + 1] + up_r[j] + down_r[j]);
                d += fabs(v - mid_r[j]);    /* 局部累加, 循环外不碰共享变量 */
                out[j] = v;
            }
            mydiff += d;
        }
        /* 保持 ghost 行在缓冲区交换后仍有效(边界 rank 的 ghost 恒为边界值) */
        for (int j = 0; j < N; ++j) {
            nxt[j]                          = cur[j];
            nxt[(size_t)(rows + 1) * N + j] = cur[(size_t)(rows + 1) * N + j];
        }
        double *tmp = cur; cur = nxt; nxt = tmp;   /* 双缓冲交换 */

        /* ---------- 步骤 3: 全局收敛判定(所有 rank 必须得到同一结论) ---------- */
        double gdiff = 0.0;
        MPI_Allreduce(&mydiff, &gdiff, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
        if (gdiff / ((double)N * (double)N) < TOL) { ++iter; break; }
    }
    const double t1 = MPI_Wtime();

    /* ---------- 校验和: 只做一次, 用于结果正确性检查 ---------- */
    double local_sum = 0.0;
    for (int i = 1; i <= rows; ++i)
        for (int j = 0; j < N; ++j) local_sum += cur[(size_t)i * N + j];
    double gsum = 0.0;
    MPI_Reduce(&local_sum, &gsum, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);

    if (rank == 0) {
        const double per_it = (t1 - t0) / (iter > 0 ? iter : 1);
        const double bytes  = (double)(N - 2) * (double)(N - 2) * sizeof(double) * 2.0;
        printf("P=%d N=%d iters=%d 总时间=%.4f s  每迭代=%.3f ms  "
               "聚合有效带宽=%.1f GB/s  校验和=%.6f\n",
               size, N, iter, t1 - t0, per_it * 1e3, bytes / per_it / 1e9, gsum);
    }
    free(A); free(B);
    MPI_Finalize();
    return 0;
}

【代码做什么？】

1. MPI_Init 建立 rank/size，随后做1D 行块域分解：全局 N×N 网格按行切成 size 块，每块 rows = N/size 行，本地数组额外带上下各 1 行 ghost。
2. 初始化时把全局边界（gi<=0 \|\| gi>=N-1 \|\| j==0 \|\| j==N-1）置 1.0，其余置 0.0——这就是 Dirichlet 边界条件的固定值。
3. 每轮迭代分三段：(a) 4 个非阻塞通信（Irecv 两行 ghost、Isend 本地首末行），最后一次 Waitall；(b) 计算所有本地内点的 5 点平均并累加局部 mydiff；(c) Allreduce 求和后所有 rank 用同一判据决定是否退出。
4. 双缓冲 cur/nxt 交换避免写覆盖，最后用 MPI_Reduce 汇总校验和，rank 0 打印每迭代耗时与聚合有效带宽。

【并行机制与性能解说】

• 进程如何创建、工作如何分配：mpirun -np P 由 MPI launcher（orted/hydra）在节点上 fork P 个进程；没有线程创建，每个 rank 是一个独立地址空间，靠固定映射 rank ↔ 行区间 完成 assignment，属于纯静态分配（零运行时调度开销）。
• 共享数据如何处理：没有共享数据。每 rank 只有自己的 rows+2 行；跨 rank 的依赖靠显式复制（ghost row）解决；通信是整行 bulk 传输（N 个 double = 32 KB @ N=4096）。
• Work / Span / 并行度：设每点代价 c_pt（≈ 20 B 内存流量 + 5 flops，见第 4 节）。
  • 每迭代每 rank 的 Work：W = rows·(N−2)·c_pt（≈ 4096·512 = 2.1M 点 @ P=8）。
  • 每迭代关键路径：D = c_pt + 2·α + N·4/BW_net + log₂P·(α + 8/BW_net)。因为同一次迭代内所有内点互相独立（只读上一轮的 cur），点更新本身的 span 是 O(1)（用无限处理器），通信的 span 是”消息延迟 + 归约树深度”。
  • 并行度 = W_total/D_total：整体 T1 = ITERS·N²·c_pt，T∞ = ITERS·(α + N·4/BW + log₂P·α)。取 α = 2 µs、BW = 6 GB/s、N = 4096、ITERS = 200：T∞ ≈ 200 × (2 + 2.7 + 2·2) µs ≈ 1.7 ms，而 T1 ≈ 200 × 16.8M × 20 B / 15 GB/s ≈ 4.5 s ⇒ 并行度 ≈ 2600，远大于 P = 128。结论：MPI 场景下 span 通常不是瓶颈，瓶颈是带宽与延迟的开销项。
• 瓶颈：(1) 内存带宽——5 点 stencil 的算术强度只有 0.25 flop/byte（严重 memory-bound）；(2) 通信/计算比 ≈ 0.4/rows，rows 小时通信占比飙升（P 增大而问题不变时）；(3) Allreduce 频率——每迭代一次全局集合操作，是小消息延迟（α）的固定税；(4) 负载不均——若 N 不被 P 整除则尾部 rank 空闲。

3.2 示例二：OpenMP 调度策略、reduction 与伪共享对照

/* omp_scaling.c —— 调度策略对比(reduction) + 伪共享对照
 *
 * 编译: gcc -O3 -fopenmp -march=native omp_scaling.c -o omp_scaling -lm
 *       (ICC/Clang: icx -O3 -qopenmp ... / clang -O3 -fopenmp ...)
 * 运行: OMP_NUM_THREADS=8 OMP_PROC_BIND=close OMP_PLACES=cores ./omp_scaling 20000000
 *       (绑定线程到物理核, 避免迁移与跨 NUMA 访存)
 */
#define _GNU_SOURCE
#include <omp.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

/* 每次迭代的"真实工作": 混合 sqrt/sin/cos, 让编译器无法把它优化成常量 */
static inline double kernel(long i)
{
    const double x = (double)(i & 1023) * 1e-3;
    return sqrt(x + 1.0) * sin(x) + cos(x) * 0.5;
}

/* ---------- (1) 三种调度策略 ---------- */
static double bench_sched(int mode, long chunk, long N)
{
    double sum = 0.0;
    const double t0 = omp_get_wtime();

    if (mode == 0) {                     /* static: 进入并行区时一次性切块 */
        #pragma omp parallel for schedule(static, chunk) reduction(+:sum)
        for (long i = 0; i < N; ++i) sum += kernel(i);
    } else if (mode == 1) {              /* dynamic: 每 chunk 一次原子领号 */
        #pragma omp parallel for schedule(dynamic, chunk) reduction(+:sum)
        for (long i = 0; i < N; ++i) sum += kernel(i);
    } else {                             /* guided: chunk 递减, 不小于 chunk */
        #pragma omp parallel for schedule(guided, chunk) reduction(+:sum)
        for (long i = 0; i < N; ++i) sum += kernel(i);
    }

    const double dt = omp_get_wtime() - t0;
    if (sum < -1e300) printf("never\n");   /* 防止 sum 被优化掉 */
    return dt;
}

/* ---------- (2) 伪共享对照: 只有 padding 不同 ---------- */
struct padded { long long v; char pad[64 - sizeof(long long)]; };   /* 每计数器独占 64B */

static double bench_falsesharing(long iters, int use_padding)
{
    const int    T      = omp_get_max_threads();
    const size_t stride = use_padding ? sizeof(struct padded) : sizeof(long long);
    char *buf = (char *)calloc((size_t)T * stride + 64, 1);
    long long *base = (long long *)buf;          /* 与 padding 无关, 都按 long long 访问 */
    long long **slot = (long long **)malloc((size_t)T * sizeof(long long *));
    for (int t = 0; t < T; ++t) slot[t] = (long long *)(buf + (size_t)t * stride);

    const double t0 = omp_get_wtime();
    #pragma omp parallel
    {
        const int tid = omp_get_thread_num();
        long long *mine = slot[tid];
        for (long i = 0; i < iters; ++i) mine[0] += 1;   /* 独占: L1 命中; 共享: 乒乓 */
    }
    const double dt = omp_get_wtime() - t0;

    free(slot); free(buf);
    return dt;
}

int main(int argc, char **argv)
{
    const long N = (argc > 1) ? atol(argv[1]) : 20000000L;
    const long CH[] = { 1, 64, 4096, 262144 };   /* 四档 chunk, 观察开销/均衡权衡 */
    const int  NC   = (int)(sizeof(CH) / sizeof(CH[0]));

    printf("== 调度策略 (线程=%d, N=%ld) ==\n", omp_get_max_threads(), N);
    printf("%-22s %10s %12s\n", "schedule", "time(s)", "Miter/s");
    for (int m = 0; m < 3; ++m) {
        const char *nm = (m == 0) ? "static" : (m == 1) ? "dynamic" : "guided";
        for (int c = 0; c < NC; ++c) {
            double best = 1e30;
            for (int rep = 0; rep < 3; ++rep) {           /* 取 3 次最小值, 降噪 */
                const double dt = bench_sched(m, CH[c], N);
                if (dt < best) best = dt;
            }
            char label[64];
            snprintf(label, sizeof(label), "%s,chunk=%ld", nm, CH[c]);
            printf("%-22s %10.4f %12.1f\n", label, best, (double)N / best / 1e6);
        }
    }

    printf("\n== 伪共享对照 (线程=%d, 每线程 3e7 次自增) ==\n", omp_get_max_threads());
    const long FS = 30000000L;
    const double t_share = bench_falsesharing(FS, 0);
    const double t_pad   = bench_falsesharing(FS, 1);
    printf("共享 cache line : %.4f s (%.1f Mop/s)\n", t_share, (double)FS / t_share / 1e6);
    printf("64B 隔离        : %.4f s (%.1f Mop/s)\n", t_pad,   (double)FS / t_pad   / 1e6);
    printf("伪共享惩罚       : %.1fx\n", t_share / t_pad);
    return 0;
}

【代码做什么？】

1. bench_sched 对同一个循环分别用 schedule(static\|dynamic\|guided, chunk) 执行，reduction(+:sum) 让运行时自动分配每线程私有累加器并在结束时合并；外层取 3 次最小值降噪。
2. bench_falsesharing 分配 T 个 long long 计数器，唯一区别是 stride：无 padding 时 T 个计数器挤在同一条 64 字节 cache line 内，有 padding 时每个独占一条 cache line；每个线程只更新”自己”的计数器。
3. main 打印调度策略的时间/吞吐表，以及伪共享惩罚倍数（理论上无 padding 版本会慢几十倍）。

【并行机制与性能解说】

• 线程如何创建、工作如何分配：omp_get_max_threads() 个 worker 在首次并行区被懒创建并常驻（GCC/Clang 运行时默认在屏障处自旋/休眠）；parallel for 由编译器轮廓化成一个带 (lower, upper, stride) 参数的函数，每个线程调用它并只跑自己那段。static 的切分发生在进入并行区时（纯索引运算，无原子）；dynamic/guided 每个 chunk 需要一次原子领号（libgomp 用 __atomic_fetch_add 或带锁的回退路径）。
• 共享数据如何处理：sum 是 reduction 变量 → 每线程一份私有副本（栈/线程本地），循环结束由运行时按树形合并；slot[tid] 是每线程独占写 → 但若落在同一 cache line，写操作要求独占所有权（MESI 的 M 态），每写一次就把别的核的副本置为 Invalid，造成 cache line 在核间来回搬运。
• Work / Span / 并行度：
  • Work：W = N·c_iter（+ 调度开销）。Span：D ≈ (N/c)·t_sync + c·t_iter + log₂P·t_merge，其中 (N/c)·t_sync 是领号临界区的串行链（这是最容易忽略的串行段），c·t_iter 是动态调度固有的最大不均衡（最后一个 chunk），log₂P·t_merge 是 reduction 的合并树。
  • 并行度 = W/D = N·c_iter / ((N/c)·t_sync + c·t_iter + O(log P))。
  • 数值例：N = 2×10⁷、t_iter = 8 ns、t_sync = 25 ns、P = 8。
    • chunk = 1：D ≈ 2×10⁷ × 25 ns = 0.5 s，而 W = 0.16 s ⇒ 并行度 < 1，比串行还慢（领号临界区成为串行瓶颈，正是讲义中”fine granularity: potential for high synchronization cost / it’s serial execution (recall Amdahl’s Law)”）。
    • chunk = 4096：D ≈ (2×10⁷/4096)×25 ns + 4096×8 ns = 0.122 ms + 0.033 ms ≈ 0.155 ms。
    • chunk = 262144：D ≈ 76×25 ns + 2.1 ms ≈ 2.1 ms（不均衡项主导）。
    • 由 c* = sqrt(N·t_sync/t_iter) = sqrt(2×10⁷×25/8) ≈ 7.9×10³ 可知 chunk≈4096–8192 最优——这就是”均匀代价循环用 static、极不均匀代价用 dynamic + 中等 chunk（或 guided）”的定量依据。
  • 伪共享部分：W = T·FS 次自增，Span = FS·t_line_pingpong（每个计数器的更新链在共享 cache line 上被串行化）。FS = 3×10⁷：L1 命中自增约 1–1.5 ns（独占时 Span ≈ 30–45 ms）；跨核乒乓约 40–100 ns ⇒ Span ≈ 1.2–3 s。因此即使 Work 相同，Span 相差 50–100 倍，扩展性直接崩塌——这是”Work/Span 决定并行度”的最直观演示。
• 瓶颈：(1) 领号临界区的串行化 + cache line 争抢；(2) 伪共享；(3) 线程绑核失败的NUMA 远端访问（延迟约 1.5–2×，带宽约 1/2）；(4) 每次并行区的 fork-join（约 1–10 µs）在小循环上主导。

3.3 示例三：Cilk / OpenCilk —— cilk_for、并行快排与 work stealing 的观测

/* cilksort.cpp —— Cilk/OpenCilk: cilk_for 映射、并行快速排序、worker 数观测
 *
 * 编译(OpenCilk, 推荐):
 *     clang++ -fopencilk -O3 -march=native cilksort.cpp -o cilksort -lm
 * 编译(传统 Intel Cilk Plus / GCC <= 7; GCC 8 起已移除 -fcilkplus):
 *     g++ -fcilkplus -O3 cilksort.cpp -o cilksort -lm
 * 串行对照: g++ -O3 -DCILK_SERIAL cilksort.cpp -o cilksort_serial -lm
 * 运行: CILK_NWORKERS=8 ./cilksort 100000000
 */
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <chrono>
#include <algorithm>

#if !defined(CILK_SERIAL)
#  include <cilk/cilk.h>
#  include <cilk/cilk_api.h>
#  if defined(__has_include)
#    if __has_include(<cilk/reducer_opadd.h>)
#      define HAVE_CILK_REDUCER 1
#      include <cilk/reducer_opadd.h>
#    endif
#  endif
#endif

#define CUTOFF 4096      /* 小于该规模就串行排序: 避免 spawn 开销压过并行收益 */

static void insertion_sort(int *a, int lo, int hi)
{
    for (int i = lo + 1; i < hi; ++i) {
        const int key = a[i];
        int j = i - 1;
        while (j >= lo && a[j] > key) { a[j + 1] = a[j]; --j; }
        a[j + 1] = key;
    }
}

/* Lomuto 划分, 返回 pivot 的最终下标 */
static int partition_lomuto(int *a, int lo, int hi)
{
    const int pivot = a[hi - 1];
    int i = lo;
    for (int j = lo; j < hi - 1; ++j)
        if (a[j] < pivot) { std::swap(a[i], a[j]); ++i; }
    std::swap(a[i], a[hi - 1]);
    return i;
}

/* ---------- 串行版本: 用于测 T1 (work) ---------- */
static void qsort_serial(int *a, int lo, int hi)
{
    while (hi - lo > CUTOFF) {
        const int p = partition_lomuto(a, lo, hi);
        qsort_serial(a, lo, p);
        lo = p + 1;                       /* 尾递归改成循环, 只对左半递归 */
    }
    insertion_sort(a, lo, hi);
}

/* ---------- 并行版本: 先划分, spawn 左半, 本地继续右半 ---------- */
static void qsort_parallel(int *a, int lo, int hi)
{
    if (hi - lo <= CUTOFF) { insertion_sort(a, lo, hi); return; }
    const int p = partition_lomuto(a, lo, hi);
    cilk_spawn qsort_parallel(a, lo, p);      /* 把左半交给"可被窃取"的 child */
    qsort_parallel(a, p + 1, hi);             /* 本地继续右半(continuation 语义) */
    cilk_sync;                                /* 函数末尾本就有隐式 sync, 写出更清晰 */
}

/* ---------- 每 worker 独占 cache line 的并行求和 ---------- */
static long long checksum_padded(const int *a, long n)
{
    int W = 1;
#if !defined(CILK_SERIAL)
    W = __cilkrts_get_nworkers();
#endif
    long long *partial = nullptr;
    if (posix_memalign((void **)&partial, 64, (size_t)W * 64) != 0) exit(1);
    memset(partial, 0, (size_t)W * 64);

#if defined(CILK_SERIAL)
    for (long i = 0; i < n; ++i) partial[0] += a[i];
#else
    cilk_for (long i = 0; i < n; ++i) {
        const int w = __cilkrts_get_worker_number();
        partial[(size_t)w * 8] += a[i];       /* 每 worker 间隔 64B: 无伪共享 */
    }
#endif
    long long s = 0;
    for (int w = 0; w < W; ++w) s += partial[(size_t)w * 8];
    free(partial);
    return s;
}

#if defined(HAVE_CILK_REDUCER)
/* 惯用写法: reducer(hyperobject) 自动为每个 worker 维护视图, 窃取时合并 */
static long long checksum_reducer(const int *a, long n)
{
    cilk::reducer_opadd<long long> s(0);
    cilk_for (long i = 0; i < n; ++i) s += a[i];
    return s.get_value();
}
#endif

/* 讲义中的 recursive_for: 折半产生任务, 是 cilk_for 展开后的典型形态 */
static void recursive_for(int *a, int begin, int end, long long *local)
{
    while (end - begin > 1024) {
        const int mid = begin + (end - begin) / 2;
#if !defined(CILK_SERIAL)
        cilk_spawn recursive_for(a, begin, mid, local);
#endif
        begin = mid;                       /* 本地继续后半段 */
    }
    long long s = 0;
    for (int i = begin; i < end; ++i) s += a[i];
    *local += s;
}

static double now_sec()
{
    using namespace std::chrono;
    return duration<double>(steady_clock::now().time_since_epoch()).count();
}

int main(int argc, char **argv)
{
    const long n = (argc > 1) ? atol(argv[1]) : 100000000L;

    int *a = (int *)malloc((size_t)n * sizeof(int));
    int *b = (int *)malloc((size_t)n * sizeof(int));
    if (!a || !b) { fprintf(stderr, "内存不足\n"); return 1; }

    srand(12345);
    for (long i = 0; i < n; ++i) a[i] = rand();      /* 随机数据 ⇒ 划分接近平衡 */
    memcpy(b, a, (size_t)n * sizeof(int));

    /* --- 串行基准 (T1) --- */
    double t0 = now_sec();
    qsort_serial(b, 0, (int)n);
    double t1 = now_sec();
    const double T1 = t1 - t0;

    /* --- 并行版本 (T_P) --- */
    t0 = now_sec();
#if defined(CILK_SERIAL)
    qsort_serial(a, 0, (int)n);
#else
    qsort_parallel(a, 0, (int)n);
#endif
    t1 = now_sec();
    const double TP = t1 - t0;

    if (memcmp(a, b, (size_t)n * sizeof(int)) != 0) { fprintf(stderr, "结果不一致!\n"); return 2; }

    int nw = 1;
#if !defined(CILK_SERIAL)
    nw = __cilkrts_get_nworkers();
#endif
    const long long sum = checksum_padded(a, n);
#if defined(HAVE_CILK_REDUCER)
    const long long sum2 = checksum_reducer(a, n);
#else
    const long long sum2 = sum;
#endif

    printf("workers=%d  n=%ld\n", nw, n);
    printf("串行 T1   = %7.3f s\n", T1);
    printf("并行 T_P  = %7.3f s   (加速比 %.2fx, 效率 %.1f%%)\n",
           TP, T1 / TP, 100.0 * T1 / TP / nw);
    printf("校验: memcmp=OK, sum=%lld, sum(reducer)=%lld\n", sum, sum2);
    free(a); free(b);
    return 0;
}

【代码做什么？】

1. qsort_serial 是同一算法去掉 spawn 的版本，用来测 T1（Work）；尾递归改成循环只对左半递归，避免深栈（讲义中的 “run continuation first” 思路在串行代码里的对应物）。
2. qsort_parallel 先划分，然后 cilk_spawn 左半、本地继续右半——这正是讲义描述的实现选择：”Run child first: record continuation for later execution”（continuation stealing）：child 立即执行，续体压入本地 dequeue 供别人窃取。
3. checksum_padded 用 cilk_for 做数据并行映射，并用”每 worker 独占 64 字节”的局部累加器避免伪共享；checksum_reducer 展示惯用的 reducer（hyperobject）写法（若编译器提供 reducer 头文件）。
4. recursive_for 是讲义中 cilk_for 折半展开的形态：不断二分并把前半段 spawn 出去、本地继续后半段。
5. 最后比较串行/并行结果是否一致（memcmp），打印加速比与效率。

【并行机制与性能解说】

• worker 如何创建、工作如何分配：Cilk 运行时在首次 spawn 时懒创建一个 worker 线程池，大小 ≈ 机器执行上下文数（CILK_NWORKERS 可覆盖）。cilk_spawn 在编译后是”把续体 push 到本 worker 的 dequeue 尾部 + 跳转到 child 代码”的少量指令；只有真的发生窃取才会访问其他 worker 的 dequeue。cilk_for 被编译器变换为递归折半 + spawn/本地继续（类似上面的 recursive_for），因此它天然产生 O(log N) 层、每层可窃取的工作。
• 共享数据如何处理：共享数组 a[] 就地划分（in-place partition），因此没有显式通信；同步完全由 spawn/sync 表达；reducer 通过”每 worker 一个视图（view）+ 窃取/加入时合并”避免锁。
• Work / Span / 并行度：
  • Work（T1）：平衡划分下 T1 = Θ(n log n) 次比较。n = 10⁸ ⇒ 约 2.7×10⁹ 次操作。
  • Span（T∞）：由于 partition_lomuto 本身是串行的，且每次只 spawn 左半、本地继续右半，关键路径是”每一层划分 + 最后一次串行排序”：T∞ = Θ(n + g log g) = Θ(n)（主导项是顶层那次 O(n) 划分不断累加在右脊上）。
  • 并行度 = T1/T∞ = Θ(log n)。n = 10⁸ ⇒ log₂ n ≈ 27 ⇒ 并行度 ≈ 27。也就是说，即使有 64 个核，这个快排的”理论可用核数”也只有约 27（greedy 界 T_P ≤ T1/P + T∞ 在此变成 T_P ≈ T1/27 + T∞，P 超过 27 后几乎收益归零）。
  • 怎么改善：把 span 缩短——并行化 partition（并行划分 span 可到 O(log n)），并对两半都 spawn，则 T∞ = O(log² n)，并行度升到 Θ(n/log n)（10⁸/27 ≈ 3.7×10⁶）。这正是”扩展性上限由 span 而非 spawn 数量决定“的最好例证。
• 瓶颈：(1) span（串行 partition 链）——上一条已量化；(2) 窃取成本——每次窃取要跨核访问 victim 的 dequeue 头部（0.1–5 µs），若 cutoff 太小（spawn 过密），窃取/推送次数暴涨；(3) recursive_for 与 spawn-per-iteration 的差异——讲义指出，在循环里逐次 cilk_spawn 会让”剩余迭代的续体”在每次迭代后都变成可窃取项，因此更快把机器填满；而折半递归需要深入若干层才有足够可窃取的工作，机器填充更慢（但队列条目更少、粒度更粗）；(4) 伪共享——若图省事用 partial[worker] 而不做 64 字节隔离，求和会被 cache line 乒乓主导。

4. 性能模型与复杂度分析

4.1 四个定量工具

1. Amdahl 定律：设固有串行比例 S，则 S_P ≤ 1/S。讲义给出的图像是 S = 0.01/0.05/0.1 对应的上限 100/20/10。关键推论：哪怕只有 5% 的串行段，P = 32 时的上限也只有约 17×。
2. Work-Span（T1/T∞）：T_P ≥ max(T1/P, T∞)，贪心/工作窃取调度器 T_P ≤ T1/P + T∞；并行度 = T1/T∞；加速比上限 = min(P, T1/T∞)。
3. 算术强度与 Roofline：AI = FLOPs / Bytes；可达性能 P ≤ min(峰值算力, AI × 内存带宽)。
4. 延迟-带宽模型：t_msg = α + n/BW（α = 启动延迟，n = 消息字节数，BW = 有效带宽）；集合通信 t_coll ≈ log₂P·α + (P−1)/P·n/BW。

表 4：典型延迟/带宽量级（用于定量估算，来自公开硬件与实现资料）

4.2 数值算例 A：MPI stencil 的性能上界（与 3.1 代码对应）

假设参数

步骤 1：算术强度与 Roofline 判定

• AI = 5 flops / 20 B = 0.25 flop/byte。
• 单核机器平衡点 = 48 GFLOPS / 12.5 GB/s = 3.84 flop/byte。
• AI (0.25) ≪ 3.84 ⇒ 严重 memory-bound；单核可达性能 = 0.25 × 12.5 GB/s = 3.1 GFLOPS，仅为峰值的 6.5%。
• 全机聚合上界 = 0.25 × (128 × 12.5 GB/s) = 400 GFLOPS，而全机峰值 = 128 × 48 = 6144 GFLOPS ⇒ 无论怎么优化，这个 stencil 也只能跑到约 400 GFLOPS（6.5%）。

步骤 2：每迭代每 rank 的计算时间

• 每 rank 内点 = 8192 × 64 = 524288 点。
• 内存流量 = 524288 × 20 B = 10.49 MB。
• t_comp = 10.49 MB / 12.5 GB/s = 839 µs（若实际有效流量按 8 B/点计则约 336 µs，这里取保守的 839 µs）。

步骤 3：每迭代每 rank 的通信时间

• halo 数据 = 2 行 × 8192 float × 4 B = 65.5 KB。
• 传输时间 = 65.5 KB / 6 GB/s = 10.9 µs；加上 2 条消息的 α = 2 × 2 µs = 4 µs；MPI_Allreduce（约 2·log₂128 = 14 步 × 2 µs ≈ 28 µs，实际实现用树形会低一些）。
• t_comm ≈ 10.9 + 4 = 14.9 µs（halo，可与计算重叠/接近）；Allreduce ≈ 10–28 µs。

步骤 4：结论

• 每迭代 t_iter ≈ 839 + 15 + 20 ≈ 874 µs，其中通信约占 4%。
• 表面-体积比公式：通信字节 / 计算字节 = 0.4/rows = 0.4/64 = 0.00625；换成时间来算，t_comm/t_comp ≈ 14.9/839 = 1.8%，量级一致。
• 强扩展（fixed problem, P↑）：P 从 128 → 1024 时 rows = 8，通信占比按 1/rows 放大 8 倍 ⇒ 通信从 ~4% 涨到 ~20%，剩余可扩展性被通信吃掉。
• 弱扩展（fixed rows/rank, P↑）：rows 恒为 64，通信占比恒定 ⇒ 弱扩展效率接近 100%（这正是 stencil 应用在超算上都用弱扩展报告的原因）。
• 提升方向：因为 memory-bound，唯一有效手段是降低每点内存流量（时间分块/缓存分块：沿 N 方向分成条带，使 5 个邻点尽量命中 L1/L2），否则加核无用。这就是 Roofline 的直接指导。

   Roofline (单核: 峰值 48 GFLOPS, 带宽 12.5 GB/s)

   GFLOPS
    48 ├───────────────────────────────●════════════  ← 峰值算力天花板 (compute bound)
       │                            ╱
       │                         ╱  ← 斜线斜率 = 带宽 12.5 GB/s
       │                      ╱
       │                   ╱
       │                ╱
       │             ╱
       │          ╱
       │       ╱
       │    ╱
     3.1 ├─● (AI=0.25, stencil: 5 flops / 20 B)          ← memory bound 区域
       │
       └────┬────────┬─────────┬─────────┬───────────►  AI (flop/byte)
            0.25     1        3.84       8
                        ↑ 机器平衡点 = 48/12.5 = 3.84 ⇒ 只有 AI > 3.84 才可能 compute bound

4.3 数值算例 B：OpenMP 动态调度的最优 chunk 与并行度上限（与 3.2 代码对应）

假设参数：N = 10⁹ 次迭代；t_iter = 1 ns（每次迭代算一点工作）；动态领号临界区 t_sync = 20 ns（原子 + cache line 传输）；P = 64 线程。

• 总 Work：T1 = 10⁹ × 1 ns = 1.0 s。
• 串行段（领号）总时间：T_serial = (N/c) × t_sync。
• 最大不均衡：T_imb = c × t_iter（最后一个 chunk 可能只有一个线程在跑）。
• 目标函数（关键洞察：领号是串行的，按 Amdahl 直接相加）： T_P ≈ T1/P + T_serial + T_imb = T1/P + (N/c)·t_sync + c·t_iter
• 最优 chunk（对 c 求导置零）：c* = sqrt(N·t_sync / t_iter) = sqrt(10⁹ × 20 / 1) = sqrt(2×10¹⁰) ≈ 1.41×10⁵。
• 最优开销：2·sqrt(N·t_sync·t_iter) = 2 × 1.41×10⁵ ns = 283 µs = 0.283 ms。

结论：

1. schedule(dynamic,1) 在 10⁹ 次迭代下把程序变成 20 秒的串行程序——这正是讲义反复强调的”细粒度划分的同步开销本身就是串行执行（recall Amdahl’s Law）”。
2. chunk 的选取遵循 c* = sqrt(N·t_sync/t_iter)；schedule(guided) 用递减 chunk 逼近这个折中，不需要手工调参。
3. 若每次迭代代价均匀，static 是零开销且 100% 效率的最优解 ⇒ 先试 static，再按需升级到 dynamic/guided（讲义 TIP #1：”Always implement the simplest solution first, then measure performance”）。

4.4 数值算例 C：Cilk 的并行度上限与 parallel slack（与 3.3 代码对应）

• 假设：n = 10⁸ 个 int；比较/交换总代价 T1 = n·log₂n·t_cmp = 10⁸ × 27 × 1 ns = 2.7 s；span T∞ = n·t_part = 10⁸ × 1 ns = 0.1 s（顶层划分链主导）。
• 并行度 = T1/T∞ = 27。
• greedy 界给出的加速比上限（P = 64）： T_64 ≤ T1/64 + T∞ = 0.042 + 0.1 = 0.142 s ⇒ S ≤ 2.7/0.142 ≈ 19×（效率 30%）。
• 若 P = 27：T_27 ≈ 0.1 + 0.1 = 0.2 s ⇒ S ≈ 13.5×（效率 50%）。
• 若把 partition 并行化（span → O(log²n)·t ≈ 27² × 1 ns ≈ 0.73 µs）并对两半都 spawn： T∞' ≈ 0.73 µs，并行度 ≈ 3.7×10⁶ ⇒ 在 P = 64 时受限于 P 本身，T_64 ≈ T1/64 + 1 µs ≈ 42 ms ⇒ S ≈ 64×（效率接近 100%）。
• parallel slack 视角（讲义准则 ≈ 8）：T1/T∞ = 27 对 P = 8 的机器来说 slack = 27/8 ≈ 3.4（偏低，容易受调度抖动影响）；把 cutoff g 调小能提高 slack（任务更多），但窃取次数随之增加：每次窃取的代价 0.1–5 µs，若窃取次数达到 10⁵ 量级，额外开销就是 0.01–0.5 s，足以吃掉加速收益。因此 cutoff 的调优目标是”让每次窃取的收益（被偷走的计算量）≫ 窃取成本”：让被偷走的子树至少包含约 10–100 µs 的计算。

表 5：三种”任务分配/调度机制”的定量对比（P = 64，每次任务代价 1 µs 的情形）

5. 关键要点

1. 三种系统的差别本质是”责任归属”，不是”能力”。MPI 把分解/分配/编排/映射几乎全交给程序员（控制力最强、第一条正确程序最难）；OpenMP 由编译器 outlining + 运行时线程池接管分配与编排（schedule 子句是唯一的旋钮）；Cilk 把分配与映射全部交给 work-stealing 运行时。选择时先问”这个问题的瓶颈是通信、同步频率还是 span”，再选系统，而不是反过来。

2. 机制开销必须被每次迭代的计算量摊薄，这可以用一条不等式检验：要让 P 个 worker 有收益，必须满足 T1/P ≫ 单次机制开销（MPI 一次消息 α、OpenMP 一次屏障/fork-join、Cilk 一次窃取）。算例：10⁹ 次迭代用 schedule(dynamic,1) 时，20 ns 的领号被串行化 10⁹ 次 = 20 s，程序比串行还慢。

3. 扩展性的天花板是 span（关键路径），不是 spawn / 任务的数量。greedy 界 T_P ≤ T1/P + T∞ 说明只要 T∞ 不小，加核就无用；快排的 T∞ = Θ(n) 使并行度只有 Θ(log n) ≈ 27（n = 10⁸）。先并行化串行瓶颈（例如 partition），再增加并行度。

4. 数据搬运与同步的”面积”决定实际性能：MPI 的通信/计算比 ≈ 0.4/rows（表面-体积比 ⇒ 弱扩展友好、强扩展受限）；OpenMP 的伪共享把本该 1 ns 的 L1 自增变成 30–100 ns 的 cache line 乒乓（50–100× 惩罚）；Cilk 的 dequeue 头尾必须分处不同 cache line，否则本地 push 被窃取者拖慢。

5. 先测后优，机制按需升级：讲义 TIP #1——”Always implement the simplest solution first, then measure performance”。工程路径通常是：static/blocked → 发现负载不均 → dynamic/guided 或 work stealing → 发现同步/窃取开销 → 调整 chunk/cutoff（c* = sqrt(N·t_sync/t_iter)，cutoff 让每次窃取收益 ≫ 窃取成本），并用 OMP_PROC_BIND/--bind-to core 保证 NUMA 局部性。

6. 常见陷阱与注意事项

• 用同步阻塞 Send/Recv 做环形 halo 交换会死锁：所有 rank 都先 Send 后 Recv，而同步 Send 要等接收方确认，于是全员互相等待——这正是讲义中”big problem with our message passing solver if it uses synchronous send/recv”的场景。修复方式有三：(a) 奇偶错开（偶数 rank 先发后收、奇数 rank 先收后发，即讲义中的修正版）；(b) 用配对的 MPI_Sendrecv；(c) 用非阻塞 Isend/Irecv + Waitall（最推荐，同时还能重叠通信与计算）。注意方向搞错（把”我先发”当成”我先收”）只是把死锁转移到另一侧。
• 把 OpenMP 的”隐式屏障”当成免费的：parallel for 末尾、single、critical 的隐含语义都可能插入屏障，而屏障把”最慢的线程”变成相位长度。常见后果：一个循环里嵌套两次并行区 = 两次 fork-join + 两次屏障；用 nowait 去屏障却忘了随后的数据依赖，直接引入数据竞争。辨析：critical 是互斥（保正确性），barrier 是相位同步（保顺序），atomic 只保证单个内存操作的原子性（不保证复合语句的原子性）。
• 伪共享（false sharing）把”每线程私有变量”变成通信：int partial[MAX_THREADS] 只要两个元素落在同一条 64 字节 cache line 内，两个核心就会反复争夺该行的独占权。修复：按 64 字节对齐/填充（示例 3.2 中的 struct padded），或用 reduction / reducer / 线程局部变量后一次性合并。注意 alignas(64) 只解决”跨对象”共享，若数组元素密集仍需 stride 填充。
• 忽略内存带宽与 Roofline，只盯着 FLOPS 和核心数：5 点 stencil 的算术强度只有 0.25 flop/byte，机器平衡点约 3.84 flop/byte，因此无论多少核都只能拿到峰值算力的约 6.5%；MPI 场景下更常见的错误是”P 翻倍就期望时间减半”，而实际被 0.4/rows 的通信占比与网络对分带宽卡住。先算 AI，再判断是 compute bound 还是 memory bound，然后才决定优化方向。
• 任务粒度过细（过度 spawn / chunk=1）：Cilk 里对每个循环迭代都 cilk_spawn 会让 dequeue 条目与窃取次数暴涨（窃取 0.1–5 µs 一次）；OpenMP 里 schedule(dynamic,1) 让计数器临界区成为串行段；MPI 里逐元素 Send 让小消息延迟 α 主导。准则：单个任务的执行时间应至少是机制开销的 10–100 倍（Cilk 的 CUTOFF、OpenMP 的 chunk、MPI 的整行 bulk 通信都是这条准则的实例）。
• 内存序与可见性假设错误：讲义中的 flag 示例（x = 1; flag = 1; 与 while (flag == 0); print x;）在弱一致性机器上会因编译器/硬件重排序而失败，必须使用原子操作 + 正确的内存序（release/acquire）。同理，Cilk 的 cilk_sync、OpenMP 的 barrier/flush 语义之外，不要假设”时间上后写就一定后可见”；reducer 看似”无锁共享变量”，但它依赖运行时的同步语义，把普通共享变量当 reducer 用是数据竞争。

7. 思考题（带答案）

问题 1（MPI / 实现层）：某程序用 128 个 rank 做 1D 行块分解的 5 点 stencil，每 rank 负责 64 行、行宽 8192（float）。它用同步阻塞 Send/Recv 交换 halo，但运行 4 个 rank 的小规模测试时完全正常，扩到 128 个 rank 时偶发挂死。请解释原因、给出修复方案，并定量说明修复后每迭代的通信时间构成（假设 α = 2 µs、单 rank 有效带宽 6 GB/s）。

【答案】 挂死的根因是同步阻塞通信的循环等待：每个 rank 都先 Send 给上邻居、再 Send 给下邻居，然后才 Recv；同步 Send 只有在接收方确认数据已到达其地址空间后才返回。于是形成”所有人都已发出、所有人都在等对方接收”的环状依赖，构成死锁。小规模测试”通过”往往只是因为消息足够小、库走 eager 路径（发送方直接复制到网络缓冲立即返回），掩盖了问题；rank 数或消息大小变化后切换到 rendezvous 路径（需要接收方参与）才暴露死锁——这也解释了为什么”小规模正常、放大就挂”。修复方案（满足其一即可）：(a) 奇偶错开：偶数 rank 先 Send 后 Recv，奇数 rank 先 Recv 后 Send（讲义 06 中给出的修正版，T0→T1→T2→… 的流向因此在每个方向上都有接收方就绪）；(b) 用 MPI_Sendrecv（库内部保证不会自锁）；(c) 改用非阻塞 MPI_Isend/Irecv + MPI_Waitall（最推荐：既避免死锁，又能把通信与计算重叠）。定量上：halo 数据 = 2 行 × 8192 × 4 B = 65 536 B = 64 KB,传输时间 = 64 KB / 6 GB/s ≈ 10.9 µs,两条消息的启动延迟 = 2 × 2 µs = 4 µs（非阻塞时两条消息的延迟可部分重叠，紧下界约 2 µs）,所以每迭代通信 ≈ 11–15 µs,而同迭代计算时间约为 8192 × 64 × 20 B / 12.5 GB/s ≈ 839 µs（见 4.2 节），通信占比约 1.3%–1.8%。结论：修复后通信不是瓶颈；但如果把问题固定而把 P 增到 1024（rows = 8），通信字节数不变而计算量下降 8 倍，通信占比升到 ~10%–14%，这就是强扩展的真实上限来源（通信/计算比 ≈ 0.4/rows）。

问题 2（OpenMP / 调度实现）：一个循环有 N = 10⁹ 次迭代，每次迭代约 1 ns，动态领号临界区一次约 20 ns，机器有 64 个核。请推导最优 chunk 大小、给出该点的加速比，并说明为什么 schedule(guided) 往往能自动接近这个效果。

【答案】 动态调度的总时间可写成三项：并行计算项 T1/P、串行领号项 (N/c)·t_sync（领号必须串行，按 Amdahl 直接相加到关键路径）、不均衡项 c·t_iter（最后一个 chunk 只能由一个线程做，其间其他线程空闲）。即 T_P(c) ≈ T1/P + (N/c)·t_sync + c·t_iter。 对 c 求导：−N·t_sync/c² + t_iter = 0 ⇒ c* = sqrt(N·t_sync/t_iter) = sqrt(10⁹ × 20 ns / 1 ns) = sqrt(2×10¹⁰) ≈ 1.41×10⁵。此时两项开销相等：(10⁹/1.41×10⁵) × 20 ns = 1.41×10⁵ ns = 141.4 µs，c·t_iter = 141.4 µs，合计 283 µs。于是 T_P ≈ 1.0 s/64 + 283 µs = 15.625 ms + 0.283 ms ≈ 15.91 ms， 加速比 S = 1.0/0.01591 ≈ 62.9×，效率 62.9/64 ≈ 98.3%。对照 c = 1：串行领号 10⁹ × 20 ns = 20 s，T_P ≈ 20 s，加速比 0.05×（比串行慢 20 倍）；c = 1024：T_P ≈ 15.6 + 19.5 + 0.001 ≈ 35.1 ms，加速比仅 28.5×。所以 chunk 的选择直接决定了可达到的并行度上限：并行度 = T1/((N/c)t_sync + c·t_iter) = 1/(2×sqrt(t_sync·t_iter/N)) = sqrt(N/(4 t_sync t_iter)) ≈ sqrt(10⁹/(4×20)) ≈ 3536，远大于 64，因此在大 chunk 区间内 P 是限制项，在小 chunk 区间内串行项是限制项。schedule(guided) 的策略是”chunk 从大约 剩余迭代数/线程数 开始、随剩余量递减、并设一个下限”，前期大 chunk 控制领号次数（把 (N/c)t_sync 压到 O(P·t_sync·log) 量级），后期小 chunk 消除尾部不均衡（把 c·t_iter 压到与单次迭代同量级），因此无需精调就接近 c* 给出的折中点。补充说明：若迭代代价完全均匀且可预测，schedule(static) 的 T_P = T1/P（效率 100%）比上述任何动态方案都好，动态调度只在代价不均时才值得。

问题 3（Cilk / work stealing 实现）：讲义强调 Cilk 采用 “run child first”（continuation stealing，续体窃取），并指出对 for (i=0;i<N;i++) cilk_spawn foo(i); 这种写法，continuation stealing 相对 child stealing 有明显优势。请说明 (1) 两种策略在执行顺序与空间占用上的差别；(2) 为什么窃取要从 dequeue 的头部而不是尾部；(3) 若把 N 个 cilk_spawn 报告为”并行度 = N”，这个说法错在哪里——请用快排的例子定量说明。

【答案】 (1) 执行顺序与空间：continuation stealing 下，遇到 cilk_spawn foo(i) 时，当前线程立即执行 child foo(i)，并把”剩余工作的续体”（从 i+1 继续的循环控制）压入本地 dequeue，随后弹出（后进先出）。若无人窃取，执行顺序与删掉 cilk_spawn 的串行程序完全一致（深度优先遍历调用图），而且调用者线程一次只创建 1 个可被窃取的条目（代表剩余全部迭代），所以 dequeue 空间是 O(1) 量级；可以证明 work queue 的存储量不超过单线程栈存储量的 T 倍（T = 线程数）。child stealing 则相反：先把 child 记录起来、先把所有迭代的 child 都创建出来再执行，相当于对调用图做广度优先遍历，需要 O(N) 空间，且无窃取时的执行顺序与串行程序显著不同（局部性更差）。 (2) 为什么要从头部窃取：本地线程在尾部做 push/pop（后进先出 ⇒ 深度优先 ⇒ 缓存与栈局部性好），窃取者在头部取（先进先出 ⇒ 拿到的是调用树中”更早搁置、覆盖更多后续工作”的大块任务）。三个收益：(a) 头尾分离 ⇒ 两端落在不同 cache line ⇒ 与本地线程的竞争最小（这也是实现上必须把 head/tail 用 64 字节隔开的原因）；(b) 偷到的任务更大 ⇒ 一次窃取（0.1–5 µs）的成本被更长的后续计算摊薄，避免”偷得比偷的成本还少”；(c) 与 run-child-first 配合，本地线程始终留在自己调用树的局部，窃取者拿走远端子树，最大化双方的数据局部性。 (3) “并行度 = 任务数”是错的：并行度定义为 T1/T∞，即 Work 比关键路径（span），与 spawn 数量无关，只与最长依赖链有关。以并行快排（示例 3.3）为例：cilk_spawn 的次数是 Θ(n)，看起来”并行度 = 10⁸”，但 partition_lomuto 本身是串行的，且代码是”spawn 左半、本地继续右半”，关键路径 = 逐层 partition + 最后一次串行排序 = T∞ = Θ(n)。于是 T1 = Θ(n log n)、T∞ = Θ(n)，并行度 = Θ(log n)；在 n = 10⁸ 时 log₂n ≈ 27 ⇒ 并行度只有约 27，greedy 界 T_P ≤ T1/P + T∞ 在 P = 64 时给出 T_64 ≤ 0.042 + 0.1 ≈ 0.142 s（T1 = 2.7 s），加速比上限约 19×、效率 30%。正确的改法是缩短 span：把 partition 也并行化（span 降到 O(log²n)）并对两半都 spawn，此时 T∞ 只有微秒级，并行度升到 ~3.7×10⁶，P = 64 时加速比可接近 64×。这正说明：任务多 ≠ 扩展性好；决定扩展性的是关键路径，多 spawn 只增加了调度开销。