从Amdahl定律到GPU加速：并行计算性能优化路线图

简介：这是一份由Manning于2021年出版的《Parallel and High Performance Computing》电子书PDF，面向具备一定编程基础、希望系统掌握并行与高性能计算的开发者、科研人员及高年级学生。内容围绕多核CPU与GPU并行展开，从为何并行、并行化规划、性能上限与剖析，到数据设计、性能模型、并行算法与模式，并延伸到向量化、OpenMP、MPI等实践主题，帮助读者建立从问题分析到性能调优的完整思路。压缩包仅含1个PDF文件，大小约35.83MB，便于在电脑或平板端检索阅读。该资料已有692人学习，适合需要补足并行计算知识体系、优化科学计算或工程程序性能的读者作为案头参考。

1. 一张图看清全局：0.8% 的串行决定了 125 倍的天花板

先记住一组数字：Serial 0.8%，Parallel 99.2%。这是并行高性能计算里最该先看的一张饼图，也是整条路线最反直觉的结论——按 Amdahl 定律，就算把 99.2% 的代码都并行化，加速比天花板也只有 1/0.008 = 125 倍，堆多少 CPU 核、插几块 GPU 都不改变这个上限。parallel 从来不是「并行比例高就线性变快」的游戏。

把这条链摊开，是一张四层地图：底层是性能定律与测量（Amdahl、Gustafson-Barsis、机器平衡、roofline、profiling），回答「理论上能快多少、现在卡在哪」；第二层是数据与算法（缓存三 C、AoS/SoA/AoSoA、空间哈希、前缀和），回答「数据怎么排布，硬件才吃得下」；第三层是 CPU 三件套（向量化、OpenMP、MPI），负责从单核推到整机再到多机；最上面一层是 GPU（架构、编程模型、指令式与原生语言、profiling 工具），负责另一个数量级的并行度。这一路适合已经能写 C/C++ 或 Fortran 科学计算代码、想把程序从「能跑」推到「跑满机器」的人，也适合复习并行计算时手边缺一张总图的人。

2. 上限、布局与算法：先算清楚能快多少，再决定怎么快

这一章是整张地图的地基。很多翻车不是因为代码写得差，而是因为动手之前没算过账：不知道上限在哪、不知道程序卡在算力还是带宽、不知道自己选的数据布局在 SIMD 和 GPU 上会变成什么形状。把这几个数先算出来，后面的向量化、OpenMP、MPI、GPU 才有取舍依据。

2.1 Amdahl 与 Gustafson-Barsis：两个方向相反的答案

Amdahl 定律描述的是「问题规模固定、只加并行单元」的场景：设可并行部分占比为 f，并行单元数为 p，加速比 S(p) = 1 / ((1-f) + f/p)。当 p 趋于无穷，加速比收敛到 1/(1-f)——串行比例是唯一的天花板。0.8% 的串行对应 125 倍，这个数在写第一行代码之前就该算出来，它会直接告诉你「值不值得为最后 1% 并行化投入两周」。

Gustafson-Barsis 换了个前提：不固定问题规模，而是在固定时间内把问题做大。此时 S = s + p(1-s)，其中 s 是并行执行中串行部分所占比例。它给出的结论乐观得多——并行单元越多，能解的问题越大，加速比可以随 p 线性增长。两个定律并不矛盾，一个回答「同样的活能多快干完」，一个回答「同样的时间能多干多少活」，选哪个取决于你的应用是强扩展（fixed-size）还是弱扩展（scaled-size）需求。

最容易忽略的是两个公式里符号含义不同：Amdahl 里的 f 是「可并行比例」，Gustafson 里的 s 是「串行比例」，复习时如果混用，算出来的数会差一个量级。下面这段代码把两种口径都算了一遍，可以直接拿来验算自己项目的上限。

# Amdahl：固定问题规模，看加速比上限 def amdahl(parallel_fraction, nproc): serial = 1.0 - parallel_fraction return 1.0 / (serial + parallel_fraction / nproc) # Gustafson：固定时间、放大规模，s 是并行执行时的串行比例 def gustafson(serial_fraction, nproc): return serial_fraction + (1.0 - serial_fraction) * nproc f = 0.992 # 可并行比例，对应封面那张 99.2% 的图 for p in (1, 8, 64, 1024, 10**6): print(f"p={p:>8} Amdahl={amdahl(f, p):7.2f}x Gustafson={gustafson(1-f, p):8.1f}x") print("Amdahl 极限:", 1.0 / (1.0 - f))

