Lecture 7: Parallel Programming Basics

目录 · ← l6 · l8 →

Lecture 7: Parallel Programming Basics

1. 章节标题与概述

Lecture 7: Parallel Programming Basics（并行编程基础）

• 本讲核心问题：拿到一个串行程序（或一个问题描述），怎样系统地把它变成一个正确的、可扩展的并行程序？讲义给出的答案是四个必须依次回答的问题——Decomposition（分解：把问题切成可并行执行的 task）→ Assignment（分配：把 task 交给 worker）→ Orchestration（编排：通信、同步、数据布局、调度）→ Mapping（映射：把 worker 落到物理执行单元上）——并反复强调这条链条上真正的”拦路虎”是依赖（dependencies）：依赖决定并行度上限（Amdahl 定律），依赖决定你需要哪种同步机制（锁 / 屏障 / 消息 / 标志位），依赖决定你把数据切成什么形状（blocked 还是 interleaved，一维还是二维分块）才能既负载均衡又少通信。

• 涉及的主要硬件/软件机制：
  • 软件侧：三种编程模型的协同使用（共享地址空间 shared address space、消息传递 message passing、数据并行 data parallel）、SPMD（Single Program Multiple Data，单程序多数据）执行模型、for_all / ISPC foreach 数据并行循环、ISPC programCount/programIndex 与 launch[N] 任务、pthreads 的静态块分配、Cilk 式工作-跨度（work-span）模型、锁 lock / 屏障 barrier / 原子操作（atomic read-modify-write）/ 标志位 + 自旋（flag-based spin，等价于长度为 1 的消息队列）。
  • 硬件侧：多核 CPU 的缓存层次与 SIMD lane（ISPC 的 program instance 被编译器映射到向量 lane）、GPU 的 SM（CUDA thread block 由硬件映射到 SM）、集群节点的私有地址空间 + NIC + 网络（延迟 α、带宽 β）、以及决定一切的内存带宽与同步代价。
  • 抽象与实现的接口：同一份并行算法在三种模型下写出三种结构完全不同的代码（for_all / SPMD + 锁与屏障 / send-recv），而”哪种更好”取决于机器——讲义明确说 “it depends on the system this program is running on”。

• 在并行计算知识体系中的角色：这是全课程从”硬件/模型导论”（第 1–6 讲）转向”如何动手写并行程序”的分水岭，也是后续 Performance Optimization、同步与一致性、锁与无锁、异构计算等讲座的词汇表来源。第 8 讲起的性能优化（测量、局部性、伪共享）、第 10–17 讲的缓存一致性 / 内存一致性 / 同步 / 无锁数据结构，全部建立在本讲引入的 Decomposition/Assignment/Orchestration 框架与”依赖 → 同步”这条因果链上；作业 2（Assignment 2）在此后发布，正是要求把这些概念落到真实机器上。

• 配套材料：
  • lectures/06_progbasics.pdf（对应抽取文本 extracted/06_progbasics.txt，共 63 页幻灯片，首页明确写作 “Lecture 7: Parallel Programming Basics”，页脚学期字样为 Fall 2026）：已公开，位于 Fall 2026 公开讲义目录 https://www.cs.cmu.edu/~418/lectures/ 之下，可公开下载。Fall 2026 日程表 https://www.cs.cmu.edu/~418/schedule.html 中本讲位于 Sep 9。
  • 讲义第 2 页的课程通知（原文）：Assignment 1 due tonight（+ late days）、Assignment 2 released tonight、周五有一次 ~5 分钟的 GPUs/CUDA 小测验、所有人已加入 ed / Autolab / Gradescope。
  • 历史学期的消息传递专题讲义（如 f22_26-msgpassing、s23_20a_msgpassing 对应的 PDF）在本地仅有登录页抽取结果（内容为空），位于 /afs/cs/academic/class/15418-*/public/ 之下，需要 CMU 登录，属未公开。
  • 录制视频（Panopto/YouTube）：Fall 2026 日程表中被注释隐藏，属 未发布。Ed 讨论区、Autolab、Canvas、Gradescope 均需登录，非公开。
  • cs149_supp/thoughtprocess.txt（Stanford CS149 Fall 2025 Lecture 4 “Parallelizing Code: The Programming Thought Process”，74 页）：已公开的姊妹课程补充读物，内容与本讲高度重叠（同一套来自 Culler/Singh/Gupta 的网格求解器案例、Amdahl 定律、ISPC/pthreads 分配示例），本笔记在用到其中的额外数据（如 Summit 的 27,648 GPU × 5,376 ALU 规模算例）时会显式标注来源。
  • Fall 2026 授课教师为 Brian Railing 与 Dimitrios Skarlatos；课程由 Kayvon Fatahalian 创建。讲义首页的学期页脚沿用历史版本属正常现象。

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

2.1 起点：三种编程模型的复习，与”混合模型”的工程现实

• 定义与目的：并行编程模型是”程序员描述并行性、通信与同步的约定”。讲义复习了三种：
  1. 共享地址空间（shared address space）：通信隐式地发生在 load/store 里，非结构化。最自然，但”很容易打中自己的脚”——程序可能正确却毫无性能（因为隐式的通信/一致性代价看不见）。
  2. 消息传递（message passing）：所有通信都必须结构化地表达为消息（send/recv）。写出第一个正确版本的难度高于共享地址空间，但结构本身常常帮助你先写出正确且可扩展的程序（因为你被迫显式面对通信与所有权）。
  3. 数据并行（data parallel）：把计算组织成对一个集合的大 map；它假设存在共享地址空间来 load 输入 / store 结果，但严格限制 map 各次迭代之间的通信（目标是保持迭代的独立处理）。现代实现（CUDA、OpenCL 的 kernel 内同步原语）鼓励但不强制这种结构，允许有限的迭代间通信。
• 直观解释（”它是什么？”）：把三种模型想成同一个办公室里的三种协作方式。共享地址空间 = 大家共用一块白板，谁想写就写、谁想看就看——效率极高，但两个人同时改同一行数字就会互相覆盖。消息传递 = 每人有自己带锁的笔记本，要共享信息必须写便条塞进对方信箱——麻烦，但绝不会有人偷改你的笔记。数据并行 = 一份流水线作业单：把一叠文件平均发给每个人，每人按同一条规则各改各的文件，互不通信。
• 工程现实（现代实践：混合模型）：在集群的一个多核节点内部用共享地址空间编程，在节点之间用消息传递——这是”非常非常常见”的做法：在能高效实现共享地址空间的地方（节点内）享受便利，在别处要求显式通信。同时，数据并行风格的模型允许 kernel 内使用共享内存风格的同步原语（CUDA 的 __syncthreads()、OpenCL 的 barrier）。

2.2 从问题到机器：四个阶段（Decomposition / Assignment / Orchestration / Mapping）

• 定义与目的：讲义把”创建一个并行程序”的思维过程形式化为一条流水线。它是本讲的骨架，也是后续所有性能问题定位的坐标系：性能不好时，先问”是分解出了错（依赖没找对 / 任务太少），还是分配不均，还是编排（同步/通信）太贵，还是映射破坏了局部性”。
• 直观解释（”它是什么？”）：类比盖楼。Decomposition = 把整栋楼拆成一张张图纸（任务）；Assignment = 决定哪个施工队做哪张图纸；Orchestration = 安排塔吊、材料堆放、以及”哪道工序必须等哪道工序”的交接；Mapping = 把人真正派到工地上（甚至把相关工序安排在同一层楼以减少往返）。
• 架构/机制图解：

                    ┌──────────────────────────────────────────────┐
                    │   问题：解 PDE / 处理 N×N 图像 / N 体引力      │
                    └───────────────────────┬──────────────────────┘
                                            │
              (1) Decomposition  分解        ▼
   ┌──────────────────────────────────────────────────────────────────┐
   │ 目标：任务数 ≫ 执行单元数（"create at least enough tasks to      │
   │       keep all execution units busy"）                            │
   │ 核心动作：识别依赖（dependency）——依赖 = 并行度的上限            │
   │ 注意：分解不必是静态的，可在执行中动态产生新任务                   │
   │ 负责者：绝大多数情况是程序员；自动并行化编译器只在简单 loop nest  │
   │         上有 modest success，"魔法并行编译器"至今不存在           │
   └───────────────────────────────┬──────────────────────────────────┘
                                   │ (2) Assignment 分配
                                   ▼
   ┌──────────────────────────────────────────────────────────────────┐
   │ 把 task 绑到 worker（pthread / ISPC 实例 / warp / MPI 进程 / 向量lane）│
   │ 两条路：静态（blocked 连续块 / interleaved 交错）                  │
   │        动态（任务队列 + next task ptr / 工作窃取）                 │
   │ 双目标：balance workload（负载均衡）+ reduce communication cost    │
   │ 负责者：ISPC foreach → 系统；pthreads → 程序员；CUDA block → 硬件  │
   └───────────────────────────────┬──────────────────────────────────┘
                                   │ (3) Orchestration 编排
                                   ▼
   ┌──────────────────────────────────────────────────────────────────┐
   │ 内容：① 组织通信结构 ② 为保持依赖而加同步 ③ 组织内存中的数据结构   │
   │       ④ 调度任务                                                  │
   │ 目标：降低通信/同步成本、保持数据引用的局部性、降低开销            │
   │ 机器细节会反过来影响决策：同步贵 → 就用得更稀疏（coarser）         │
   └───────────────────────────────┬──────────────────────────────────┘
                                   │ (4) Mapping 映射
                                   ▼
   ┌──────────────────────────────────────────────────────────────────┐
   │ 逻辑 worker → 物理执行单元                                        │
   │ 例 1：OS 把 pthread 映射到某核上的硬件执行上下文                    │
   │ 例 2：编译器把 ISPC program instance 映射到向量指令的 lane          │
   │ 例 3：硬件把 CUDA thread block 映射到 GPU 的 SM                     │
   │ 决策：相关的线程放同一处理器（最大化局部性/数据共享，降低通信同步）   │
   │       不相关的线程放同一处理器（一个带宽受限 + 一个计算受限 → 互补） │
   └──────────────────────────────────────────────────────────────────┘

性能特征：四个阶段对性能的影响是乘性而非加性的。分解决定了理论上限（Amdahl），分配决定了负载不均衡造成的”最长工人”时间，编排决定了每次迭代的固定同步开销（屏障延迟 ~1–3 µs 量级、一次锁 100 ns 量级、一条网络消息 α ~1–2 µs + m/β），映射决定了缓存命中率与有效带宽（20 GB/s 是节点级共享资源，多线程争抢后每个线程分到的更少）。吞吐量 = min(算力上限, 带宽 × 算术强度) / (1 + 同步开销占比) 是这一讲的隐含公式。

2.3 Decomposition 与依赖：Amdahl 定律

