做底层性能优化的人早晚都要跟矩阵乘法正面碰一次。原因很简单它计算量巨大、访存模式典型、优化空间从内存布局一路延伸到CPU指令集几乎把一台PC的性能瓶颈全部串起来了。这次记录的是一个正经的GEMM优化过程——输入输出都是普通双精度稠密矩阵运行环境就是一台再普通不过的八核桌面PC没有GPU、没有异构加速最终把朴素三重循环从个位数GFLOPS拉到了两百多GFLOPS折算下来大约30倍的性能提升。这篇文章会沿着我实际动手的路径展开先建立基线然后依次做循环重排、缓存分块、寄存器分块、SIMD FMA、数据打包、多线程每一步都有代码片段、性能数据和踩坑记录每个阶段为什么这么做、收益从哪来我都会尽量讲透。不管你平时做图形学、数值计算、机器学习推理还是刚刚接触高性能计算这条优化路线都很值得完整走一遍——它教你的不是某个固定技巧而是当一个程序慢怎么系统地找到瓶颈并把效率榨出来的思路。整个项目最核心的思路其实只有一句话矩阵乘法本身是计算密集的但绝大多数朴素实现会把它变成访存密集的问题。优化GEMM本质上就是让CPU的每一个时钟周期都在做乘加运算而不是在等数据。1. 为什么要跟矩阵乘法死磕它是一切性能优化的“试金石”1.1 问题定义与峰值账本先明确我们要解决的问题。给定三个矩阵A(M×K)、B(K×N)、C(M×N)计算C A×B这里C的每个元素是A的对应行与B的对应列做K次乘加累加的结果。数据布局采用C语言风格的行主序M、N、K三个维度可以任意变化不利用稀疏性、对称性等任何结构。这就是所谓通用矩阵乘GEMM。为了后面讨论方便先把算力账算清楚。两个N×N矩阵相乘计算量是2×N³次浮点运算其中N³次乘法、N³次加法。我的测试矩阵取MNK2048也就是单次运算量是2×2048³ ≈ 17.2 GFLOP。这台测试机的CPU是一颗八核心桌面处理器主频接近4GHz支持AVX2与FMA指令。理论浮点峰值怎么算AVX2的向量寄存器是256位能同时处理4个双精度浮点数FMA指令一次完成乘法和加法等于8次浮点运算现代CPU通常每个核心有两个FMA执行单元每周期可以发射两条FMA。所以单核每周期是16 flops八核心全开就是8×4GHz×16 flops/cycle ≈ 512 GFLOPS。这个512 GFLOPS就是天花板但实际因为睿频降频、内存带宽、指令发射限制等因素能达到理论峰值的一半就已经很强了。我在优化结束时跑到248 GFLOPS大约是理论峰值的48%这是合理且可复现的数字。1.2 朴素三重循环为什么这么慢大部分人第一次写的矩阵乘法长这样void naive_gemm(const double *A, const double *B, double *C, int M, int N, int K) { for (int i 0; i M; i) { for (int j 0; j N; j) { double sum 0.0; for (int k 0; k K; k) { sum A[i * K k] * B[k * N j]; } C[i * N j] sum; } } }逻辑完全正确但性能几乎是灾难。原因在最内层循环访问B[k * N j]时k每增加1地址跳跃N×8字节。N2048时这个跨度是16KB等于每读一个double就把一个缓存行浪费掉——64字节的缓存行只取其中8字节剩下的56字节根本不用然后立刻跳到下一个远得不能再远的内存地址。这会导致什么结果B矩阵的访问模式几乎没有空间局部性每一个B元素都可能触发一次缓存缺失CPU大量的时间在等内存响应而不是在算乘法。A矩阵的访问是连续的情况好一些但整体上CPU的浮点计算单元大多数时候处于空转状态。我实测这个基线版本的性能大约只有7.8 GFLOPS相对于512 GFLOPS的理论峰值利用率连2%都不到。这不是GCC不够聪明是访存模式决定了CPU根本没有足够的喂给计算单元的数据。1.3 先搭好测试与正确性框架动手优化前先把如何测量这件事确定下来否则后面每一步都无法判断是变快还是变慢。性能测量我采用这样的方式矩阵尺寸固定为2048×2048用随机数填充范围[-1, 1]。每次测试前先跑一次目标函数做预热然后用clock_gettime计时连续跑5次取中位数。这样能排除频率波动、页错误和第一次触发冷缓存的干扰。正确性校验用朴素版本的输出作为基准。因为浮点运算的累加顺序不同会导致结果不完全一致所以比较的是最大相对误差阈值放在1e-10量级。FMA指令参与计算后误差会比朴素累加略大但只要在这个量级内就算通过。double max_rel_err 0.0; for (int i 0; i M * N; i) { double denom fmax(fabs(ref[i]), fabs(opt[i])); max_rel_err fmax(max_rel_err, fabs(ref[i] - opt[i]) / denom); } printf(max relative error %e\n, max_rel_err);这一步看起来平平无奇但后面做寄存器分块和打包时会遇到很多性能变好、结果错乱的瞬间没有一个可靠校验框架排查会极其痛苦。2. 第一板斧循环重排让访存从“爬楼梯”变成“逛超市”2.1 ikj循环交换背后发生了什么解决朴素实现最大的问题——B矩阵列访问——不需要高深理论只要把循环顺序换一下。把最内层从k循环改成j循环void ikj_gemm(const double *A, const double *B, double *C, int M, int N, int K) { for (int i 0; i M; i) { for (int k 0; k K; k) { double r A[i * K k]; for (int j 0; j N; j) { C[i * N j] r * B[k * N j]; } } } }这个版本里最内层的j循环同时遍历B矩阵的一整行和C矩阵的一整行两个都是连续内存访问。A矩阵的A[i*K k]被提取到一个标量寄存器r中在整个j循环期间复用。内存访问模式从每步跨16KB变成了每步走8字节顺路扫完一整片配合CPU的硬件预取器效果立竿见影。为什么方案奏效关键在于循环重排改变了数据在内存中的流动方向。原版的ijk顺序B和C的访问模式都糟糕改成ikj后B和C都是行连续硬件预取器可以按顺序把一整块数据提前拉进缓存CPU不再为单个数据等待内存。2.2 为什么编译器开-O3还是不理想有朋友会问GCC开-O3不是能自动优化吗这里有个常见的误解。GCC确实可以对某些简单循环做循环交换但不是总能做对它要保守地分析别名、浮点结合律、循环依赖等面对二维数组指针和未知的矩阵尺寸时很多优化宁可不动。即便做了向量化也会因为C[i][j]的反复读写生成额外的load/store指令无法真正逼近峰值。更重要的是靠编译器自动优化很难实现后面要讲的寄存器分块——那是把一组C值长期保留在寄存器里循环跑完才写回内存。这个模式需要从结构上改变代码组织方式编译器不会替你设计。2.3 循环重排的实测提升循环重排这一步的实测大约从7.8 GFLOPS提升到19.6 GFLOPS接近2.5倍。这个代价几乎为零只是改变了三重循环的顺序没有任何一条额外的机器指令被手动添加。它充分说明很多性能问题的根子不在CPU不够快而在访存模式把CPU饿死了。不过19.6 GFLOPS仍然远低于理论峰值。原因有两个一是C矩阵在j循环里被反复读和写内层每迭代一次都要把C[i][j]从内存取出再写回二是还没有利用向量寄存器每条指令一次只算一个浮点运算。这两点恰好是接下来几步要解决的。3. 第二板斧缓存分块把“远水”变成“近水”3.1 分块思想与经典五重循环循环重排解决的是访问连续性问题但没有解决数据重用问题。前面ikj版本中内层j循环结束后下一次k迭代会把B的新一行重新读一遍C的一整行也要从缓存里重新取。当N2048时C的一行是16KB现代CPU的L1缓存通常32KB勉强能装下但B的一行8KB加上A的数据L1已经被塞满。更糟糕的是外层的i循环每轮都要重新扫描B的全部数据意味着B矩阵被反复从主存读入。缓存分块的思想是不把整个矩阵装进缓存而是切成小方块让每个小方块在某一层缓存中完成大部分计算后再换下一块。这不是什么神秘方法本质上是利用局部性原理——把一个大循环拆成多个嵌套的小循环每次只处理能放进缓存的一小块数据。经典的实现是那个著名的五重循环void blocked_gemm(const double *A, const double *B, double *C, int M, int N, int K) { for (int i0 0; i0 M; i0 MC) for (int k0 0; k0 K; k0 KC) for (int j0 0; j0 N; j0 NC) for (int i i0; i i0 MC; i) for (int k k0; k k0 KC; k) { double r A[i * K k]; for (int j j0; j j0 NC; j) C[i * N j] r * B[k * N j]; } }这层分块后A的一个MC×KC块和B的一个KC×NC块会被固定在某个缓存层级反复使用而不是每访问一次就失效。五重循环本身还有一些排列变体外层是i0、k0、j0的六种排列不同排列对应的缓存命中情况不同。我实测下来i0→k0→j0的排列比较稳妥C块尺寸MC×NC能留在L2A块和B块的流动方向也符合从L3读大块、在L2算小块的预期。3.2 块尺寸选择不是随便试出来的分块尺寸的选择是整个缓存优化中最玄学也最关键的一步。我的初始参数是MC256、KC256、NC4096理由如下NC×KC的B块是256×4096×8字节8MB正好压在L3缓存容量附近让B块在L3里被反复利用MC×KC的A块是256×256×8字节512KB可以放在L2级别MC×NC的C块是256×4096×8字节8MBC块会被多次k迭代复用这个尺寸也适合L3。这不是我拍脑袋定的是参考现代GEMM库的分块策略再结合自己机器缓存层级调整的结果。如果你用的CPU主频和缓存不同最靠谱的做法是写一个参数扫描脚本遍历MC∈{64,128,256,512}、KC∈{64,128,256,512}、NC∈{1024,2048,4096}的笛卡尔积每组跑一遍找最优。这类自动调优在现代库中已经变成标配了人肉调参容易错过好组合。3.3 分块后的性能变化加入缓存分块后性能从19.6 GFLOPS提升到45.3 GFLOPS。这一倍的收益不是来自某个魔法指令而是数据待在正确的位置上A块和B块的重用次数大幅增加主存访问减少计算单元终于有稳定数据流可以消费了。但45.3 GFLOPS离理论峰值还差一个数量级。原因很明确最内层仍然是标量double乘法每个时钟周期最多执行一两条浮点指令向量寄存器完全没有被利用。接下来要动真格的了。4. 第三板斧寄存器分块与FMA逼近峰值4.1 为什么必须手动向量化从标量到向量是GEMM优化里收益最明显的一步。AVX2寄存器是256位一个寄存器能装4个doubleFMA指令一条同时完成4次乘法和4次加法等于8个flops。单条指令的算力是标量的8倍这是硬件白送的不利用就是暴殄天物。但直接对原来的内层循环做向量化有个问题编译器自动生成的向量代码还是会把C反复load和store因为这些操作在语义上是必须的。真正高性能的写法是寄存器分块手动把多个C元素作为局部变量一直留在向量寄存器里内层循环只更新这些寄存器循环结束后一次性写回内存。4.2 4×8微内核代码我选择的基本计算单元是4×8微内核一次计算C矩阵中4行×8列的小块。为什么是4×8AVX2一次能处理4个double所以4行恰好对上向量宽度8列则可以让C在寄存器里同时维护两个256位向量。4×8块需要维护8个向量寄存器用于C累加加上A的广播向量和B的加载向量总共大约12~14个ymm寄存器不影响后续循环展开。#include immintrin.h static inline void micro_kernel_4x8(int KC, const double *A_pack, // 4 x KC连续 const double *B_col, // 指向 B_pack 中 KC x NC 块的列偏移 int strideB, // B_pack 每行步长按double计 double *C, int ldC) { __m256d c[4][2]; for (int i 0; i 4; i) { c[i][0] _mm256_loadu_pd(C i * ldC); c[i][1] _mm256_loadu_pd(C i * ldC 4); } for (int k 0; k KC; k) { __m256d b0 _mm256_loadu_pd(B_col k * strideB); __m256d b1 _mm256_loadu_pd(B_col k * strideB 4); for (int i 0; i 4; i) { __m256d a _mm256_broadcast_sd(A_pack i * KC k); c[i][0] _mm256_fmadd_pd(a, b0, c[i][0]); c[i][1] _mm256_fmadd_pd(a, b1, c[i][1]); } } for (int i 0; i 4; i) { _mm256_storeu_pd(C i * ldC, c[i][0]); _mm256_storeu_pd(C i * ldC 4, c[i][1]); } }这段代码里首尾各有一轮load和store中间的KC次循环全部是FMA指令更新寄存器。_mm256_broadcast_sd把A的单个元素广播到4个向量lane里_mm256_fmadd_pd一次完成c a×b c。四个循环迭代共享b0和b1等于每轮k循环只要加载A的4个标量和B的2个向量就能完成4×832次乘加load和计算的比例降到了很健康的水平。4.3 为什么不再用更大或者更小的块微内核尺寸不是越大越好。如果NR从8增大到16每个i需要维护4个C向量4行就是16个向量寄存器还没算A和B的临时向量寄存器必然溢出到栈。溢出意味着额外的内存访问性能反而崩坏。反过来如果MR和NR太小比如2×4加载A和B的指令开销在总指令中占比过高FMA单元的利用率会下降。还有个经常被忽略的点微内核外层需要保证矩阵的边界能被整除。我在测试矩阵里选择M、N都为2048刚好被4×8整除。实际工程中遇到任意尺寸需要先计算边界块用标量版本补完或者对尾部做掩码处理。这部分代码虽然繁琐但逻辑上不复杂最难的是别把边界判断写进内层循环导致性能回落。4.4 引入微内核后的效果把5重循环的最内层换成4×8微内核后单线程场景下性能直接跳到108 GFLOPS。相比之前45.3 GFLOPS又翻了一倍多。顺带一提编译时一定要带上-marchnative否则GCC不会生成AVX2和FMA指令intrinsic会退化成函数调用或者编译报错。走到这一步单核性能已经比最开始的朴素版本提升了大约14倍。但离标题里的30倍还差得远剩余空间主要来自两个方面一是让微内核的访存更规整二是把八个核心全部用起来。5. 第四板斧Pack数据布局把微内核喂饱5.1 打包到底解决什么问题微内核的性能上来了新的瓶颈也随之出现每次微内核要从原始的A和B矩阵中取数据但原始数据的行跨度是K或N和微内核需要的连续小块并不一致。直接取数会导致大量非连续访问和地址计算开销更麻烦的是缓存利用率不稳定。解决办法是打包packing在进入微内核之前把当前需要计算的A块和B块拷贝到连续内存中拷贝后的布局完全按照微内核的读取模式设计。A块按照MC×KC连续存放B块按照KC×NC连续存放。虽然多了一次拷贝但这个拷贝是一次性的之后整个块可以在微内核中被反复读取代价远远小于收益。void pack_A(const double *A, int M, int K, int k0, double *A_pack) { for (int i 0; i MC; i) memcpy(A_pack i * KC, A (i0 i) * K k0, sizeof(double) * KC); } void pack_B(const double *B, int K, int N, int k0, int j0, double *B_pack) { for (int k 0; k KC; k) memcpy(B_pack k * NC, B (k0 k) * N j0, sizeof(double) * NC); }这里有个细节B_pack的每行宽度是NC而不是8。为什么要保留这么大的跨度因为一次pack整个KC×NC的B块后微内核在jj方向上每移动一个NR列只需要从B_pack中取不同列偏移即可所有微内核共享同一个打包好的B块。这个设计比每8列打包一次减少了几十倍的重复拷贝。5.2 对齐、跨页与TLB三个容易被忽略的坑Pack之后性能能再涨一截但前提是别踩几个暗坑。第一个坑是内存对齐。AVX2的对齐load指令_mm256_load_pd要求指针32字节对齐不对齐会直接崩溃。我前面代码用了_mm256_loadu_pd不对齐版本这是为了示例代码的普适性。实测证明如果你的数据能做到32字节对齐用对齐load会比不对齐load快几个百分点。分配内存时用aligned_alloc(32, ...)即可。第二个坑是跨页。当pack后的数组跨内存页边界时访问一次会触发两次TLB miss。解决方案是让每行的长度尽量避开页边界或者在pack时预留填充padding。比如B_pack每行NC×8字节如果NC4096刚好每行32KB连续两行正好跨在不同的4KB页上这其实问题不大因为一行有32KB的连续数据TLB能覆盖。更反直觉的是过度对齐到64字节甚至128字节可能出现缓存体冲突具体得看你CPU内部的缓存切片方式这类问题只能靠实际测量。第三个坑是数据初始化的页错误。大块内存首次访问时会有page fault开销我的做法是在计时前先遍历一遍两个pack缓冲区把它们在内存中落地否则首次运行的耗时会被路外干扰。5.3 Pack带来的实际性能变化加上Pack之后单线程性能从108 GFLOPS提升到158 GFLOPS。这接近50 GFLOPS的收益来自两方面访存局部性提高了地址计算也简化了——微内核里不再需要根据原始矩阵的lda/ldb计算偏移所有读取都是连续的固定步长。这一步看起来平平无奇但对稳定性很重要。Pack后的代码在反复运行中耗时波动明显变小因为数据布局固定之后缓存行为变得可预测了。6. 第五板斧多线程把八颗核都拉进来干活6.1 OpenMP并行不要在最内层加pragma到了这一步单核已经达到158 GFLOPS剩下的提升潜力在于并行化。很多人第一反应是给最内层循环加#pragma omp parallel for这是大忌微内核循环每次迭代才计算32个乘加线程同步和调度开销会彻底淹没收益。正确的并行位置在最外层的分块循环。我最开始尝试的是并行i0块循环#pragma omp parallel for schedule(static) for (int i0 0; i0 M; i0 MC) { for (int k0 0; k0 K; k0 KC) { pack_A(...); for (int j0 0; j0 N; j0 NC) { pack_B(...); for (...) { micro_kernel_4x8(...); } } } }每个线程完整地遍历k0和j0自己负责自己的MC块。好处是负载均衡极好——每个块的计算量完全一样线程之间没有依赖。坏处是B_pack会被每个线程重复打包浪费了一点带宽但相对计算量来说可以接受。实测在8线程下性能从158 GFLOPS提升到248 GFLOPS加速比大约1.57倍没有达到理论上的8倍这是正常的内存带宽和L3带宽成为多线程争抢的资源加上AVX2高负载下的频率回落能有这个提升已经不错了。6.2 并行粒度、负载均衡与超线程陷阱并行框架搭好之后还有三个影响最终性能的细节要处理。第一个是调度策略。我用schedule(static)而不是dynamic因为每个i0块的计算量完全均等static调度不需要运行时任务分配开销最小。只有当矩阵尺寸不是MC的整数倍、块与块之间负载不均时才考虑dynamic。第二个是超线程。我的测试机有8个物理核心、16个逻辑线程。实测开16线程比开8线程反而慢了一截。原因不难想两个逻辑线程共享同一个物理核心的FMA执行单元超线程不仅不增加浮点算力反而因为排队和缓存争抢拖慢整体速度。所以最后的结论是线程数设成物理核数。第三个是绑核。我在启动时通过环境变量把线程固定到物理核上尽量避免线程在核心间迁移导致缓存失效。这个操作在Windows和Linux下的做法不同但概念是一样的让每个线程稳定拥有自己的L1和L2。多线程这一步之后最终性能定格在248 GFLOPS对应标题里近30倍的说法——比最初的7.8 GFLOPS快了31.8倍。7. 踩坑与测试为什么跑不出30倍别急着怀疑优化无效7.1 性能测试要避开的坑如果读者复现我的流程很可能发现自己的提升倍数和我不一样。这不一定是哪里做错了多数情况下是测试方式有差异。有几点经验值得单独拿出来说。浮点性能的计算公式是FLOPS 2×M×N×K / 耗时。很多人忘记乘2算出来的数字直接打五折跟别人对比时就会被误导。计时时矩阵尺寸太小的坑也很严重。如果MNK256整个矩阵可能全部塞进L3缓存访存不再是瓶颈任何优化看起来都不明显。要让测试有意义矩阵至少要远超L3容量我建议从1024开始最好直接上2048甚至4096。还有频率问题。现代CPU频率波动极大同样一段代码这次跑0.07秒下次能跑到0.1秒。我的做法是每个优化版本都跑5次取中位数实测中这种方法能有效压掉异常值。如果条件允许在BIOS里关掉睿频或锁频数据会更稳定。7.2 精度问题FMA让误差变大了但不代表算错细心的读者应该已经注意到微内核版的计算结果和朴素版并不完全比特级一致。原因有两点累加顺序变了FMA把乘法结果先舍入再和累加项相加。这些差异在浮点数运算中完全正常通常只有最后几位不同。我在实战中遇到的更隐蔽的问题是用-O3 -ffast-math编译时编译器会激进地假设不存在NaN和Inf并重排浮点运算这可能导致原本应该出现的NaN被悄悄吞掉或者错误传播路径改变。对于数值计算项目我建议不要贪图那点性能开-ffast-math除非你完全接受它带来的语义变化。相比之下-marchnative和-O3已经足够。7.3 完整性能汇总把整个优化链条的每一步数据汇总成一张表看趋势一目了然优化阶段性能GFLOPS相对朴素提升朴素ijk三重循环7.81×ikj循环重排19.62.5×缓存分块45.35.8×4×8微内核AVX2 FMA10813.8×数据Pack布局15820.3×8线程并行24831.8×需要强调的是这些绝对值只对这台测试机有参考意义。换一颗CPU、换一套内存每个阶段的收益比例都会变化但优化顺序和思路是通用的先解决访存再做向量化最后上并行。8. 优化到30倍之后的一点感想这套优化流程走完我最大的体会是性能优化没有玄学只有找到瓶颈的决心。很多时候性能上不去不是因为你不知道某个高级指令而是你根本没有把问题拆成访存模型、计算模型、线程模型三个层次逐一检查。如果你的场景是生产环境直接用成熟的开源BLAS库会比手写更快更稳。我写这个项目的过程本身就是学习而不是为了造一个替代品。通过亲手实现一遍你会理解为什么那些BLAS库内部有那么多奇怪的循环嵌套和打包步骤理解为什么它们在不同CPU上需要调优参数也才能在碰到库解决不了的边缘情况时有能力去修改和定制。这个优化路径继续走下去还有很多空间比如针对AVX-512做32位浮点版本的微内核、用SVE适配ARM平台、把微内核里的循环再做展开和双发射调度、引入编译期自动调优等。我后面的计划是把手写的微内核重构一下塞进一个简单的自动调优框架里让它能在不同机器上自己搜出最优的分块参数和微内核尺寸然后对比一下和开源BLAS还差多少。那又会是一篇值得写的东西。