参数说明：parallel_fraction是可并行部分占比，nproc是并行单元数（进程或线程），serial_fraction是 Gustafson 口径下的串行比例。输出会看到 Amdahl 那一列很快贴着 125 停住，而 Gustafson 那一列随 p 线性上扬——这就是两条曲线分道扬镳的直观样子。实际测加速比时还要注意：分母必须是你实际会用的串行版本，不是未经优化、带调试符号、跑在单核上的那个版本，否则加速比会虚高。

2.2 机器平衡与 Roofline：算力和带宽谁先见底

机器平衡（machine balance）的定义很简单：峰值浮点算力除以峰值内存带宽，单位是 flop/byte。它回答的是「每从内存搬一个字节，机器能做多少次浮点运算」。某台 64 核机器，每核每周期 16 个浮点操作、频率 2.5 GHz，峰值约 2.56 TFlops；实测 STREAM 带宽 200 GB/s，平衡点就是 12.8 flop/byte。应用的算术强度低于这个数，它就是带宽受限，加核没用，改数据布局、做分块、提高复用率才有用。

Roofline 模型把这个判断画成一张图：性能 = min(峰值算力, 带宽 × 算术强度)。左侧斜线是带宽墙，右侧水平线是算力墙，拐点就是机器平衡点。稠密矩阵乘的算术强度随矩阵边长 N 线性增长（N³ 次运算搬 N² 个数据），N=1000 时约 1000 flop/byte，远在拐点右侧，所以它吃算力；而 STREAM triad 这类 a[i] = b[i] + s*c[i] 的循环，每搬 12 字节做 2 次浮点运算，算术强度不到 0.2，怎么优化都只能贴着带宽墙走。

# 机器平衡点：决定你的应用是算力受限还是带宽受限 cores, freq, flop_per_cycle = 64, 2.5e9, 16 peak_flops = cores * freq * flop_per_cycle # 峰值浮点算力，单位 flop/s peak_bw = 200e9 # STREAM 实测带宽，单位 B/s balance = peak_flops / peak_bw # 机器平衡点，单位 flop/byte print(f"峰值 {peak_flops/1e12:.2f} TFlops, 平衡点 {balance:.1f} flop/byte") for name, ai in (("stream triad", 0.17), ("3D 七点模板", 1.0), ("稠密矩阵乘 N=1000", 1000)): roof = min(peak_flops, peak_bw * ai) # Roofline 的两条边取小 print(f"{name:<20} 算术强度 {ai:>8.2f} 可达上限 {roof/1e12:>6.2f} TFlops")

这段的关键在于flop_per_cycle要按你机器实际支持的 SIMD 宽度和 FMA 能力去填，填大了算出来的平衡点会偏乐观。ai是算术强度，用应用的浮点运算次数除以访存字节数估算即可，粗估就够用——判断「该往哪个方向优化」不需要精确到小数点后两位。

2.3 数据布局：AoS、SoA、AoSoA 与缓存三 C

缓存失效分三类，业内叫三 C：强制失效（compulsory，第一次访问必然缺失，躲不掉）、容量失效（capacity，工作集超过缓存容量）、冲突失效（conflict，多个地址映射到同一组）。对应的补法也很固定：强制失效靠预取和分块掩盖；容量失效靠分块（tiling）让工作集塞进 L1/L2；冲突失效靠填充（padding）错开地址，比如把数组长度补成 2 的幂加一个缓存行的偏移。

AoS 与 SoA 的选择直接决定向量化和 GPU 合并访问的效率。AoS（Array of Structures）把同一对象的字段放在一起，标量代码一次取全字段很方便，但 SIMD 指令要的是「同一字段的连续多个值」，AoS 下这些值跨着结构体，需要 gather，性能大打折扣。SoA 把每个字段单独成数组，正好匹配向量寄存器和 GPU 的合并访存，代价是标量代码的局部性变差。AoSoA 是折中：每 8 或 16 个对象打成一包，包内按字段连续排列，包与包之间连续，既保住向量宽度又保住缓存行利用率。