• 定义与目的：Decomposition 把”要解决的问题”切成 subproblem（task）。其关键方面是识别依赖（或者依赖的缺失）。只要存在依赖，两个 task 就不能同时执行；因此依赖链的长度（span，见第 4 节）直接给出并行度的上界。Amdahl 定律把这件事量化：设 S 为总工作量中本质上串行的比例，则并行加速比上限为 1/S。
• 直观解释（”它是什么？”）：类比做年夜饭。你有 8 个灶眼（8 核），但”腌肉必须等 2 小时”这道工序无法被 8 个灶眼拆开。哪怕其他 100 道菜都能并行，总时间仍被那 2 小时钉死。S=0.05 意味着”不管你有多少核，最多快 20 倍”。
• 一个具体的两步算例（讲义原例）：对一张 N×N 图像做两步处理——第 1 步把所有像素亮度翻倍（元素间完全独立），第 2 步求所有像素的平均值（对所有元素的归约）。
  • 串行实现：两步各 ~N²，总时间 ~2N²，并行度 = 1。
  • 第一次并行尝试：第 1 步并行（N²/P），第 2 步仍串行（N²）→ 总时间 N²/P + N² → 加速比 ≤ 2（P→∞ 时也超不过 2）。这就是”小串行区域钉死大机器”的经典图示。
  • 第二次改进：第 2 步先并行求部分和（N²/P），再串行合并 P 个部分和（P）→ 总时间 2N²/P + P，当 N ≫ P 时加速比 → P。代价是”合并部分和”这一新增开销。

   时间                                       时间
   ▲  串行程序                                 ▲  并行程序（第 1 步并行 + 第 2 步串行）
   │  ┌──────────┐ 并行度 1                    │  ┌────┐ 并行度 P
 N²│  │ 步骤 2    │                            │  │步骤2│ 并行度 1（Amdahl 的毒药）
   │  ├──────────┤                            │  ├────┤
 N²│  │ 步骤 1    │                            │  │步骤1│
   │  └──────────┘                            │  └────┘
   └──────────────────────────────►            └──────────────────────────►
       总时间 2N²，加速比上限 = 2              总时间 N²/P + N²，S = 0.5 → ≤ 2

   改进后：步骤 2 也并行求部分和 + 串行合并
   ┌────┐ 并行度 P          ← N²/P
   ├────┤ 并行度 P          ← N²/P
   ├────┤ 串行合并          ← P        ← "overhead: combining the partial sums"
   └────┘

谁来负责 Decomposition？ 大多数情况是程序员。自动分解串行程序仍是困难的研究问题：编译器必须先做依赖分析，而很多依赖是数据相关（编译期不可知）的；研究者只在简单 loop nest 上取得有限成功。

2.4 Assignment：静态（blocked / interleaved）与动态（任务队列）

• 定义与目的：把 task 分配给 thread（”worker”）。目标是平衡负载与降低通信成本。可以静态做（编译期/启动时确定），也可以在执行中动态做。讲义强调：虽然分解通常由程序员负责，但很多语言/运行时接管了分配。
• 直观解释（”它是什么？”）：类比发传单。静态 blocked = 把街区切成 4 段，每人包一段；静态 interleaved = 每人隔三张发一张；动态 = 大家排在一个待发清单前，谁发完手里的就回来领下一批。街区每户门牌密度不同时（负载不均），只有动态方式能自动补齐。
• ISPC 中的两种 Assignment（讲义原例；ISPC = Intel SPMD Program Compiler，SPMD 方言，编译器把 program instance 映射到 SIMD lane）：

  • 程序员管理的静态交错分配：显式用 programCount 步长循环，把迭代交错地分给各 program instance（interleaved）。

    // 编译: ispc --target=avx2-i32x8 sinx.ispc -o sinx.o -O2
    export void sinx(uniform int N, uniform int terms,
                     uniform float* uniform x, uniform float* uniform result)
    {
        // 假设 N % programCount == 0
        for (uniform int i = 0; i < N; i += programCount) {
            int idx = i + programIndex;                 // 每个实例处理一个 lane
            float value = x[idx];
            float numer = x[idx] * x[idx] * x[idx];
            uniform int denom = 6;                      // 3!
            uniform int sign  = -1;
            for (uniform int j = 1; j <= terms; j++) {
                value += sign * numer / denom;
                numer *= x[idx] * x[idx];
                denom *= (2*j + 2) * (2*j + 3);
                sign  *= -1;
            }
            result[idx] = value;
        }
    }
    

  • 系统管理分配的 foreach：foreach 把独立的工作（迭代空间）交给系统，由系统决定”哪个迭代交给哪个 program instance”。抽象上允许动态分配，当前 ISPC 实现是静态的（交错）。

    export void sinx(uniform int N, uniform int terms,
                     uniform float* uniform x, uniform float* uniform result)
    {
        foreach (i = 0 ... N) {                        // 系统接管 assignment
            float value = x[i];
            float numer = x[i] * x[i] * x[i];
            uniform int denom = 6;
            uniform int sign  = -1;
            for (uniform int j = 1; j <= terms; j++) {
                value += sign * numer / denom;
                numer *= x[i] * x[i];
                denom *= (2*j + 2) * (2*j + 3);
                sign  *= -1;
            }
            result[i] = value;
        }
    }
    

• pthreads 的静态 blocked 分配（讲义原例）：启动一个线程处理数组前半，主线程处理后半，pthread_join 收尾。这是”程序员管理、静态、blocked”。
• ISPC 动态任务分配（讲义原例）：launch[100] my_ispc_task(...) 产生 100 个任务构成任务表，ISPC 运行时把任务发给 worker 线程池；分配策略：worker 完成当前任务后查看任务表，把下一个未完成的任务分给自己（next task ptr）——这正是”动态分配”的原型。

   launch[100] my_ispc_task(...) 产生的任务表（所有 worker 共享）
   ┌────┬────┬────┬────┬────┬────┬─────┬─────┐
   │ T0 │ T1 │ T2 │ T3 │ T4 │ T5 │ ... │ T99 │
   └────┴────┴────┴────┴────┴────┴─────┴─────┘
     ▲
     └── next task ptr（原子地取出并前移；谁空了谁来取）

   Worker 0 ──领 T0──► 算完 ──领 T4──► 算完 ──领 T8 ──► ...
   Worker 1 ──领 T1──► 算完 ──领 T5──► 算完 ──领 T9 ──► ...
   Worker 2 ──领 T2──►（慢，T2 遇到长尾）──────────► 算完 ──领 T10 ──►
   Worker 3 ──领 T3──► 算完 ──领 T6──► 算完 ──领 T7 ──► 算完 ──┐
                                                              │ 空闲时继续领
                                                              ▼
   性能特征：任务粒度↓ → 负载越均衡，但每次领取的原子操作 + 任务表争用（吞吐量）↑；
             任务粒度↑ → 争用↓，但尾部的负载不均↑（"最后一个大任务"决定结束时间）

2.5 Orchestration：锁、屏障与”用数据换依赖”

• 定义与目的：Orchestration 负责组织通信、加入同步以保持依赖、组织内存中的数据结构、调度任务。目标是降低通信/同步成本、保持数据访问的局部性、降低开销。讲义反复强调一句话：机器细节会影响这些决策——如果同步很贵，就应该用得更稀疏。
• 直观解释（”它是什么？”）：类比十字路口的交通。锁 = 单车道桥梁的红绿灯（一次只放一辆车）；屏障 = 整队集合点（所有人到齐才继续）；而”用数据换依赖” = 与其在集合点等所有人，不如给每个人一份自己的计数器，大家各记各的，最后再汇总。
• 案例：共享地址空间版网格求解器（SPMD，程序员负责同步）。依赖是”红格更新完才能更新黑格”，所以每个颜色相位之间必须有屏障：

   SPMD 版网格求解器的时间线（P = 4 个线程，SPMD = 所有线程执行同一个函数，
   threadId 不同 ⇒ 计算的网格区域不同）

   时间 ──────────────────────────────────────────────────────────────────────►
   T0  ▓▓▓ 红格 ▓▓▓│░░░ 黑格 ░░░│◆lock 归约◆│▒▒ 检查收敛 ▒▒│▓▓▓ 红格 ▓▓▓│ ...
   T1  ▓▓▓ 红格 ▓▓▓│░░░ 黑格 ░░░│◆lock 归约◆│▒▒ 检查收敛 ▒▒│▓▓▓ 红格 ▓▓▓│ ...
   T2  ▓▓▓ 红格 ▓▓▓│░░░ 黑格 ░░░│◆lock 归约◆│▒▒ 检查收敛 ▒▒│▓▓▓ 红格 ▓▓▓│ ...
   T3  ▓▓▓ 红格 ▓▓▓│░░░ 黑格 ░░░│◆lock 归约◆│▒▒ 检查收敛 ▒▒│▓▓▓ 红格 ▓▓▓│ ...
       └─ barrier ─┘└─ barrier ─┘           └── barrier ──┘
        (1) 相位边界        (2) 归约前必须            (3) 判定后必须广播 done
            红 → 黑            等所有线程写完 diff        否则有的线程读到旧的 done

   为什么是三个屏障？（讲义原问）
     屏障 1：红格读黑格 ⇒ 不同颜色相位之间必须全局同步
     屏障 2：diff 由所有线程用锁累加，收敛判定只由"读到完整 diff"的那个线程做出
     屏障 3：done 是共享变量，判定完成后必须让所有线程看到同一个值
   性能特征：每个 barrier 的代价是 O(1) 同步延迟（~1–3 µs 量级，随核数缓增），
             与"本轮计算量"无关 ⇒ 网格越小，屏障占比越高（见 §4.3 的数值算例）

• 互斥为什么必要（讲义的核心小例子）：diff += myDiff 在机器上展开为三条指令 r1 ← diff / r1 ← r1 + r2 / diff ← r1。若 T0 与 T1 交错执行：两者都读到 0，各自算出 1，各自写回 1 → 一次加法丢失（丢失更新）。因此这一组指令必须原子（atomic）。
• 保持原子性的三类机制：
  1. 临界区加锁/解锁：LOCK(myLock); /* critical section */ UNLOCK(myLock);
  2. 硬件支持的原子读-改-写指令（intrinsics，如 fetch_add、atomicAdd、atomicCAS）；
  3. 语言级原子块：atomic { ... }、atomicAdd(x, 10)。
• 性能陷阱与修复（讲义的原例）：把锁放在内层 (i,j) 循环里 → 每更新一个格子就加一次锁（还很可能造成 cache line 在核间反复弹跳）。修复：每个 worker 先在本地累加部分和 myDiff，每轮只加一次锁把部分和并入全局 diff。
• 用数据换依赖：把 3 个屏障减到 1 个（讲义 slide 50 的技巧）：把单一的 diff 换成 diff[3]，用 index = (index+1) % 3 轮转——相邻两次迭代使用不同的累加变量，从而消除”必须先把 diff 清零”的那条依赖。因为某线程清零 diff[(index+1)%3] 时，别的线程正在读的是 diff[index]，两者模 3 不同余，所以3 个副本是最小值（2 个不够：将会清掉别人正在读的缓冲区）。这就是”trade off footprint for removing dependencies“这一常见并行编程手法的实例。

   单个 diff：                           diff[3] 轮转：
   轮 k:  清零→算→归约→判→(轮 k+1 必须等清零完成)
          轮 k+1 的清零与轮 k 的判定冲突 ⇒ 需要额外屏障
                                       轮 k:  归约到 diff[k%3]；清零 diff[(k+1)%3]
                                              屏障；判 diff[k%3]；index++
                                       归约与"清零下一格"互不冲突 ⇒ 少一个屏障
                                       （代价：footprint ×3，且三者最好别放同一 cache line）

2.6 Barrier 是”保守”的：更细粒度的依赖表达

