从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(fp{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² 个数据N1000 时约 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), (稠密矩阵乘 N1000, 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 合并访问的效率。AoSArray 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; /* AoSoA8 个粒子一包包内按字段排匹配 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、把块间结果加回各块偏移。全局求和看着最简单坑却最深浮点加法不满足结合律(ab)c 和 a(bc) 的结果在末位就不一样。线程数、进程数、归约树形状一变结果就变调试时会出现「同样的输入两次跑出不同结果」的玄学现象。补法是固定归约顺序比如按固定树形做两阶段归约、用补偿求和Kahan/Neumaier降低误差或者直接把结果按可复现的要求对齐到某个精度。2.5 从 benchmark 到 profiler把理论上限和实测对上顺序不能反先测硬件再测应用最后才优化。硬件侧要拿到四个数——核数与频率、峰值浮点算力、实测内存带宽、缓存层级与容量应用侧要拿到三个数——热点函数占比、向量化比例、通信或数据传输占比。两边一对就知道自己离 roofline 的那条线还差多远也就知道优化空间到底在算力、带宽还是同步开销上。profiling 工具分三类采样式低开销看热点分布、硬件计数器看 IPC、缓存命中率、向量化指令占比、专项检查线程竞争、内存越界。采样先跑一遍拿到粗粒度画像再用硬件计数器确认「热点函数到底是被算力卡住还是被访存卡住」两者结论经常不一致以计数器为准。容易忽略的是编译选项也要记录同一份代码-O2和-O3 -marchnative能差出几倍优化前后如果不锁定编译选项测出来的加速比没有意义。3. CPU 主战场向量化、OpenMP 与 MPI 的三层分工CPU 这一侧的并行是三层嵌套的最内层是单核里的 SIMD 向量化中间层是节点内的多线程OpenMP最外层是跨节点的消息传递MPI。三层各管一段混着用才对只做其中一层通常只能拿到一小部分性能。这一章的复习重点是「每层解决什么问题、边界在哪、什么时候该往下一层走」。3.1 向量化从自动向量化到 intrinsics 的五级台阶向量化的收益是白捡的——同一颗核宽度从 4 个单精度变成 16 个单精度理论算力直接乘 4不增加任何线程开销。代价是代码约束变多。按投入从低到高可以排成五级级别做法改动量典型收益适用场景1换优化库BLAS/LAPACK 类极小高标准线性代数、FFT2编译器自动向量化无中到高简单连续访存循环3编译指示omp simd、ivdep小中编译器保守、存在可证明无依赖的循环4intrinsics 内建函数大高复杂分支、需要精确控制数据打包5手写汇编极大有限极端热点、库函数级优化绝大多数项目停在第二、三级就够了。自动向量化失败的原因很集中循环里有跨迭代依赖、指针别名让编译器不敢假设不重叠、访问步长不是 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[i1]于是放弃向量化。#pragma omp simd是可移植的写法比编译器私有的指示更稳。循环后面的逻辑说明很简单这段代码做的是 y a*x y但它的价值不在算术而在于它是一个「编译器能证明可以打包成向量指令」的形状。验证是否真的向量化要看编译器的优化报告-fopt-info-vec、-Rpassloop-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 混合 MPIOpenMP什么时候真的值得上混合编程的收益来自两点一是减少通信域规模把节点内的线程协作交给共享内存只有节点间才走网络二是减少内存副本每个节点只存一份数据被多线程共享而不是每个进程各存一份。代价是复杂度上升线程安全的 MPI 需要MPI_Init_thread指定级别FUNNELED只允许主线程调 MPIMULTIPLE允许多线程同时调但要承担实现质量风险多数项目选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 256grid (n 255) / 256向上取整保证覆盖全部元素这段代码值得逐行看blockIdx.x * blockDim.x threadIdx.x是固定套路if (i n)是防止最后一组越界删掉它跑小数组时可能看不出问题跑大数组必然崩相邻线程访问相邻地址硬件才能合并成一次大事务如果改成y[i * stride]合并访存立刻退化性能掉一个数量级。4.3 三种落地方式指令式、原生语言与可移植层给已有 CPU 代码加速有三条路可走投入和上限差别很大方式代表改动量可移植性性能上限适合谁指令式OpenACC、OpenMP target小加编译指示较好中已有大段循环、想快速验证 GPU 收益原生语言CUDA、HIP大重写热点 kernel差绑定厂商高热点集中、需要精细控制寄存器与共享内存开放标准OpenCL、SYCL中大好中到高需要跨多种加速器部署可移植层Kokkos、RAJA中改写成模板化循环很好中到高长期维护、要同时支持多后端实操路径通常是先用指令式把代码跑起来量一下端到端能快多少如果收益够就停在这儿如果不够再把最热的 kernel 用原生语言重写其余部分保留指令式。一上来就全量重写是最常见的翻车方式——花两周重写完发现瓶颈其实在数据搬运上。4.4 数据搬运PCIe 带宽、异步队列与重叠GPU 加速最容易被低估的成本是数据搬运。CPU 到 GPU 的总线带宽Gen3 x16 单向理论约 16 GB/s实测通常只有理论值的六七成而 GPU 自身的显存带宽是数百 GB/s 量级两者差了一个数量级以上。这意味着一个算术强度低的 kernel搬数据的 time 比算数据的 time 还长加速比自然是负数。三条缓解手段一是把数据尽量留在显存里只在初始化和输出时搬运二是用页锁定内存pinned memory提高传输带宽普通可分页内存的传输要走中转缓冲三是用异步拷贝加多队列让第 k 批的计算和第 k1 批的传输重叠。重叠能不能生效取决于你的 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. 易错点五处最容易翻车的细节前面四章讲的是「应该怎么做」这一章讲「实际会怎么错」。这些本文还有配套的精品资源点击获取