/* AoS：一次访问拿到一个粒子的 x/y/z/w，适合标量逻辑 */ typedef struct { float x, y, z, w; } Particle; Particle particles[N]; /* SoA：同一字段连续，向量化和 GPU 合并访存的最爱 */ typedef struct { float x[N], y[N], z[N], w[N]; } ParticleSoA; /* AoSoA：8 个粒子一包，包内按字段排，匹配 32 字节的 SIMD 宽度 */ #define LANE 8 typedef struct { float x[LANE], y[LANE], z[LANE], w[LANE]; } Pack; Pack blocks[N / LANE]; /* 遍历时先走 blocks，再走 LANE */

三个布局不是「谁更先进」的问题，而是「谁匹配你的内层循环」。判断方法很土但有效：看内层循环每一轮实际读了哪些地址，如果这些地址跨着结构体跳，就该换 SoA 或 AoSoA；如果内层循环本来就要一次用完整对象的四个字段，AoS 反而更省带宽。改布局之前先 profile，别凭直觉重写——这是血泪经验里最常见的一条。

2.4 算法模式：空间哈希、前缀和与全局求和

并行算法的核心不是把串行算法逐行加#pragma，而是换一种「局部算、边界换、整体归约」的组织方式。空间哈希是典型例子：粒子或网格点找邻居，暴力法是 O(N²)，按相互作用半径划 cell、把 cell 索引哈希到桶里，就降到接近 O(N)。哈希函数选得差会让桶退化成链表，所以要么用完美哈希（键空间小且已知），要么用紧凑哈希加质数取模；cell 大小取相互作用半径，保证只需要查相邻 27 个桶。

前缀和（scan）是并行算法的另一块基石，排序、流压缩、稀疏矩阵装配都靠它。串行 scan 是 O(N) 的简单循环，并行版有两条路：step-efficient 版本做 log N 轮、每轮全量读写，总工作量 O(N log N)；work-efficient 版本用上下扫的两趟树形结构，总工作量 O(N)。数组大到装不进一个工作组时，还要分三段：块内 scan、块间 scan、把块间结果加回各块偏移。

全局求和看着最简单，坑却最深：浮点加法不满足结合律，(a+b)+c 和 a+(b+c) 的结果在末位就不一样。线程数、进程数、归约树形状一变，结果就变，调试时会出现「同样的输入两次跑出不同结果」的玄学现象。补法是固定归约顺序（比如按固定树形做两阶段归约）、用补偿求和（Kahan/Neumaier）降低误差，或者直接把结果按可复现的要求对齐到某个精度。

2.5 从 benchmark 到 profiler：把理论上限和实测对上

顺序不能反：先测硬件，再测应用，最后才优化。硬件侧要拿到四个数——核数与频率、峰值浮点算力、实测内存带宽、缓存层级与容量；应用侧要拿到三个数——热点函数占比、向量化比例、通信或数据传输占比。两边一对，就知道自己离 roofline 的那条线还差多远，也就知道优化空间到底在算力、带宽还是同步开销上。

profiling 工具分三类：采样式（低开销，看热点分布）、硬件计数器（看 IPC、缓存命中率、向量化指令占比）、专项检查（线程竞争、内存越界）。采样先跑一遍拿到粗粒度画像，再用硬件计数器确认「热点函数到底是被算力卡住还是被访存卡住」，两者结论经常不一致，以计数器为准。容易忽略的是编译选项也要记录：同一份代码，-O2和-O3 -march=native能差出几倍，优化前后如果不锁定编译选项，测出来的加速比没有意义。

3. CPU 主战场：向量化、OpenMP 与 MPI 的三层分工

CPU 这一侧的并行是三层嵌套的：最内层是单核里的 SIMD 向量化，中间层是节点内的多线程（OpenMP），最外层是跨节点的消息传递（MPI）。三层各管一段，混着用才对，只做其中一层通常只能拿到一小部分性能。这一章的复习重点是「每层解决什么问题、边界在哪、什么时候该往下一层走」。

3.1 向量化：从自动向量化到 intrinsics 的五级台阶