• 定义与目的：barrier(num_threads) 把所有线程切成相位（phase）：所有线程在屏障之前的计算全部完成后，任何线程在屏障之后的计算才能开始。它简单但保守——粒度太粗：本可以只等一个生产者，却等所有人。
• 直观解释（”它是什么？”）：类比电梯 vs 传呼机。屏障 = 所有人都要在一楼集合齐全才上楼（哪怕你只需要等一个人）；flag + 自旋 = 传呼机（对方一按你就收到，其余人照常工作）。
• 更细粒度的依赖（讲义原例）：两个线程，一个生产结果 x，一个消费它。T0 写 x = 1; flag = 1; 然后继续做与 x 无关的工作；T1 自旋 while (flag == 0); 然后 print x;。这实际上就是实现了一个长度为 1 的消息队列——说明消息传递与共享地址空间在表达能力上等价，差别只在”谁负责维护这个结构”。

• 同步原语对比表：

  原语	语义	典型延迟量级（8 核节点，示意）	优点	风险
  锁（uncontended lock/unlock）	互斥进入临界区	~20–100 ns	表达力强、易理解	争用时串行化；循环内加锁会毁掉性能
  原子 RMW（fetch_add/atomicAdd）	单次读-改-写不可分割	~10–50 ns（无争用）；争用同一地址时线性变差	无需临界区、可组合	同一地址的高度争用（全局计数器）成为瓶颈
  屏障 barrier	相位边界，全员到齐	~1–3 µs（随核数、NUMA 缓增）	简单、保守但正确	频繁小相位下占比爆炸（见 §4.3）
  标志位 + 自旋（flag）	一对一生产者-消费者	一次 cache line 传输（~100 ns）+ 轮询	细粒度、无全局同步	必须配内存序（acquire/release）否则编译器/CPU 重排可致错
  消息 send/recv	通过消息通信与同步	α ~1–2 µs + m/β	跨节点唯一选择；结构清晰	阻塞语义易死锁；小消息受延迟支配

2.7 Mapping：谁把 worker 放到硬件上？

• 定义与目的：Mapping 把逻辑 worker 映射到物理执行单元。讲师给的三条路径：OS（pthread → 某核上的硬件执行上下文）、编译器（ISPC program instance → 向量指令 lane）、硬件（CUDA thread block → GPU 的 SM）。
• 直观解释（”它是什么？”）：类比会议室安排。同一个项目的同事（相关线程）安排在同一间会议室（同一核/同一节点）→ 资料共享方便（缓存命中）；不相干的人安排在同一间（一个带宽受限、一个计算受限）→ 房间利用率最高。
• 架构/机制图解（节点级硬件结构——共享地址空间与消息传递在同一台机器上的分界线）：

   节点 0（多核 CPU + 私有/共享缓存 + 本地内存）              节点 1
 ┌──────────────────────────────────────────┐        ┌──────────────────────────┐
 │  Core0      Core1      Core2     Core3   │        │  Core0 ... Core3         │
 │  ┌────┐     ┌────┐     ┌────┐    ┌────┐  │        │  ┌────┐      ┌────┐      │
 │  │L1D │     │L1D │     │L1D │    │L1D │  │        │  │L1D │      │L1D │      │
 │  │32KB│     │32KB│     │32KB│    │32KB│  │        │  └──┬─┘      └──┬─┘      │
 │  └─┬──┘     └─┬──┘     └─┬──┘    └─┬──┘  │        │     │           │        │
 │  ┌─▼─────────▼─┐     ┌───▼────────▼──┐  │        │  ┌──▼───────────▼──┐     │
 │  │ L2 (256KB)  │     │ L2 (256KB)    │  │        │  │  L2 / LLC        │     │
 │  └──────┬──────┘     └───────┬───────┘  │        │  └────────┬─────────┘     │
 │         └────────┬───────────┘          │        │           │               │
 │            ┌─────▼─────┐                │        │     ┌─────▼─────┐         │
 │            │ LLC (共享) │                │        │     │ LLC (共享) │         │
 │            └─────┬─────┘                │        │     └─────┬─────┘         │
 │            ┌─────▼─────┐   NIC          │        │     ┌─────▼─────┐  NIC    │
 │            │ 本地 DRAM │───┴──┐         │        │     │ 本地 DRAM │──┴──┐   │
 │            └───────────┘      │         │        │     └───────────┘     │   │
 └───────────────────────────────┼─────────┘        └─────────────────────┼───┘
                                 └───────────  网络（α 延迟 + β 带宽）────────┘
     ▲ 共享地址空间在这里结束                    ▲ 消息传递从这里开始
     │ 隐藏的代价：一致性流量、伪共享、带宽争抢    │ 显式的代价：打包/拷贝、协议栈、延迟

   性能特征（数量级示意，用于第 4 节的模型）：
     L1 命中      ~4 周期   (~1.3 ns @3GHz)      片内带宽 >100 GB/s
     L2 命中      ~12 周期                       
     LLC/DRAM     ~60–100 ns                     节点内存带宽 ~20 GB/s（共享！）
     跨节点消息   α ~1–2 µs + m/β（β ~10 GB/s）  比 DRAM 访问慢 1–2 个数量级
   ⇒ 结论：跨节点只做"大块、少次"的通信（bulk transfer，整行整行地传），
           绝不逐元素发消息；节点内则要靠分块提升局部性、避免伪共享

2.8 案例研究：二维网格求解器（Gauss-Seidel → 对角波前 → 红黑重排）

• 问题定义：在 (N+2)×(N+2) 的网格上解偏微分方程（PDE），用迭代法：不断做 Gauss-Seidel sweep 直到收敛。更新式为 A[i,j] = 0.2 * (A[i,j] + A[i,j-1] + A[i-1,j] + A[i,j+1] + A[i+1,j])。
• 第一步：找依赖。原位（in-place）Gauss-Seidel 中，每个格子读左邻与上邻的新值，写自己的值（再被右邻与下邻读）。所以：
  • 沿一行有链式依赖（左 → 右）；沿列也有依赖（上 → 下）；
  • 但同一条反对角线 k = i+j 上的格子互不依赖 → 存在对角波前（wavefront）并行；
  • 坏消息：计算刚开始和快结束时并行度很小（波前长度从 1 涨到 N 再降回 1），而且每完成一条对角线就要同步一次（一次 sweep 需要 2N-1 次同步）。
• 第二步：为了并行而改变算法——红黑着色（red-black coloring）。改变网格单元的更新顺序：先并行更新所有红格，再并行更新所有黑格（黑格依赖红格的新值），重复直到收敛。注意这是允许的，因为它需要领域知识：新算法以不同的方式收敛到同一个解；浮点计算的中间值不同，但最终仍收敛到误差阈值之内。
• 架构/机制图解：

 (a) 原位 Gauss-Seidel 的依赖（一次 sweep 之内）
              ┌───────────────┐
              │  (i-1, j)     │   "每个格子依赖左邻与上邻的新值"
              └──────┬────────┘
                     ▼
   (i, j-1) ──────► (i, j) ──────► (i, j+1)
                     ▲
              ┌──────┴────────┐
              │  (i+1, j)     │
              └───────────────┘

 (b) 只有反对角线独立 ⇒ 波前并行（波前编号 k = i+j，同一 k 的格子可并行）
   ┌────┬────┬────┬────┬────┐     波前总数 = 2N-1
   │ 2  │ 3  │ 4  │ 5  │ 6  │     每推进一条波前同步一次  ⇒ Span(sweep) = Θ(N)
   ├────┼────┼────┼────┼────┤     波前长度 ~ min(k, N) ⇒ 并行度 ≈ N/2
   │ 3  │ 4  │ 5  │ 6  │ 7  │     两端"三角形"并行度低、同步频繁
   ├────┼────┼────┼────┼────┤
   │ 4  │ 5  │ 6  │ 7  │ 8  │
   └────┴────┴────┴────┴────┘

 (c) 红黑重排：把依赖"降维"成两个无依赖的相位
   ┌────┬────┬────┬────┬────┐     红格 (i+j 偶) 只读黑格 ⇒ 全体红格可同时更新
   │ R  │ B  │ R  │ B  │ R  │     黑格 (i+j 奇) 只读红格 ⇒ 全体黑格可同时更新
   ├────┼────┼────┼────┼────┤
   │ B  │ R  │ B  │ R  │ B  │     依赖： 红相位 ──barrier──► 黑相位 ──barrier──► 红相位 ...
   ├────┼────┼────┼────┼────┤             ▲                              │
   │ R  │ B  │ R  │ B  │ R  │             └──────────────────────────────┘
   └────┴────┴────┴────┴────┘
   Span(一个相位) = Θ(1)（相位内零依赖）；一次 sweep 只需 2 次相位同步
   并行度 = Θ(N²) —— 用"多迭代几次"换"每个相位内近乎无限的并行度"

 (d) 数据并行表达（交给系统编排）
   for_all (red cells (i,j)) {                    // 分解：单个网格元素 = 独立工作
       prev = A[i,j];                             // 分配：??（系统决定）
       A[i,j] = 0.2f*(A[i-1,j]+A[i,j-1]+A[i,j]+A[i+1,j]+A[i,j+1]);
       reduceAdd(diff, abs(A[i,j]-prev));         // 编排：系统提供的内置通信原语
   }                                              // for_all 末尾 = 隐式屏障（等所有 worker）

2.9 消息传递：私有地址空间、Ghost Cell、阻塞与非阻塞

• 定义与目的：消息传递模型下每个线程有自己的地址空间——没有共享变量，线程通过发送/接收消息来通信与同步。它把”数据所有权”变成显式概念：网格数据被切成 4 份，分别住在 4 个地址空间里（4 个私有数组）。
• 直观解释（”它是什么？”）：类比四个各自记账的仓库。仓库 2 要算自己的边界，就必须知道邻居仓库 1 和 3 的最后一排货物——于是邻居把那一排抄一份送过来贴在边界上，这一份复制就叫 ghost cell（影子/虚拟单元），其”所有权”仍在邻居手里。
• 数据复制是正确性的必要条件：因为不再有共享内存，边界数据必须在本地留一份拷贝（ghost cells）。每轮红格处理完后，线程 1 与线程 3 各发一行数据给线程 2（线程 2 下一相位需要最新的红格信息）。批量传输是准则：一次发整行，而不是逐元素发。
• 架构/机制图解（rank 的本地数组与 ghost row）：

        rank r-1                rank r                 rank r+1
   ┌──────────────┐       ┌──────────────┐        ┌──────────────┐
   │ row rows     │──────►│ row 0 (ghost)│        │              │
   │  (属于 r-1)   │  send │ row 1        │        │              │
   │   ...        │       │   ...        │        │              │
   │ row 1        │       │ row rows     │───────►│ row 0 (ghost)│
   └──────────────┘       │ row rows+1   │◄───────│ row 1        │
                          │  (ghost)     │  send  └──────────────┘
                          └──────────────┘
   代码形态：
     float* localA = allocate(rows_per_thread + 2, N + 2);
     recv(&localA[0,0],            sizeof(float)*(N+2), tid-1, MSG_ID_ROW);  // 上 ghost
     recv(&localA[rows_per_thread+1,0], sizeof(float)*(N+2), tid+1, MSG_ID_ROW); // 下 ghost
   性能特征：每个 rank 每轮通信量 = 2 行 × (N+2) × 4 B，
             而计算量 = (N/P) × N × k FLOP ⇒ 表面/体积比 ∝ P/N
             ⇒ 强扩展（固定 N）时 P 增大到通信与计算相当时就撞墙；
               弱扩展（N ∝ P，每 rank 行数固定）能保持效率

