缓存分块（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字节）：

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

3.2 实际考虑

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

4. 分块的局限

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

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

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

5.1 Intel ICC/ICX

#pragma omp parallel for
for (int i = 0; i < N; i++)
  #pragma unroll
  for (int j = 0; j < N; j++)
    // 编译器自动分块

5.2 LLVM-Polly^[2]

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要匹配缓存容量：
2. 实际B值通常取64或128（考虑SIMD对齐）
3. 多层分块（L1/L2/L3）能进一步提升性能
4. 现代编译器能自动分块，但手工调优仍有价值

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

参考

1. Intel Advisor Cookbook. Optimize Memory Access Patterns using Loop Interchange and Cache Blocking. 1.68x speedup with cache blocking. ↩

2. LLVM Dev Meeting. Cache-aware Scheduling and Performance Modeling with LLVM-Polly. Tile size selection algorithm. ↩