向量化的收益是白捡的——同一颗核，宽度从 4 个单精度变成 16 个单精度，理论算力直接乘 4，不增加任何线程开销。代价是代码约束变多。按投入从低到高可以排成五级：

绝大多数项目停在第二、三级就够了。自动向量化失败的原因很集中：循环里有跨迭代依赖、指针别名让编译器不敢假设不重叠、访问步长不是 1、循环体里有函数调用或分支。前三类用restrict加omp simd就能解决大半。

// 三个关键点：连续访存、restrict 消歧义、omp simd 给编译器明确许可 void axpy(int n, float a, const float * restrict x, float * restrict y) { #pragma omp simd for (int i = 0; i < n; ++i) { y[i] = a * x[i] + y[i]; // 单条 FMA 完成一次乘加 } }

restrict向编译器承诺 x 和 y 指向不重叠的内存，没有这个承诺，编译器必须假设写 y[i] 可能改到 x[i+1]，于是放弃向量化。#pragma omp simd是可移植的写法，比编译器私有的指示更稳。循环后面的逻辑说明很简单：这段代码做的是 y = a*x + y，但它的价值不在算术，而在于它是一个「编译器能证明可以打包成向量指令」的形状。

验证是否真的向量化，要看编译器的优化报告（-fopt-info-vec、-Rpass=loop-vectorize之类），别只看运行时间变快了就以为成了——变快可能只是缓存效应。

3.2 OpenMP：循环级、函数级与高层并行的取舍

循环级 OpenMP 是最常见的入门用法：在 for 前面加#pragma omp parallel for，一行搞定。它的问题是每次进入循环都要 fork/join 一次线程组，循环体太短时，线程创建和同步的开销会吃掉全部收益。经验阈值是循环体工作量要足够大，粗略估计：循环体耗时要远大于一次 fork/join 的微秒级开销。

函数级（高层）OpenMP 换个思路：把整个时间步或整个迭代包在一个parallel区里，内部用for、single、barrier分工，只在最外层 fork 一次。这样每步省下的 fork/join 累积起来相当可观，代价是代码结构要重写，parallel区内的变量作用域要重新梳理。

// 循环级：写法最简单，但每进入一次循环就 fork/join 一次 #pragma omp parallel for for (long i = 0; i < n; ++i) a[i] = b[i] + scalar * c[i]; // 函数级：整个时间步只 fork 一次，适合迭代上千步的模拟 #pragma omp parallel { for (int step = 0; step < nsteps; ++step) { #pragma omp for for (long i = 0; i < n; ++i) a[i] = b[i] + scalar * c[i]; #pragma omp single { swap(a, b); } // 只有一条线程做交换，其余线程等 barrier #pragma omp barrier } }

两段代码的差别在 fork 次数上：第一段若放在时间步循环里，就是 nsteps 次 fork/join；第二段只有一次。single后面必须跟隐式或显式的 barrier，否则其他线程可能在交换完成前就进入下一轮，读到半新半旧的数组。这个顺序错一次，结果就是概率性错误，跑十次才对八次。

3.3 变量作用域、reduction 与伪共享：正确性的第一道关

OpenMP 的作用域规则不复杂，但违反它的后果很严重：parallel区外定义的变量默认共享，区内定义的变量默认私有；for的循环迭代变量自动私有；需要跨迭代累加的变量必须用reduction子句声明。最典型的错法是把累加变量写成共享。

// 错误：sum 是共享变量，多个线程同时读改写，结果随机偏小 double sum = 0.0; #pragma omp parallel for for (int i = 0; i < n; ++i) sum += x[i]; // 正确：reduction 给每个线程一份私有副本，最后按子句指定的运算符合并 double total = 0.0; #pragma omp parallel for reduction(+:total) for (int i = 0; i < n; ++i) total += x[i];

reduction(+:total)的语义是：每个线程各自累加到私有副本，退出时按 + 合并。要注意两点：合并顺序不确定，浮点结果不可复现；如果你的结果需要逐位可复现，得自己做固定顺序的两阶段归约。

比数据竞争更隐蔽的是伪共享：两个线程写同一缓存行里的不同变量，硬件层面反复失效，性能掉得莫名其妙但结果全对。补法是给每个线程的计数器加填充，或者用#pragma omp parallel for配合线程本地的对齐数组，让每个线程的私有数据独占缓存行。