• 阻塞（synchronous）send/recv 的语义与死锁：
  • send() 在发送方收到”数据已驻留在接收方地址空间”的确认后才返回；recv() 在数据被拷入接收方缓冲区并回送 ack 后返回。
  • 于是”所有线程同时 sendDown“会形成循环等待：每个线程都阻塞在自己的 send 上，没有人执行 recv → 死锁。
  • 修复（讲义原方案）：让相邻两条线程的 (send, recv) 顺序相反——偶数号线程先 send 后 recv，奇数号线程先 recv 后 send。这样每一对相邻线程中总有一方率先发出消息，环上的等待被打破。
  • 非阻塞（asynchronous）send/recv：send() 立即返回（返回句柄 h1），但调用线程在 checksend(h1) 确认发送完成前不得修改发送缓冲区；recv() 只是”登记接收意图”并立即返回（句柄 h2），用 checkrecv(h2) 查询是否真的收完。好处：通信与应用线程执行重叠（消息处理与线程执行并发进行），线程可以在等消息时做别的工作。
  • 消息传递库通常还在 send/recv 之上提供高层原语：reduce_add(0, &my_diff, sizeof(float))（把所有 my_diff 加到 0 号线程）、broadcast(0, &done, sizeof(bool), MSG_DONE)（0 号线程把 done 广播给所有人）。

 (a) 阻塞 send/recv + 全员同时 send：循环等待 → 死锁（hang）
   T0 ──send↓──►[等 T5 recv]        ┌────────────────────────────────────────┐
   T1 ──send↓──►[等 T0 recv]        │ send() 语义：数据进入对方地址空间才返回    │
   T2 ──send↓──►[等 T1 recv]        │ 每个线程都在等对方先 recv                │
   ...                              │ 等待环无人打破 → 程序永久挂起（有时因      │
   T5 ──send↓──►[等 T4 recv]        │ 消息小而被 library 缓冲，于是"偶尔能跑"）  │
                                    └────────────────────────────────────────┘

 (b) 奇偶序修正（仍用阻塞 send/recv）
   偶数号 T0:  sendDown ─► recvDown ─► sendUp ─► recvUp
   奇数号 T1:  recvUp   ─► sendUp   ─► recvDown ─► sendDown
        ▲ 每一对相邻线程 (偶, 奇) 中：偶数先发、奇数先收
          ⇒ 每对里必有一方的 send 能被对方的 recv 匹配 ⇒ 无环 ⇒ 无死锁

 (c) 非阻塞（asynchronous）版本：把"等待"从时间轴上抹掉
   T0  Isend/send(h1) ─┐
   T0  Irecv/recv(h2) ─┤ 立即返回
   T0  计算不依赖 ghost 的内部行 ◄──── 通信与计算重叠（重叠期可能数百 µs）
   T0  Wait/checkrecv(h2) ──► 计算需要 ghost 的边界行
       ▲ 前提：在 checksend(h1) 之前绝不能改写发送缓冲区

3. 代码示例与性能分析

3.1 示例一：pthreads 静态块分配 vs 动态任务队列（sinx 的 Taylor 展开）

// 编译: g++ -O3 -march=native -pthread -std=c++17 sinx_pthreads.cpp -o sinx_pthreads
// 运行: ./sinx_pthreads 67108864 5 8 4096      # N=2^26, terms=5, 8 线程, 动态块 4096
#include <pthread.h>
#include <atomic>
#include <chrono>
#include <cstdio>
#include <cstdlib>
#include <cmath>

static int              g_N;              // 元素个数
static int              g_terms;          // Taylor 项数（<=5，避免 int denom 溢出）
static const float*     g_x;
static float*           g_result;
static std::atomic<int> g_next(0);        // 动态分配用的"next task ptr"
static int              g_chunk = 4096;   // 动态分配的块粒度

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

// 串行内核：对 [begin,end) 的每个元素做 sin(x) 的 Taylor 展开
//   sin(x) = x - x^3/3! + x^5/5! - x^7/7! ...
static void sinx_range(int begin, int end, int terms, const float* x, float* result) {
    for (int i = begin; i < end; i++) {
        float value = x[i];                      // 第 1 项
        float numer = x[i] * x[i] * x[i];        // x^3
        int   denom = 6;                         // 3!
        int   sign  = -1;
        for (int j = 1; j <= terms; j++) {
            value += (float)sign * numer / (float)denom;
            numer *= x[i] * x[i];                // 升幂：x^3 -> x^5 -> x^7
            denom *= (2*j + 2) * (2*j + 3);      // 3! -> 5! -> 7! ...
            sign  *= -1;
        }
        result[i] = value;
    }
}

// ---------- 静态 blocked 分配：线程 t 负责连续的一块 ----------
struct BlockArgs { int tid; int P; };

static void* blocked_worker(void* arg) {
    BlockArgs* a = (BlockArgs*)arg;
    long long  N = g_N;
    int begin = (int)(N * a->tid       / a->P);
    int end   = (int)(N * (a->tid + 1) / a->P);
    sinx_range(begin, end, g_terms, g_x, g_result);
    return nullptr;
}

// ---------- 动态分配：谁空了谁从共享计数器领一块 ----------
static void* dynamic_worker(void*) {
    for (;;) {
        int begin = g_next.fetch_add(g_chunk, std::memory_order_relaxed);
        if (begin >= g_N) break;                 // 领完了
        int end = begin + g_chunk;
        if (end > g_N) end = g_N;
        sinx_range(begin, end, g_terms, g_x, g_result);
    }
    return nullptr;
}

int main(int argc, char** argv) {
    g_N     = (argc > 1) ? atoi(argv[1]) : (1 << 26);
    g_terms = (argc > 2) ? atoi(argv[2]) : 5;
    int P   = (argc > 3) ? atoi(argv[3]) : 8;
    if (argc > 4) g_chunk = atoi(argv[4]);
    if (P < 1 || P > 64) { fprintf(stderr, "P 需在 1..64 之间\n"); return 1; }

    float* x      = (float*)malloc(sizeof(float) * (size_t)g_N);
    float* result = (float*)malloc(sizeof(float) * (size_t)g_N);
    for (int i = 0; i < g_N; i++)
        x[i] = -0.75f + 1.5f * (float)(i % 1000) / 999.0f;   // |x| <= 0.75，Taylor 收敛
    g_x = x; g_result = result;

    // ---- 串行基线 ----
    double t0 = now_sec();
    sinx_range(0, g_N, g_terms, x, result);
    double t_seq = now_sec() - t0;

    pthread_t th[64];
    BlockArgs args[64];

    // ---- 并行 A：静态 blocked ----
    t0 = now_sec();
    for (int t = 1; t < P; t++) { args[t] = BlockArgs{t, P};
        pthread_create(&th[t], nullptr, blocked_worker, &args[t]); }
    args[0] = BlockArgs{0, P};
    blocked_worker(&args[0]);                    // 主线程也当 worker（避免 P+1 路并行）
    for (int t = 1; t < P; t++) pthread_join(th[t], nullptr);
    double t_blk = now_sec() - t0;

    // ---- 并行 B：动态（next task ptr + atomic fetch_add） ----
    g_next.store(0);
    t0 = now_sec();
    for (int t = 1; t < P; t++) pthread_create(&th[t], nullptr, dynamic_worker, nullptr);
    dynamic_worker(nullptr);
    for (int t = 1; t < P; t++) pthread_join(th[t], nullptr);
    double t_dyn = now_sec() - t0;

    // ---- 正确性检查 ----
    double maxerr = 0.0;
    for (int i = 0; i < g_N; i++)
        maxerr = fmax(maxerr, fabs((double)result[i] - sin((double)x[i])));

    const double bytes = 8.0 * g_N;      // 每元素读 4B + 写 4B
    const double flops = 23.0 * g_N;     // 每元素约 23 次浮点运算（见下文分析）
    printf("N=%d terms=%d P=%d chunk=%d\n", g_N, g_terms, P, g_chunk);
    printf("serial  : %8.2f ms   %7.2f GFLOP/s  %7.2f GB/s\n",
           t_seq*1e3, flops/t_seq/1e9, bytes/t_seq/1e9);
    printf("blocked : %8.2f ms   %7.2f GFLOP/s  %7.2f GB/s   speedup = %5.2fx\n",
           t_blk*1e3, flops/t_blk/1e9, bytes/t_blk/1e9, t_seq/t_blk);
    printf("dynamic : %8.2f ms   %7.2f GFLOP/s  %7.2f GB/s   speedup = %5.2fx\n",
           t_dyn*1e3, flops/t_dyn/1e9, bytes/t_dyn/1e9, t_seq/t_dyn);
    printf("max |approx - sinf| = %.3e\n", maxerr);
    free(x); free(result);
    return 0;
}

• 【代码做什么？】
  1. 初始化 x[i] ∈ [-0.75, 0.75]（Taylor 展开在小 (|x|) 上收敛），result 作为输出。
  2. 串行跑一遍得到基线时间。
  3. 静态 blocked：把 [0, N) 按 N*t/P 切成 P 个连续块，线程 t 拿第 t 块（主线程也干活，避免只开 P-1 个真实 worker）。
  4. 动态：所有线程循环执行 g_next.fetch_add(g_chunk)，原子地领取下一块（块大小 4096 个元素），直到领到的起点越过数组末尾。这就是 ISPC launch[N] 任务表 + next task ptr 的手写版本。
  5. 与 sinf() 比较最大误差做正确性验证。
