GEMM 案例研究：四步优化方法论
本文所有数字来自
examples/09-matrix-multiply/（Ryzen 7 5800H / WSL2 实测），代码在规范库include/hpc/gemm.hpp。GEMM（C = A·B）是性能优化的标准载体：内存层次、SIMD、多核三个瓶颈依次显现，每一步优化都有可验证的物理原因。
为什么选矩阵乘法
- 所有机器层次都会暴露：naive 版卡在缓存，tiled 版暴露 SIMD 宽度，多核版暴露并行开销——一套代码讲完整个优化栈；
- 有清晰的代价模型：n×n 乘法固定 2n³ flops，GFLOPS 可直接对照理论峰值；
- 有工业参照物：BLAS 是"同类最优"的标尺，能看清教学实现与工业实现的差距在哪。
阶段 0：代价模型
n×n GEMM 需要 2n³ 次浮点运算。naive 的 ijk 循环里，内层 k 循环按列访问 B（步长 n 个 float）：每读一个 B[k][j] 只用了缓存行里 1/16 的数据（64B 行 = 16 float），其余 15 个被浪费——这就是"内存层次不友好"的定量含义。
阶段 1 → 2：缓存分块（13–14x）
把矩阵切成 block_size×block_size 的块，块内改用 i,k,j 序：B 的行和 C 的行都变成连续访问。块尺寸的选择是容量算术：三个 64×64 float 块 = 3×16 KB = 48 KB，恰好压在 32 KB L1 与 512 KB L2 之间——块再大就溢出 L1/L2，再小则分块开销占比升高。
实测（512×512）：0.82 → 11.5 GFLOPS，一个数量级来自纯粹的访存修复，没有增加任何计算资源。这是"布局/分块优化通常比指令级优化收益大"的最直观证据。
阶段 2 → 3：SIMD 微内核（+50%，以及为什么不是 8x）
tiled 的内层循环天然适合向量化：固定 A[i][k]，沿 j 方向 C[i][j..j+8] += A[i][k] * B[k][j..j+8]，一条 FMA 处理 8 个 float（AVX2）或 4 个（NEON）。
实测只有 ~1.5x（11.5 → 17.4），远不是"宽度 8 倍"的天真预期。原因：微内核是访存受限的——每个 k 都要把 C 的整行加载进来、FMA 后存回去，2n² 次行级 load+store 的流量压住了 FMA 单元的吞吐。SIMD 只加速了算术部分，搬数据的部分没变。
要真正兑现宽度，需要寄存器分块：一次在寄存器里累加 mr×nr 的 C 子块（如 6×8），把同一份 B 数据复用 mr 次、C 的读写摊薄到 1/(mr·nr)。这正是 BLAS 微内核（GotoBLAS 方法论）的核心，也是教学版与工业版的主要差距。
阶段 3 → 4：OpenMP（以及负优化的尺寸）
按行 #pragma omp parallel for schedule(static) 划分：每行工作量相等，静态划分无调度开销；各线程写不同的行，天然无共享冲突。
实测的两面：
- 512×512：17.4 → 40.0 GFLOPS（~2.3x）——受限于微内核访存，不是核数不够；
- 256×256：18.9 → 5.7 GFLOPS，并行反而慢了 3 倍——16 个线程的创建/同步开销超过了 34 MFLOP 的工作量。
教训可以量化：并行收益 = 工作量 ÷ (并行开销 + 受带宽限制的部分)。工作量盖不住开销时，单线程更快。1024×1024 时 parallel 回到 ~108 GFLOPS，趋势符合预期。
与 BLAS 的距离
教学版刻意不做的事（以及 BLAS 做了的事）：
| 技术 | 作用 |
|---|---|
| Packing（把 A/B 块拷贝进连续缓冲） | 消除 TLB/缓存冲突，保证微内核输入连续 |
| 寄存器分块微内核（mr×nr） | 摊薄 C 的读写，兑现 FMA 吞吐 |
| 对齐分配与对齐 load/store | 避免跨行惩罚 |
| 动态核间划分 | 大小核/负载不均场景 |
实际工程中，矩阵乘法应该直接调用 BLAS（OpenBLAS/MKL）；本案例的价值在于让你看懂这些技术为什么存在，以及在自己的非标准内核里如何复用同一套诊断顺序。
方法论迁移
四步顺序（访存 → 向量化 → 并行 → 算法）适用于任何计算内核，与优化决策手册的诊断清单一一对应。每一步都遵守同一纪律：改一处，测一处，对照参考实现验证正确性——examples/09-matrix-multiply/src/matrix_multiply.cpp 对每个阶段都对照 naive 校验数值，这正是"可验证的性能工程"在案例层的体现。