3.4 MPI：六个函数、集合通信与自定义数据类型

任何 MPI 程序都由六个函数撑起来：MPI_Init、MPI_Comm_size、MPI_Comm_rank、MPI_Send、MPI_Recv、MPI_Finalize。点对点通信的难点是死锁：两个进程都用阻塞式MPI_Send互相发送大消息，双方都在等对方收，程序就挂住。补法是调整收发顺序（先奇数进程发、偶数进程收）、用MPI_Sendrecv合并成一次调用，或者直接改用非阻塞的MPI_Isend/MPI_Irecv加MPI_Waitall。

集合通信是 MPI 真正省代码的地方：Bcast分发小规模输入、Reduce汇总单值、Gather/Scatter做数据分发与收集、Barrier同步计时点。最常踩的坑是把集合通信写在条件分支里——集合通信要求通信域内所有进程都调用，少一个就死锁。任何if (rank == 0) { MPI_Bcast(...) }的写法都是错的。

自定义数据类型（MPI_Type_create_struct、MPI_Type_vector）是另一个性能点：把「结构体数组里只发某个字段」这种非连续数据打包成一次通信，比手动拷进临时缓冲区再发要快，代码也更短。规则是先用MPI_Get_address算偏移、建类型、MPI_Type_commit提交，用完MPI_Type_free释放。三维模板计算里的 ghost cell 交换就是它的典型应用：六个面、每个面是非连续的一块，用向量类型描述一次就能发完。

3.5 混合 MPI+OpenMP：什么时候真的值得上

混合编程的收益来自两点：一是减少通信域规模，把节点内的线程协作交给共享内存，只有节点间才走网络；二是减少内存副本，每个节点只存一份数据被多线程共享，而不是每个进程各存一份。代价是复杂度上升：线程安全的 MPI 需要MPI_Init_thread指定级别，FUNNELED只允许主线程调 MPI，MULTIPLE允许多线程同时调但要承担实现质量风险，多数项目选FUNNELED就够了。

什么时候不值得：纯 MPI 版本已经能扩展到目标规模、通信占比不高、或者团队里没人能维护混合代码。混合编程的调优空间很大但也很玄学，绑核方式、线程亲和性、通信线程与计算线程的争抢，任何一项没配对都可能比纯 MPI 还慢。判断标准很简单：先用纯 MPI 跑出 scaling 曲线，如果通信占比超过 30% 再考虑混合，否则先优化算法和通信模式。

4. GPU 加速：SM、占用率与三种落地方式的选择

GPU 的并行度和 CPU 不在一个量级，但它的脾气也完全不同：CPU 靠大缓存和乱序执行隐藏延迟，GPU 靠海量线程切换隐藏延迟；CPU 的核数是个位到百位，GPU 的线程数是十万级。理解这一点，才能理解为什么 GPU 上的「优化」大部分是在讲「访存是否合并」和「占用率够不够」。

4.1 硬件模型：SM、SIMT 与占用率

GPU 的计算单元叫流式多处理器（SM，不同厂商叫法有差异，有的叫子切片），一个 SM 里有一批处理单元（PE），同一时刻以 SIMT 方式执行同一段指令的一批线程称为一个 warp（或 wavefront），宽度通常是 32 或 64。warp 内的线程如果走不同分支，会串行执行各分支，这是 GPU 上分支代价高的根本原因。

占用率（occupancy）是核心指标：一个 SM 上实际驻留的活跃 warp 数除以硬件支持的最大 warp 数。它受三个资源限制——每线程寄存器数、每工作组共享内存字节数、工作组大小。寄存器用得越多，能同时驻留的 warp 越少，占用率越低，隐藏访存延迟的能力越差。这不是说占用率越高越好，而是说低占用率必须由高指令级并行或高数据复用率来补偿，否则流水线就在空转。

峰值算力的估算方法和 CPU 一致：PE 数 × 每 PE 每周期操作数 × 频率。拿这个数和实测带宽一比，照样能画出 GPU 的 roofline，判断某个 kernel 是算力受限还是带宽受限。千万别拿厂商标称的峰值直接对比实测性能，中间隔着占用率、访存模式和指令混合三重损耗。