• 【并行机制与性能解说】
  • Work / Span / 并行度（本示例的核心分析）：
    • 记每个元素的内层递推（Taylor 迭代）代价为 c(T)（T = terms），则 Work（总工作量）W = N · (c₀ + c(T)) = Θ(N·T)。
    • Span（关键路径）：不同元素之间完全没有依赖（result[i] 只依赖 x[i]），所以并行性全部来自”N 个独立的串行链”；每个链的长度就是内层 for j 递推的长度 → Span = Θ(T)（与 N 无关）。
    • 并行度 = W / Span = Θ(N)。取 N = 2²⁶、T = 5：并行度 ≈ 6700 万，远远超过任何机器的线程数 → 并行度不是限制，内存带宽和除法吞吐量才是。
    • 对比：若把内层递推误写成跨元素依赖（例如累加器写在共享变量里），Span 立刻变成 Θ(N·T)，并行度掉到 1——这就是 Decomposition 阶段”识别依赖”的意义。
  • 硬件上如何并行执行：pthread_create 由 OS 把线程映射到不同 CPU 核（Mapping），每个线程执行 sinx_range 的各自区间；x[i] 与 result[i] 是私有工作集（不同线程写不同下标）→ 无共享写、无需同步、也无伪共享（除非 P 很大且块边界落在同一 cache line —— 4 字节元素、4096 元素/块，可以忽略）。
  • 瓶颈分析（Roofline 视角）：每元素搬 8 B（读 4 B + 写 4 B），做浮点运算的数量级为：内层每轮 ≈ 5 个浮点操作（numer 两次乘、x*x 一次乘、一次除、一次加）→ 5 轮 ≈ 25 次，加上收尾的几次，代码里按 23 FLOP/元素 估算 → 算术强度 AI ≈ 23/8 ≈ 2.9 FLOP/byte。而单核的 machine balance = 48 GFLOP/s ÷ 2.5 GB/s = 19.2 FLOP/byte → AI < balance ⇒ 每核也是带宽受限的：单核可期望 ≈ 2.5 GB/s × 2.9 ≈ 7.2 GFLOP/s，只有峰值的 15%。N = 2²⁶ 时总内存流量 = 8 B × 67.1 M = 537 MB，全节点带宽 20 GB/s 下界 ≈ 26.8 ms。
  • 额外杀手：浮点除法。numer / denom 在 x86 上是 ~5–10 倍于 FMA 的高延迟/低吞吐操作；每元素 5 次除法意味着 5 × 67.1 M = 3.35 亿次除法，仅除法吞吐就可能吃掉数百毫秒。优化方向：预计算 1.0f/denom 或使用 numer * inv_denom，把除法变成乘法；也可用快速多项式近似（这正是”性能优化”一讲的主题）。
  • 静态 blocked 与动态的差异：本内核每个元素成本完全均匀，因此两者应几乎相同（动态多出 fetch_add 的原子开销：2²⁶/4096 = 16384 次领取，可忽略）。把 g_chunk 调成 1 会看到动态版本变慢（原子争用 + 每块启动开销）；把元素成本改成不均匀（例如 terms 随 i 变化，模拟 Barnes-Hut 中”有的粒子周围邻居多、有的少”），静态 blocked 就会出现尾部不均，动态才扳回一局。
  • 可扩展性上限：由带宽给出，而非并行度。8 核全部投入时若每核分到 2.5 GB/s，则时间下界 26.8 ms；串行版本（单核带宽受限约 10 GB/s，且除法吞吐被单核独占）约 60–200 ms 量级 → 实测加速比会明显低于 8×，原因不是并行度不足，而是共享的内存带宽已饱和。

3.2 示例二：OpenMP 红黑网格求解器（屏障、部分和、收敛判定）

// 编译: g++ -O3 -march=native -fopenmp redblack_omp.cpp -o redblack_omp
// 运行: OMP_NUM_THREADS=8 ./redblack_omp 1024 200     # N=1024, 200 次 sweep
#include <omp.h>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cstring>

static int     N      = 1024;    // 内部网格 N x N（外圈一圈是固定边界）
static long    MAXIT  = 200;     // 最大 sweep 次数（默认跑满，便于测量稳态时间）
static double  TOL    = 1e-12;   // 收敛阈值：g_diff/(N*N) < TOL
static float*  A      = nullptr; // 工作网格 (N+2) x (N+2)
static float*  REF    = nullptr; // 串行参考网格
static double  g_diff = 0.0;     // 全局残差累加器（由各线程局部和归并）
static int     g_done = 0;       // 收敛标志（在 single 中写，在 barrier 后广播）
static int     g_iter = 0;       // 已完成的 sweep 数

#define IDX(i, j) ((i) * (N + 2) + (j))

// 串行参考：同一算法（红黑相位顺序）跑 iters 次，用于校验与测串行基线
static void run_ref(int iters) {
    memcpy(REF, A, sizeof(float) * (size_t)(N + 2) * (N + 2));
    for (int it = 0; it < iters; it++) {
        for (int color = 0; color < 2; color++) {
            for (int i = 1; i <= N; i++) {
                int jstart = ((i + color) & 1) ? 1 : 2;   // 使 (i+j)%2 == color
                for (int j = jstart; j <= N; j += 2) {
                    int p = IDX(i, j);
                    float prev = REF[p];
                    float v = 0.2f * (REF[p] + REF[IDX(i-1,j)] + REF[IDX(i+1,j)]
                                    + REF[IDX(i,j-1)] + REF[IDX(i,j+1)]);
                    REF[p] = v;
                    (void)prev;
                }
            }
        }
    }
}

int main(int argc, char** argv) {
    if (argc > 1) N     = atoi(argv[1]);
    if (argc > 2) MAXIT = atol(argv[2]);

    const size_t n2 = (size_t)(N + 2) * (N + 2);
    A   = (float*)aligned_alloc(64, sizeof(float) * n2);
    REF = (float*)aligned_alloc(64, sizeof(float) * n2);
    for (int i = 0; i < N + 2; i++)
        for (int j = 0; j < N + 2; j++)
            A[IDX(i,j)] = (i == 0 || j == 0 || i == N + 1 || j == N + 1) ? 1.0f : 0.0f;

    // ================= 并行版本：SPMD + 屏障 =================
    double t0 = omp_get_wtime();
    #pragma omp parallel
    {
        while (!g_done) {
            for (int color = 0; color < 2; color++) {          // 0 = 红, 1 = 黑
                float local = 0.0f;                            // 线程私有的部分和
                #pragma omp for schedule(static) nowait
                for (int i = 1; i <= N; i++) {
                    int jstart = ((i + color) & 1) ? 1 : 2;    // 只遍历该色格子
                    for (int j = jstart; j <= N; j += 2) {
                        int p = IDX(i, j);
                        float prev = A[p];
                        float v = 0.2f * (A[p] + A[IDX(i-1,j)] + A[IDX(i+1,j)]
                                        + A[IDX(i,j-1)] + A[IDX(i,j+1)]);
                        A[p] = v;
                        local += fabsf(v - prev);              // 局部累加，避免每格加锁
                    }
                }
                #pragma omp atomic
                g_diff += (double)local;                       // 每线程每相位只归并一次
                #pragma omp barrier                            // 屏障 1/2：相位边界（红<->黑）
            }
            #pragma omp barrier                                // 屏障 3：等所有局部和归并完毕
            #pragma omp single
            {                                                  // 单线程做收敛判定 + 清账
                g_iter++;
                if (g_diff / ((double)N * (double)N) < TOL || g_iter >= MAXIT)
                    g_done = 1;
                g_diff = 0.0;
            }                                                  // single 末尾隐式屏障：广播 g_done
        }
    }
    double t_par = omp_get_wtime() - t0;

    // ================= 串行基线（同算法、同迭代数） =================
    double t1 = omp_get_wtime();
    run_ref(g_iter);
    double t_ref = omp_get_wtime() - t1;

    // ================= 正确性校验：并行结果应逐位等于串行参考 =================
    float maxdiff = 0.0f;
    for (size_t k = 0; k < n2; k++) {
        float d = fabsf(A[k] - REF[k]);
        if (d > maxdiff) maxdiff = d;
    }

    const double cells = (double)N * (double)N;
    const double bytes = 24.0 * cells;        // 每格 5 次读 + 1 次写 = 24 B
    const double flops = 5.0  * cells;        // 4 次加法 + 1 次乘法
    printf("N=%d  sweeps=%d  threads=%d\n", N, g_iter, omp_get_max_threads());
    printf("parallel: %8.2f ms  (%6.3f ms/sweep)  %6.2f GB/s  %6.2f GFLOP/s\n",
           t_par*1e3, t_par*1e3/g_iter, bytes*g_iter/t_par/1e9, flops*g_iter/t_par/1e9);
    printf("serial  : %8.2f ms  (%6.3f ms/sweep)  speedup = %5.2fx\n",
           t_ref*1e3, t_ref*1e3/g_iter, t_ref/t_par);
    printf("max |parallel - serial| = %.3e  (红黑相位内部无依赖，期望为 0)\n", maxdiff);
    free(A); free(REF);
    return 0;
}

• 【代码做什么？】
  1. 分配 (N+2)×(N+2) 网格：外圈固定为 1.0，内部初值 0.0（这就是 Laplace 方程的边界条件，收敛解是内部全 1）。
  2. 并行区在 while 循环之外（#pragma omp parallel 只进入一次），线程在整个求解过程中常驻——避免每轮 fork/join 的巨大开销，这是”orchestration 要降低开销”的直接体现。
  3. 每轮分两个相位：color = 0 更新红格、color = 1 更新黑格。jstart = ((i+color)&1) ? 1 : 2 保证只遍历 (i+j)%2 == color 的格子，步长 2 沿行进。
  4. 每相位内每个线程本地累加残差 local，相位结束后用一次 #pragma omp atomic 并入 g_diff——而不是每更新一个格子就加一次锁（讲义 slide 46→47 的改进）。
  5. 三个屏障：相位边界、归约完成、done 广播（single 末尾自带隐式屏障，因此第三个屏障由 single 提供）。
  6. 收敛判定与 g_diff 清零放在 single 里（只有一个线程执行）；run_ref 用完全相同的更新表达式跑相同的 sweep 次数做逐位校验（红黑相位内无依赖，因此与串行结果应当完全一致）。
• 【并行机制与性能解说】
  • Work / Span / 并行度：
    • 一次 sweep 的 Work = N² 次单元更新（红黑各约 N²/2），记单格代价 c（含 5 次浮点与 24 B 访存）→ W_sweep = c·N²。
    • Span：红相位内部零依赖（红格只读黑格），span = c；黑相位同理 span = c；加上 3 次屏障，记屏障代价 B → Span_sweep = 2c + 3B。
    • 并行度 = W/Span = c·N²/(2c + 3B)。代入 N = 1024、c ≈ 24 B ÷ 20 GB/s = 1.2 ns、B ≈ 1.5 µs：并行度 = (1.05 M × 1.2 ns)/(2.4 ns + 4.5 µs) ≈ 281。这才是真正的可扩展性上限——不是 N²（≈ 10⁶），而是被屏障开销压到几百。要支持 1024 核，必须让 c·N² ≫ 3B，即网格规模要再大 3–4 个数量级或每格计算量大得多。
    • 对照讲义原版的 Gauss-Seidel 原位版本：一次 sweep 的 span = Θ(N)（沿行与沿列的两条长依赖链）→ 并行度只有 Θ(N/常数) ≈ 几百，且两端并行度低；红黑重排用”多做几次迭代”换来了 Θ(N²) 的相位内并行度，这是”改算法换并行度”的典型交易。
  • 硬件上如何执行：#pragma omp for schedule(static) 把行区间静态分给线程（blocked 划分，nowait 去掉隐式屏障，因为后面紧跟一个显式 barrier）；每个线程在自己的行块内按行扫描，行内访问连续 → 硬件预取器有效、SIMD 友好（gcc -O3 会把内层 j += 2 循环向量化）。相邻线程在行边界处共享边界行：这是真共享（只读），命中率高，不同于伪共享。
  • 瓶颈分析：
    • 内存带宽：每格 24 B，N = 1024 时一次 sweep = 25.2 MB，20 GB/s 下界 1.26 ms/sweep；200 sweeps ≈ 252 ms（计算只需 200 × 5.24 MFLOP = 1.05 GFLOP，8 核 384 GFLOP/s 下仅 2.7 ms）→ 该程序是纯带宽受限的（Roofline 见 §4.4）。
    • 屏障开销：3 × 1.5 µs = 4.5 µs/sweep，相对于 1.26 ms/sweep 仅 0.36%——N 大时同步可以忽略。但把 N 换成 128：一次 sweep 只有 393 KB → 19.7 µs，同步占总时间的 19%（相对纯计算时间则是 23%），且并行度降到 c·N²/(2c+3B) = 19.7 µs/4.5 µs ≈ 4.4 → 8 核都喂不饱。
    • 原子/锁：每线程每相位一次 atomic（≈ 2P 次/sweep），相对 1.26 ms 可忽略；但若把 atomic 放进内层循环（讲义 slide 46 的错误版本），就是 10⁶ 次/扫并伴随 cache line 争抢，性能会坍塌 1–2 个数量级。
  • 进一步优化：① 用 diff[3] 轮转把 3 个屏障降为 1 个（讲义 slide 50）；② 用 2D 分块 + 缓存分块（tiling）把 5 点模板的 24 B/格降到接近 8 B/格（cache 命中）；③ 用 #pragma omp for schedule(static,1) 或 dynamic 缓解不规则 workload；④ 把收敛检查降频（每 k 轮检查一次），进一步减少屏障。

