矩阵乘法是科学计算的核心但 naive 实现性能惨不忍睹。问题出在缓存——三个大矩阵来回折腾L1缓存根本装不下。缓存分块Cache Blocking/Tiling通过把大矩阵切成小块让数据在缓存里多待一会儿性能能提升几倍。1. 问题传统矩阵乘法的缓存噩梦标准的三层循环矩阵乘法for(i0;iN;i)for(j0;jN;j){r0;for(k0;kN;k)ry[i][k]*z[k][j];// z按列访问x[i][j]r;}问题在哪空间局部性差z[k][j]按列访问但C数组是行主序。z[0][j]和z[1][j]在内存中相隔N个元素大概率不在同一缓存行。时间局部性差计算x[i][j]需要y的第i行和z的第j列。下一个元素x[i][j1]又要重新加载y的第i行——虽然刚刚用过但可能被踢出缓存了。总访问量假设三个N×N矩阵计算量是2N³次操作但内存访问量也是O(N³)级别。如果N1000缓存装不下每次都要从内存读性能暴跌。2. 分块优化把大矩阵切成小块核心思想把矩阵分成B×B的小块确保三个块能同时驻留缓存。for(jj0;jjN;jjB)// 分块列循环for(kk0;kkN;kkB)// 分块行循环for(i0;iN;i)for(jjj;jmin(jjB,N);j){r0;for(kkk;kmin(kkB,N);k)ry[i][k]*z[k][j];// 块内访问x[i][j]r;}关键变化最内层循环只在B×B的块内操作如果3B² ≤ 缓存容量三个块都能驻留块内数据复用减少内存访问3. 分块因子的选择分块因子B不是越大越好要匹配缓存容量。3.1 理论计算假设L1缓存32KBfloat类型4字节3 B 2 × 4 ≤ 32768 3B^2 \times 4 \leq 327683B2×4≤32768B ≤ 32768 / 12 ≈ 52 B \leq \sqrt{32768 / 12} \approx 52B≤32768/12​≈52所以B≈52取整64方便SIMD对齐。3.2 实际考虑因素影响建议缓存关联性8路组相联需避免Bank冲突B取2的幂次寄存器压力B太小循环展开效率低B≥16SIMD宽度AVX-512一次算16个floatB是16的倍数TLB容量B太大可能跨页B≤512Intel Advisor的实测建议1对于矩阵乘法B64是甜点区。4. 分块的局限小块开销当N很大但B固定时分块引入的循环开销可以忽略。但如果N本身很小如N100分块反而增加开销。不规则矩阵非方阵或稀疏矩阵分块效果打折扣。5. 现代编译器的自动分块5.1 Intel ICC/ICX#pragmaomp parallelforfor(inti0;iN;i)#pragmaunrollfor(intj0;jN;j)// 编译器自动分块5.2 LLVM-Polly2Polly是LLVM的多面体优化框架能自动进行循环分块clang-O3-mllvm-polly-mllvm-polly-tile ./matmul.cPolly的tile size选择算法考虑缓存大小缓存行大小循环迭代次数5.3 自动调优Auto-tuningATLAS和OpenBLAS采用empirical tuning编译多个版本的kernel不同B值在目标机器上实测选择最快的版本这比理论计算更准确因为考虑了缓存替换策略TLB行为预取器影响6. 总结缓存分块是矩阵运算优化的核心技术优化效果复杂度Loop Interchange解决空间局部性低Cache Blocking解决时间局部性中SIMD向量化提升单周期算力中多层分块利用整个缓存层次高关键认知分块因子B要匹配缓存容量3 B 2 ≤ C L 1 3B^2 \leq C_{L1}3B2≤CL1​实际B值通常取64或128考虑SIMD对齐多层分块L1/L2/L3能进一步提升性能现代编译器能自动分块但手工调优仍有价值理解分块就能理解为什么OpenBLAS/GotoBLAS比naive实现快10倍以上。参考Intel Advisor Cookbook. Optimize Memory Access Patterns using Loop Interchange and Cache Blocking. 1.68x speedup with cache blocking. ↩︎LLVM Dev Meeting. Cache-aware Scheduling and Performance Modeling with LLVM-Polly. Tile size selection algorithm. ↩︎