4.2 编程模型：NDRange、工作组与三层内存

GPU 编程模型的抽象很统一：把问题切成一个 NDRange（网格），网格由工作组（block）组成，工作组由工作项（线程）组成。全局索引 = 工作组编号 × 工作组大小 + 组内编号。这个映射关系是所有 GPU 代码的第一行，也是所有越界错误的源头。

内存分三层：每个工作项私有的寄存器与局部内存、工作组内共享的局部内存（片上）、全局内存（显存）。工作组内可以用屏障同步，工作组之间没有 barrier，跨组只能靠「拆成两个 kernel 分两次启动」或者原子操作。归约模式因此必须两阶段：组内归约到共享内存，写出到全局，再启一个 kernel 做组间归约。

// 每个线程处理一个元素；全局索引 = 组号 × 组大小 + 组内编号 __global__ void saxpy(int n, float a, const float *x, float *y) { int i = blockIdx.x * blockDim.x + threadIdx.x; // 拼出全局下标 if (i < n) { // 边界判断不能省 y[i] = a * x[i] + y[i]; // 相邻线程访存连续，可合并 } } // 主机侧启动：block = 256，grid = (n + 255) / 256，向上取整保证覆盖全部元素

这段代码值得逐行看：blockIdx.x * blockDim.x + threadIdx.x是固定套路；if (i < n)是防止最后一组越界，删掉它跑小数组时可能看不出问题，跑大数组必然崩；相邻线程访问相邻地址，硬件才能合并成一次大事务，如果改成y[i * stride]，合并访存立刻退化，性能掉一个数量级。

4.3 三种落地方式：指令式、原生语言与可移植层

给已有 CPU 代码加速，有三条路可走，投入和上限差别很大：

实操路径通常是：先用指令式把代码跑起来，量一下端到端能快多少；如果收益够，就停在这儿；如果不够，再把最热的 kernel 用原生语言重写，其余部分保留指令式。一上来就全量重写是最常见的翻车方式——花两周重写完，发现瓶颈其实在数据搬运上。

4.4 数据搬运：PCIe 带宽、异步队列与重叠

GPU 加速最容易被低估的成本是数据搬运。CPU 到 GPU 的总线带宽，Gen3 x16 单向理论约 16 GB/s，实测通常只有理论值的六七成；而 GPU 自身的显存带宽是数百 GB/s 量级，两者差了一个数量级以上。这意味着一个算术强度低的 kernel，搬数据的 time 比算数据的 time 还长，加速比自然是负数。

三条缓解手段：一是把数据尽量留在显存里，只在初始化和输出时搬运；二是用页锁定内存（pinned memory）提高传输带宽，普通可分页内存的传输要走中转缓冲；三是用异步拷贝加多队列，让第 k 批的计算和第 k+1 批的传输重叠。重叠能不能生效，取决于你的 kernel 是否分成了足够多的批次，以及队列是否真的并发。

判断「值不值得搬」有个粗略公式：计算时间 / 传输时间 > 3 再考虑上 GPU。传输时间 = 数据量 × 2 / 实测带宽（进和出各一次），计算时间按 roofline 估。这个比值算出来小于 1 的时候，任何 kernel 优化都救不回来。

4.5 GPU profiling：先找热点，再决定搬不搬

GPU 侧的 profiling 工作流和 CPU 侧一致，只是指标换了：kernel 耗时占比、占用率、显存吞吐利用率、warp 停滞原因（等内存、等指令、等同步）。顺序是：先在 CPU 上 profile 找到热点，判断它是否适合 GPU（并行度高、分支少、算术强度够）；再上指令式版本跑通；再加数据搬运指令；最后做引导式分析看瓶颈在哪。

最容易犯的错是跳过 CPU profiling 直接上 GPU：把整个程序搬上去，发现端到端只快了 5%，然后开始怀疑 GPU 不行。实际原因往往是 90% 的时间还在 CPU 段的串行代码里，GPU 那部分根本不在关键路径上。先 profile 再搬，能省掉大量无用功。

5. 易错点：五处最容易翻车的细节

前面四章讲的是「应该怎么做」，这一章讲「实际会怎么错」。这些

本文还有配套的精品资源，点击获取