3.3 示例三：MPI 非阻塞 Ghost 行交换的网格求解器

// 编译: mpicxx -O3 -march=native mpi_solver.cpp -o mpi_solver
// 运行: mpirun -np 4 ./mpi_solver 1024 200      # 内部网格 1024x1024, 200 次迭代
#include <mpi.h>
#include <cstdio>
#include <cstdlib>
#include <cmath>

int main(int argc, char** argv) {
    MPI_Init(&argc, &argv);
    int rank, nprocs;
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    MPI_Comm_size(MPI_COMM_WORLD, &nprocs);

    const int   n       = (argc > 1) ? atoi(argv[1]) : 1024;  // 内部网格 n x n
    const int   maxIter = (argc > 2) ? atoi(argv[2]) : 200;
    const float TOL     = 1e-9f;
    if (n % nprocs != 0) {
        if (rank == 0) fprintf(stderr, "要求 n 能被进程数整除\n");
        MPI_Finalize(); return 1;
    }

    const int rows = n / nprocs;   // 本 rank 拥有的行数（blocked 一维划分）
    const int W    = n + 2;        // 含左右边界的行宽
    const int lr   = rows + 2;     // 本地行数（含上下 ghost 行）

    float* A = (float*)malloc(sizeof(float) * (size_t)lr * W);
    for (int i = 0; i < lr; i++)
        for (int j = 0; j < W; j++) {
            int gi = rank * rows + i;      // 本地第 i 行对应的全局行号 0..n+1
            A[i*W + j] = (gi == 0 || gi == n + 1 || j == 0 || j == n + 1) ? 1.0f : 0.0f;
        }

    MPI_Barrier(MPI_COMM_WORLD);
    double t0 = MPI_Wtime();
    int   iter  = 0;
    float gdiff = 1e30f;

    while (iter < maxIter && gdiff / ((float)n * (float)n) > TOL) {
        // ---------- 1) 非阻塞交换 ghost 行（避免阻塞 send/recv 的死锁） ----------
        MPI_Request req[4];
        int nreq = 0;
        if (rank > 0) {                    // 与下方邻居 rank-1
            MPI_Isend(&A[1 * W],        W, MPI_FLOAT, rank-1, 0, MPI_COMM_WORLD, &req[nreq++]);
            MPI_Irecv(&A[0 * W],        W, MPI_FLOAT, rank-1, 1, MPI_COMM_WORLD, &req[nreq++]);
        }
        if (rank < nprocs - 1) {           // 与上方邻居 rank+1
            MPI_Isend(&A[rows * W],     W, MPI_FLOAT, rank+1, 1, MPI_COMM_WORLD, &req[nreq++]);
            MPI_Irecv(&A[(rows+1) * W], W, MPI_FLOAT, rank+1, 0, MPI_COMM_WORLD, &req[nreq++]);
        }
        MPI_Waitall(nreq, req, MPI_STATUSES_IGNORE);
        // 注：更激进的做法是先算"不需要 ghost 的内部行"，再 Waitall，最后算边界行，
        //     从而把通信完全藏在计算后面。

        // ---------- 2) 计算本 rank 拥有的 rows 行 ----------
        float mydiff = 0.0f;
        for (int i = 1; i <= rows; i++)
            for (int j = 1; j <= n; j++) {
                float prev = A[i*W + j];
                float v = 0.2f * (A[i*W + j] + A[(i-1)*W + j] + A[(i+1)*W + j]
                                + A[i*W + j-1] + A[i*W + j+1]);
                A[i*W + j] = v;
                mydiff += fabsf(v - prev);
            }

        // ---------- 3) 全局收敛判定：Allreduce = reduce + broadcast ----------
        MPI_Allreduce(&mydiff, &gdiff, 1, MPI_FLOAT, MPI_SUM, MPI_COMM_WORLD);
        iter++;
        if (rank == 0 && (iter % 50 == 0 || iter == 1))
            printf("  iter %4d  residual = %.4e\n", iter, gdiff / ((float)n * n));
    }
    double t1 = MPI_Wtime();

    if (rank == 0) {
        const double cells = (double)n * n * (double)iter;
        printf("n=%d  procs=%d  rows/rank=%d  iters=%d\n", n, nprocs, rows, iter);
        printf("time = %.2f ms  (%.3f ms/iteration)\n", (t1-t0)*1e3, (t1-t0)*1e3/iter);
        printf("aggregate: %.2f GB/s  %.2f GFLOP/s\n",
               24.0*cells/(t1-t0)/1e9, 5.0*cells/(t1-t0)/1e9);
    }
    free(A);
    MPI_Finalize();
    return 0;
}

• 【代码做什么？】
  1. 把 n×n 内部网格按行块（一维 blocked 划分）分给 nprocs 个进程；每个进程分配 (rows+2) × (n+2) 的私有数组，其中第 0 行与第 rows+1 行是要从邻居收到的 ghost 行。
  2. 每轮开始：用 MPI_Isend/MPI_Irecv 与上下邻居双向交换各一行（本进程的第 1 行送给下方、第 rows 行送给上方），MPI_Waitall 等齐。
  3. 对本地 rows 行做 5 点原地更新，累加本地残差 mydiff。
  4. 用 MPI_Allreduce 求全局残差并判定收敛（等价于讲义里的 reduce_add + broadcast 两个高层原语），然后进入下一轮。
  5. 在 rank 0 打印每次迭代时间、聚合带宽与算力。
• 【并行机制与性能解说】
  • Work / Span / 并行度：
    • 全局 Work = c·n² 每迭代（c 为单格代价）。
    • Span：每个 rank 内部按行递增的原地更新形成串行链：第 i 行依赖第 i-1 行的新值 → 链长 rows × n 步（行间还要等 ghost 到达，加上消息延迟 α）。因此 Span_iter ≈ c·(rows·n) + (消息延迟与 Allreduce 的链长)， 并行度 = W/Span ≈ c·n² / (c·n²/P + sync) → 趋近 P（而不是 n²）。
    • 这是非常重要的结论：“块内 Gauss-Seidel”的并行度等于进程数 P，即”用几个进程就最多只能用几个进程的并行度”，与红黑版本（并行度 Θ(n²)）形成鲜明对比。如果希望 MPI 版本也能在块内使用多线程并行（混合编程模型），必须先把块内的 GS 换成红黑或 Jacobi 更新。
  • 硬件上如何执行：mpirun 启动的进程由 OS 映射到不同节点/核（Mapping）；通信走 NIC + 网络（延迟 α、带宽 β）；数据所有权是显式的——ghost 行是拷贝，所有权在邻居手里，所以每轮都必须在关键路径之前把它们”刷新”一遍。这就是把讲义中的共享内存版（隐式通信、锁 + 屏障）替换成显式消息后的代价：通信从”看不见的 load”变成”必须排进关键路径的阶段”。
  • 瓶颈分析（α-β 模型）：单条消息时间 T = α + m/β。取 α = 2 µs、β = 10 GB/s，并以下面的强扩展场景为例（n = 4096、P = 64、rows = 64、每节点 8 rank）：
    • 每轮 4 条 ghost 行消息（2 发 2 收，非阻塞可重叠），每行 W = n+2 = 4098 个 float → m = 16.4 KB → 单条 2 + 16.4 KB/10 GB/s = 3.64 µs；双向重叠后关键路径 ≈ 5 µs 量级。
    • MPI_Allreduce 在 P = 64 时约 2·log₂P·α = 12 × 2 = 24 µs —— 它比 ghost 交换更贵，因为它的负载极小（一个 float）却要吃满整个归约树的延迟。教训：每轮都做全局归约是延迟灾难，应降频（每 k 轮判定一次）或用分层归约。
    • 计算侧：rows = 64、n = 4096 时，每 rank 每轮触及 64 × 4096 × 24 B = 6.29 MB；若一个节点跑 8 个 rank，每 rank 可用带宽 ≈ 20/8 = 2.5 GB/s → 2.52 ms/迭代。同步 29 µs 只占 1.2% → 该配置下通信不是问题。
    • 但问题变小就会立刻撞墙：n = 256、P = 64（rows = 4）时，每 rank 每轮只触及 4 × 256 × 24 B = 24.6 KB → 9.8 µs，而同步仍是 29 µs → 开销 75%，效率跌到 ~25%。这是”强扩展必须让问题足够大”的定量说明。
  • 表面/体积比：1D 行划分下通信量 ∝ 2 行 × n = O(n)，计算量 ∝ n²/P → 比值 ∝ P/n。固定 n 增大 P 必然撞墙；弱扩展（n ∝ √P，即每 rank 的行数固定）才能保持效率。若改用 2D 分块，通信 ∝ 周长 = O(n/√P)，比值改善 √P 倍——但代价是数据要打成方块（非连续行），需要派生数据类型或打包缓冲。

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

4.1 Amdahl 定律与强/弱扩展

• Amdahl（强扩展，固定问题规模）：设 S 为固有串行比例，

  \[\text{Speedup}(P) = \frac{1}{S + \frac{1-S}{P}}, \qquad \lim_{P\to\infty}\text{Speedup} = \frac{1}{S}\]

  数值表（同一 S 下不同 P 的上限）：

  S（串行比例）	P=1	P=16	P=64	P=256	P=1024	P→∞
  0.01	1.00	13.91	39.26	72.11	91.18	100
  0.05	1.00	9.14	15.42	18.62	19.64	20
  0.10	1.00	6.40	8.77	9.66	9.91	10

  参照系（来自公开的姊妹课程补充读物）：Summit 超算有 27,648 块 GPU × 5,376 ALU = 148,635,648 个 ALU。即使 S 小到 0.001（上限 1000 倍），能用上的执行单元也只占机器的 0.0007%——这就是”小串行区域钉死大机器”的量化形态。

• Gustafson（弱扩展，问题规模随机器变大）：若并行部分随 P 增长而串行部分固定，则 Speedup_scaled = P - S(P-1)。S = 0.05、P = 64 → 64 − 0.05×63 = 60.85。两定律不矛盾：Amdahl 说”固定规模下加机器收益递减”，Gustafson 说”机器变大时应该解更大的问题”。这正是 §3.3 中”弱扩展能保持效率”的理论依据。
• 与 Work-Span 的关系：Amdahl 只考虑”串行段”，忽略了并行部分的同步与通信开销——这在 P 大、问题小时是致命遗漏。Work-Span 模型与 α-β 模型补上了这一块。

4.2 Work-Span（Cilk）模型：并行度与贪婪调度定理

