缓存分块（Cache Blocking）：矩阵乘法的救命稻草

矩阵乘法是科学计算的核心，但 naive 实现性能惨不忍睹。问题出在缓存——三个大矩阵来回折腾，L1缓存根本装不下。缓存分块（Cache Blocking/Tiling）通过把大矩阵切成小块，让数据在缓存里多待一会儿，性能能提升几倍。

1. 问题：传统矩阵乘法的缓存噩梦

标准的三层循环矩阵乘法：

for(i=0;i<N;i++)for(j=0;j<N;j++){r=0;for(k=0;k<N;k++)r+=y[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][j+1]又要重新加载y的第i行——虽然刚刚用过，但可能被踢出缓存了。

总访问量：假设三个N×N矩阵，计算量是2N³次操作，但内存访问量也是O(N³)级别。如果N=1000，缓存装不下，每次都要从内存读，性能暴跌。

2. 分块优化：把大矩阵切成小块

核心思想：把矩阵分成B×B的小块，确保三个块能同时驻留缓存。

for(jj=0;jj<N;jj+=B)// 分块列循环for(kk=0;kk<N;kk+=B)// 分块行循环for(i=0;i<N;i++)for(j=jj;j<min(jj+B,N);j++){r=0;for(k=kk;k<min(kk+B,N);k++)r+=y[i][k]*z[k][j];// 块内访问x[i][j]+=r;}

关键变化：

• 最内层循环只在B×B的块内操作
• 如果3B² ≤ 缓存容量，三个块都能驻留
• 块内数据复用，减少内存访问

3. 分块因子的选择

分块因子B不是越大越好，要匹配缓存容量。

3.1 理论计算

假设L1缓存32KB，float类型（4字节）：

3 B 2 × 4 ≤ 32768 3B^2 \times 4 \leq 327683B2×4≤32768

B ≤ 32768 / 12 ≈ 52 B \leq \sqrt{32768 / 12} \approx 52B≤32768/12​≈52

所以B≈52，取整64（方便SIMD对齐）。

3.2 实际考虑

Intel Advisor的实测建议¹：对于矩阵乘法，B=64是甜点区。

4. 分块的局限

小块开销：当N很大但B固定时，分块引入的循环开销可以忽略。但如果N本身很小（如N<100），分块反而增加开销。

不规则矩阵：非方阵或稀疏矩阵，分块效果打折扣。

5. 现代编译器的自动分块

5.1 Intel ICC/ICX

#pragmaomp parallelforfor(inti=0;i<N;i++)#pragmaunrollfor(intj=0;j<N;j++)// 编译器自动分块

5.2 LLVM-Polly²

Polly是LLVM的多面体优化框架，能自动进行循环分块：

clang-O3-mllvm-polly-mllvm-polly-tile ./matmul.c

Polly的tile size选择算法考虑：

• 缓存大小
• 缓存行大小
• 循环迭代次数

5.3 自动调优（Auto-tuning）

ATLAS和OpenBLAS采用empirical tuning：

1. 编译多个版本的kernel，不同B值
2. 在目标机器上实测
3. 选择最快的版本

这比理论计算更准确，因为考虑了：

• 缓存替换策略
• TLB行为
• 预取器影响

6. 总结

缓存分块是矩阵运算优化的核心技术：

关键认知：

1. 分块因子B要匹配缓存容量：3 B 2 ≤ C L 1 3B^2 \leq C_{L1}3B2≤CL1​
2. 实际B值通常取64或128（考虑SIMD对齐）
3. 多层分块（L1/L2/L3）能进一步提升性能
4. 现代编译器能自动分块，但手工调优仍有价值

理解分块，就能理解为什么OpenBLAS/GotoBLAS比naive实现快10倍以上。

参考