• 定义：
  • Work T₁：用 1 个处理器执行整个计算的时间（= 总工作量，所有串行化后的和）。
  • Span（关键路径）T∞：在无限多处理器上的执行时间，等于依赖图中最长路径的长度。
  • 并行度（parallelism）P∞ = T₁ / T∞：机器再大也超不过它。
  • 贪婪调度定理：任何贪婪（greedy）调度器满足 \(T_P \le \frac{T_1}{P} + T_\infty\) 且显然 T_P ≥ max(T₁/P, T∞)。因此：
    • 当 P ≪ T₁/T∞（parallel slackness 充足）时，T_P ≈ T₁/P，线性加速；
    • 当 P ≳ T₁/T∞ 时，加速比饱和在 T₁/T∞。

• 本讲各案例的 Work/Span 汇总表：

  计算	Work T₁	Span T∞	并行度 T₁/T∞	实际限制因素
  sinx（N 元素，T 项）	Θ(N·T)	Θ(T)	Θ(N)（≈6.7×10⁷）	内存带宽 + 除法吞吐（并行度早已过剩）
  两步图像处理（步 2 串行）	2N²	2N²	1	Amdahl：S = 0.5 ⇒ 加速比 ≤ 2
  两步图像处理（步 2 并行部分和 + 串行合并）	2N²	Θ(P)（串行合并的链长）	≈ 2N²/P	合并开销 P（P ≪ N 时可忽略）
  红黑网格求解器（每次 sweep）	c·N²	2c + 3B	c·N²/(2c+3B)（N=1024 时 ≈281）	屏障开销 B（不是浮点算力）
  Gauss-Seidel 原位（每次 sweep）	c·N²	Θ(c·N)	Θ(N)	对角波前依赖链 + 每波前一次同步
  MPI 块内 GS 求解器（每次迭代）	c·n²	c·n²/P + sync	≈P	通信/归约在关键路径上（延迟 α）

4.3 数值算例 A：两步图像处理的最优规模（Amdahl 的解析解）

• 假设：每处理一个像素记 1 个时间单位（忽略常数）。方案为”步 1 并行 + 步 2 并行求部分和 + 串行合并 P 个部分和”，则 \(T(P) = \frac{N^2}{P} + \frac{N^2}{P} + P = \frac{2N^2}{P} + P,\qquad T_1 = 2N^2\) \(\text{Speedup}(P) = \frac{2N^2}{2N^2/P + P} = \frac{2P}{2 + P^2/N^2},\qquad E(P)=\frac{2N^2}{2N^2+P^2}\)

• 代入数字：

  N	P	加速比	效率
  1024	64	2×64 ÷ (2 + 4096/1048576)=128/2.0039 = 63.88	99.8%
  64	64	128 ÷ (2 + 4096/4096) = 128/3 = 42.67	66.7%
  64	1024	2048 ÷ (2 + 1048576/4096) = 2048/258 = 7.94	0.78%
  1024	64（步 2 完全串行）	2 ÷ (1 + 1/64) = 1.97	3.1%

• 结论：要维持 95% 效率需要 E ≥ 0.95 ⇒ 2N² ≥ 0.95(2N²+P²) ⇒ N ≥ 3.082·P。P = 1024 时需要 N ≥ 3156（代入 N = 3155：E = 2×9,954,025/(2×9,954,025+1,048,576) = 0.94997，恰在 95% 之下；N = 3156：E = 0.95000）。也就是说：在 1024 路的机器上，图像至少要有约 3156² 个像素，串行合并的那 P 步才不碍事。若把串行合并换成树形归约（代价 log₂P 而非 P），则 T(P) = 2N²/P + log₂P，P = 1024 时 N = 1024 也能达到 E = 2N²/(2N²+P log P) = 2,097,152/(2,097,152+10240) = 99.5%。

4.4 数值算例 B：5 点模板的 Roofline（为什么网格求解器永远吃不满机器）

• 机器参数（假设一台 8 核节点）：

  参数	取值	推导
  核心数	8	—
  频率	3.0 GHz	—
  SIMD 宽度	8（AVX2，float）	—
  FMA 吞吐	2 FLOP/指令	乘加各算 1 次
  峰值算力	8 × 3.0 GHz × 8 × 2 = 384 GFLOP/s	—
  峰值内存带宽	20 GB/s（实测可达的节点级共享带宽）	DDR4-3200 理论 25.6 GB/s
  山脊点（ridge point）	384 / 20 = 19.2 FLOP/byte	峰值算力 ÷ 峰值带宽
  单核均摊带宽	20/8 = 2.5 GB/s	带宽是共享资源
  内存延迟	~80–100 ns	—
  屏障延迟	~1.5 µs（8 核）	—
  网络	α = 2 µs，β = 10 GB/s	—

• 内核的算术强度：5 点模板每格读 5 个 float（20 B）、写 1 个 float（4 B）= 24 B；浮点操作 4 加法 + 1 乘法 = 5 FLOP → AI = 5/24 = 0.208 FLOP/byte。

   性能 (GFLOP/s，对数轴)
   384 ┤●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●  峰值：8×3.0GHz×8 宽×2(FMA)
       │                                        ╱
       │                                      ╱   ← 带宽屋顶：可达性能 = 20 GB/s × AI
       │                                    ╱
       │                                  ╱
       │                                ╱
       │                              ╱
     4.2┤                            ●   ← 5 点模板：AI = 0.208 → 4.17 GFLOP/s（峰值的 1.1%）
       │                          ╱
       │                        ╱
       │                      ╱
       └────────────────────┴───────────────────────────────────────────►
                          19.2                                     AI (FLOP/byte，对数轴)
                        山脊点 = 峰值算力 / 峰值带宽
       ▲ AI = 0.208 远在山脊左侧 ⇒ 受限的是内存带宽，不是浮点算力
         提高 AI 的唯一途径：分块（tiling）让每个字节被复用更多次

• 数值结论：
  • 可达性能 = 20 GB/s × 0.208 = 4.17 GFLOP/s，只有峰值 384 GFLOP/s 的 1.09%。
  • N = 4096 时一次 sweep：16,777,216 格 × 24 B = 402.7 MB → 至少 20.1 ms；而纯计算只需 83.9 MFLOP ÷ 384 GFLOP/s = 0.22 ms → 内存受限约 92 倍。
  • 把网格缩到 N = 512（1.05 MB，装进 L2）：每格只需与内存交换 8 B（compulsory miss + 写回），AI 提到 0.625，一次 sweep 2.1 MB ÷ 20 GB/s = 105 µs，比 N=4096 的”每格 24 B”整整便宜 3 倍（相同总格数下）。这就是”orchestration 里的数据结构/分块组织”为什么值钱。
• Little’s Law（延迟隐藏所需的并发）：维持 20 GB/s 在 100 ns 延迟下需要 并发字节数 = 带宽 × 延迟 = 20 GB/s × 100 ns = 2000 B ≈ 31 条 64 B cache line。 单线程若同时只有 ~10 条未完成的 load（每条 64 B），则只能撑起 ~6.4 GB/s；因此要么多线程（≥4 个），要么让编译器向量化出 8 宽 load（一次 32 B，且 8 个元素一次取 2 条 line），要么靠硬件预取。这解释了为什么 §3.1 的 sinx 只有 8 线程时就已经摸到带宽墙：并行度不是问题，“在飞行中的字节数”才是。

4.5 数值算例 C：消息传递的 α-β 模型与同步摊销

• 模型：单条消息 T = α + m/β，取 α = 2 µs、β = 10 GB/s（跨节点）。
• 场景 1（大问题，n = 4096，P = 64，每节点 8 rank）：
  • ghost 行：每行 W = n+2 = 4098 个 float → m = 16.4 KB → 单条 = 2 + 1.64 = 3.64 µs；非阻塞双向重叠后关键路径 ≈ 5 µs。
  • Allreduce：2·log₂64 × α = 24 µs（负载只有 4 B，纯延迟开销）。
  • 每轮同步合计 ≈ 29 µs；每轮计算（带宽受限）≈ 2520 µs → 开销 1.2%。
• 场景 2（强扩展、问题变小，n = 256，P = 64，rows = 4）：
  • 每 rank 计算 = 4 × 256 × 24 B = 24.6 KB ÷ 2.5 GB/s = 9.8 µs；
  • 同步仍是 29 µs（消息变小但 α 不变，归约延迟不变）→ 开销 75%，效率 ≈ 25%。
• 结论与设计准则：
  1. 延迟项 α 不随问题缩小而缩小，所以”小问题 + 多进程”必然低效；强扩展的极限由 α / (每 rank 计算时间) 决定。
  2. 全局归约要降频（每 k 轮判定一次收敛），或用分层/树形归约，避免每轮付 2·log₂P·α。
  3. 通信要批量化：把 m 从 4 B 提到 16 KB，α 的占比从 100% 降到 55%；再大就完全被 β 支配（此时带宽成为瓶颈）。
  4. 表面/体积比：1D 划分下通信/计算 ∝ P/n；2D 分块改为 ∝ √P/n → 更抗强扩展。取 n = 4096、每格 24 B、网络 10 GB/s、内存 2.5 GB/s，令 ghost 流量折算后等于内存流量的 10%：可得 1D 划分大约在 P ≈ 4900 时才真正被通信拖住——这解释了为什么这类模板计算”很好扩展”，而它真正的墙依然是单节点的内存带宽（§4.4）。

5. 关键要点

1. 并行编程的思维过程是四个阶段，而且顺序不能颠倒：Decomposition（切出足够多、依赖尽量少的 task）→ Assignment（把 task 分给 worker，静态 blocked/interleaved 还是动态队列）→ Orchestration（通信、同步、数据布局、调度）→ Mapping（worker 落到的物理单元）。性能不佳时先定位卡在哪一阶段，而不是盲目调线程数。同一份算法在共享地址空间、消息传递、数据并行三种模型下会写出三种截然不同的代码，“哪个更好”取决于机器（”it depends on the system this program is running on”）。
2. 依赖决定一切：Amdahl 决定上限，Span 决定可扩展性天花板。Amdahl：串行比例 S 使加速比 ≤ 1/S（S = 0.05 时哪怕 148M 个 ALU 也只能用上 20 倍）。Work-Span：T_P ≤ T₁/P + T∞，可用并行度 = T₁/T∞。两者的实践含义是：先找依赖、再谈优化；能通过”改算法”（如 Gauss-Seidel 原位 → 红黑重排）把 Span 从 Θ(N) 降到 Θ(1) 时，代价往往是多做几次迭代——这是值得做的交易。
3. 同步是”用数据换依赖”的工程：屏障简单但保守（粒度粗），锁表达力强但贵，标志位+自旋最细粒度但必须配内存序。三条实用准则：① 把同步移出内层循环（本地部分和 → 每轮一次归约，而不是每格一次加锁）；② 减少同步次数（用 diff[3] 轮转消除”清零”依赖，把 3 个屏障降到 1 个）；③ 用空间换依赖（多副本、ghost cell、私有地址空间复制）是并行编程的常规武器。
4. 机器参数决定优化方向，而内存带宽通常是第一瓶颈。用 Roofline 判断受限类型：5 点模板 AI = 0.208 FLOP/byte，远低于 8 核节点 19.2 FLOP/byte 的山脊点 → 只有峰值的 1.1%，提高 AI 的唯一途径是分块复用。用 Little’s Law 判断是否够并发：带宽 × 延迟 = 20 GB/s × 100 ns = 2000 B ≈ 31 条 cache line 必须在飞行中。
5. 消息传递里”所有权”和”延迟”是显式成本：数据要复制（ghost cell），通信要排进关键路径，阻塞 send/recv 会构成循环等待而死锁（靠奇偶序或非阻塞原语打破）。α 项不随问题缩小而缩小，所以强扩展的极限是”每 rank 计算时间 ≫ α”；通信要大块、少次、能重叠（Isend/Irecv + 先算内部再算边界），全局归约要降频或分层。

6. 常见陷阱与注意事项

• 把 Amdahl 当成”加速比上界”的死结论，却忽略开销项与常数项：T(P) = T₁/P + T∞ + sync(P)——真正在小问题上钉死你的是第三项（屏障/消息的固定延迟），它不在 Amdahl 的 S 里。P = 64、n = 256 的 MPI 求解器同步占 75%，而按 Amdahl 看”串行比例很小、应该能扩展”——这就是漏掉同步项的典型误判。
• 在内层循环里加锁或做全局归约：讲义 slide 46 的版本对每一个网格单元执行一次 lock; diff += ...; unlock，而正确做法是先在寄存器里累加 myDiff，每轮只归并一次（slide 47）。前者不仅把临界区数量放大 10⁶ 倍，还会让存 diff 的 cache line 在核间反复弹跳（cache line ping-pong），性能下降一两个数量级。
• 伪共享（false sharing）：两个线程写同一 cache line 内的不同变量，硬件仍会按整条 line 传递所有权。把每线程的累加器放进一个数组（float partial[P]）时，相邻元素会落在同一 line；把 diff[3] 的 3 个副本放同一 line 也一样。修复：按 64 B 对齐/填充（padding）或每线程一个独立变量。注意 §3.2 中相邻线程共享网格边界行是真共享（只读），不是伪共享——两者要分清。
• 依赖没找全就动手并行：典型的错误是”把顺序执行的 Gauss-Seidel 直接 #pragma omp parallel for“——同一行内左邻的新值会被右邻读到，并行后会得到与串行不同的结果（通常是静默的错误答案）。要么用红黑/波前重排（需领域知识确认收敛性不受影响），要么改成 Jacobi（同时读旧值），要么加同步（那就退化成串行）。所谓”魔法并行化编译器至今不存在”说的就是这件事：依赖是数据相关的，编译期往往不可知。
• 误以为”消息传递里 send 返回就万事大吉”：阻塞 send() 的语义是”数据已到达对方地址空间”，全员同时 sendDown 会形成循环等待而死锁（小消息时库可能替你缓冲，于是表现为”有时能跑、有时挂住”——最坏的 bug 形态）。修复要么用奇偶序（偶发奇收），要么用非阻塞 Isend/Irecv。并且非阻塞版本有额外约束：在 checksend/Wait 确认之前绝不能改写发送缓冲区——这个约束被违反时，错误是间歇性的、极难复现的。
• 忽略内存带宽与”在飞行字节数”：以为”线程开得越多越快”，实际每核只能分到 20/8 = 2.5 GB/s；sinx 在 8 线程时已撞带宽墙，再加线程只会增加争抢。另一个常见误判是只数”FLOP”：sinx 每元素 5 次浮点除法（吞吐只有 FMA 的 1/5–1/10），改写为乘倒数能把这段代码从除法受限变成带宽受限——先测量再优化，不要凭直觉猜瓶颈。
• 同步不足与内存序错误：用自旋标志位做生产者-消费者（x = 1; flag = 1; / while (flag == 0); print x;）时，若 flag 不是原子且没有 acquire/release 语义，编译器和 CPU 都可以重排，导致消费者看到 flag == 1 却读到旧的 x。同理 §3.2 的并行区里，g_done 的可见性完全依赖 single 末尾的隐式屏障/flush——擅自删掉那个屏障，就可能出现”部分线程多算一轮、部分线程提前退出”的死锁或错误。

7. 思考题（带答案）

问题 1（Amdahl + 开销项的定量设计） 对 N×N 图像做两步处理：步 1 亮度翻倍（完全独立），步 2 求平均值（归约）。你采用”步 1 并行 + 步 2 并行求部分和 + 串行合并 P 个部分和（每次合并记 1 单位时间）”的方案，每个像素的步 1/步 2 各记 1 单位时间。 (a) 写出 T(P)、Speedup(P) 与效率 E(P)。 (b) 在 P = 1024 时，要让 E ≥ 95%，N 至少要多大？ (c) 如果把串行合并改成树形归约（代价 log₂P），在 P = 1024、N = 1024 时的效率是多少？

【答案】 (a) 步 1 时间 N²/P，步 2 时间 N²/P + P，故 T(P) = 2N²/P + P，T₁ = 2N²； Speedup = 2N²/(2N²/P + P) = 2P/(2 + P²/N²)；E = T₁/(P·T(P)) = 2N²/(2N² + P²)。 （注意当 N ≫ P 时 P²/N² → 0，加速比 → P，与讲义”speedup → P when N » P”一致；而当步 2 完全串行时 T(P) = N²/P + N²，加速比 = 2/(1+1/P) ≤ 2。） (b) 由 E ≥ 0.95：2N² ≥ 0.95(2N² + P²) ⇒ 0.1N² ≥ 0.95P² ⇒ N ≥ √9.5·P = 3.082P。P = 1024 → N ≥ 3156（验证 N = 3155：E = 2×9,954,025/(19,908,050+1,048,576) = 0.94997，刚好略低于 95%；N = 3156 时 E = 0.9500）。物理含义：串行合并的那 P 步必须小于总工作量的 5%，也就是像素数必须远大于处理器数。 (c) T(P) = 2N²/P + log₂P，P = 1024 → log₂P = 10；E = 2N²/(2N² + P·log₂P) = 2,097,152/(2,097,152 + 10240) = 0.9951 → 99.5%。这正是”用树形归约把 P 的线性开销换成 log P”的价值：只要肯多花一点编排（orchestration）的功夫，弱可扩展性就完全不同。

问题 2（消息传递的死锁与修复） 某个用阻塞 send/recv 实现的网格求解器里，每个进程在每轮开始时都先 sendDown()（向下邻居发自己第一行数据），然后 recvDown()。 (a) 为什么 P ≥ 3 时这个程序常常挂住？请画出等待关系。 (b) 讲义给出的修复是”偶数号线程先 send 后 recv，奇数号线程先 recv 后 send”。为什么这样就不会死锁？ (c) 改用非阻塞 MPI_Isend/MPI_Irecv 后还需要这个奇偶序吗？非阻塞版本还有哪些额外收益与新的约束？

【答案】 (a) 阻塞 send() 的返回条件是”发送方收到确认：消息数据已驻留在接收方的地址空间”（同步 send 的语义）。如果所有进程都先执行 sendDown()，那么每个进程都在等自己的下邻居执行 recvDown()；而下邻居正忙着等它自己的下邻居……等待关系构成一个环（T0→T1→T2→…→T_{P-1}→T0），环上没有任何进程会推进到 recv，因此永久挂起（hang），而不是报错。若消息足够小、MPI 实现恰好用了 eager/buffered 协议，send 会立即返回，程序”看起来正常”——这就是这类 bug 最危险的地方：依赖实现细节、时好时坏。 (b) 相邻两个进程构成的每一对 (偶, 奇) 中，偶数号先发、奇数号先收：偶数的 sendDown 立刻被奇数的 recv 匹配完成，随后奇数再发、偶数再收。于是”等待环”被打破——每一对里总有一方先开出可被对方匹配的接收动作，不会出现全员僵在 send 上的局面。（前提是交换的两侧都按同一约定排列，即偶数在它的每个邻居方向上都是”先发方”，这正是”偶/奇”这个全局一致编号带来的好处。） (c) 不再需要奇偶序：MPI_Isend/MPI_Irecv 立即返回句柄，收发请求由 MPI 进度引擎在后台完成，最后 MPI_Waitall 统一等待，不存在”必须由对方先 recv 才能返回”的循环依赖。额外收益是通信与计算重叠——可以先发/收 ghost 行，然后计算不依赖 ghost 的内部行，最后 Waitall 并计算依赖 ghost 的边界行（本讲 §3.3 示例中把 Waitall 提前了，可以进一步优化）。新的约束是：在 Wait/checksend 确认发送完成之前，绝对不能改动发送缓冲区（否则发出的可能是被改写后的数据，产生间歇性错误）；同样，在确认收到之前不能读接收缓冲区。此外还多出了请求句柄管理、进度引擎与资源占用的开销。

问题 3（为什么”改算法”是合法的并行化手段？红黑 vs 对角波前 vs 三副本技巧） (a) 对角波前（diagonal wavefront）与红黑着色（red-black coloring）分别给出了多大的并行度和多少次同步？ (b) 红黑重排改变了浮点运算的结果序列，凭什么说它是合法的？ (c) 讲义用 diff[3] 轮转把三个屏障减到一个。为什么必须是 3 个副本而不是 2 个？这个技巧的一般名字是什么？

【答案】 (a) 对角波前：一次 sweep 有 2N−1 条反对角线（k = i+j），第 k 条的长度 ≈ min(k, N)，最大并行度 ≈ N/2；但必须在每条波前之后同步一次，所以一次 sweep 要 2N−1 次同步，Span = Θ(N)，而且两端（三角形区域）并行度低、同步频繁——”好消息是并行性存在，坏消息是很难利用”。红黑着色：把单元分成两个相位，相位内零依赖 → 每个相位并行度 Θ(N²/2)，一次 sweep 只需 2 次相位同步（加收敛判定所需的同步），Span = Θ(1)（不含同步）。代价是收敛所需的迭代次数通常变多（同一容差下要多做几次 sweep，因为改变的是迭代格式），即用”多一点工作量”换”高得多的并行度和少得多的同步”。 (b) 合法性来自领域知识：红黑排序仍然是解同一个线性系统的一致迭代格式，它迭代到同一个不动点，只是收敛路径不同；浮点中间值不同、每次迭代的残差不同，但最终仍收敛到误差阈值之内。讲义明确指出：”we needed domain knowledge of Gauss-Seidel method for solving a linear system to realize this change is permissible for the application”。反过来说，如果应用对迭代序列本身敏感（例如需要与某个参考实现逐位一致、或迭代中途被外部观测），这种重排就是不合法的——所以”改算法换并行度”必须由懂该问题语义的人做决定，而不是编译器能自动完成的优化。 (c) 因为需要打破的依赖是”本轮还在读 diff，下一轮却要把它清零“。设迭代 k 使用 diff[k mod m]：线程 A 在迭代 k 里清零某个缓冲区时，线程 B 可能还停在迭代 k−1 的收敛判定上，正在读 diff[(k−1) mod m]；同时 A 自己还要往 diff[k mod m] 里累加。要求这三者互不相同，即 k、k−1、k+1 在模 m 下两两不同 → m ≥ 3。用 2 个副本时 (k+1) mod 2 = (k−1) mod 2，清零会毁掉别人正在读的值，因此仍然需要一个额外屏障，等于没省。这个技巧的一般思想是”用空间（额外副本、轮转缓冲）换依赖的消除“（trade off footprint for removing dependencies），也是并行编程中双缓冲/多缓冲（double buffering、ring buffer）的同一套路。注意副本之间最好填充到不同的 cache line（见 §6 伪共享），否则省下的屏障会被 cache line 弹跳吃回去。

附：本讲术语速查（速查表格，附在七个章节之后）