缓存与局部性：为什么换个循环顺序会快十倍

有一次，我写了一个矩阵乘法，自认为复杂度是 $O(n^3)$，应该已经「最优」了。

然后我手滑把 i-j-k 三层循环换成了 i-k-j，运行时间从 2.3 秒掉到 0.21 秒。

同样的运算量，同样的编译器，同样的 CPU。差别只有一个：内存访问的顺序。

一个极简的实验

为了公平，先准备一个朴素版本：

#define N 1024

void matmul_ijk(const double *A, const double *B, double *C) {
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            double sum = 0.0;
            for (int k = 0; k < N; k++) {
                sum += A[i * N + k] * B[k * N + j];
            }
            C[i * N + j] = sum;
        }
    }
}

再写一个只是交换了最内层循环的版本：

void matmul_ikj(const double *A, const double *B, double *C) {
    for (int i = 0; i < N; i++) {
        for (int k = 0; k < N; k++) {
            const double a = A[i * N + k];
            for (int j = 0; j < N; j++) {
                C[i * N + j] += a * B[k * N + j];
            }
        }
    }
}

在编译时记得关闭那些会改变结果的激进优化，否则编译器可能把整个循环优化掉：

gcc -O2 -fno-tree-vectorize -o bench bench.c
perf stat -e cycles,instructions,cache-references,cache-misses ./bench

在这台机器上，ijk 版本大约需要 2.1 秒，ikj 版本大约 0.25 秒。

缓存不是一层，而是一座塔

CPU 和主存之间隔着好几层缓存：

层级	典型容量	访问延迟（近似）	由谁管理
寄存器	几百字节	< 1 个周期	编译器
L1	32–64 KiB	约 4 个周期	硬件
L2	256 KiB–2 MiB	约 12–20 个周期	硬件
L3	8–64 MiB	约 40–80 个周期	硬件
主存	几 GiB 起	约 200–400 个周期	操作系统

这些数字因 CPU 而异，但数量级的关系很稳定：L1 与主存之间差了两个数量级。

缓存以 cache line 为单位工作。x86-64 上一条 cache line 通常是 64 字节，也就是 8 个 double，或 16 个 int。

当你读取 A[i] 时，硬件并不是只取那 8 个字节，而是把整条 cache line 搬进 L1。这意味着：如果你接下来访问相邻元素，它们几乎不要钱。

数组在内存里是怎么躺着的

C 的二维数组按行优先（row-major）存储。一个 $4 \times 4$ 的矩阵：

A[0][0] A[0][1] A[0][2] A[0][3]
A[1][0] A[1][1] A[1][2] A[1][3]
A[2][0] A[2][1] A[2][2] A[2][3]
A[3][0] A[3][1] A[3][2] A[3][3]

在内存里其实是：

[0][0] [0][1] [0][2] [0][3] [1][0] [1][1] ... [3][3]

一行是连续的，行与行之间也是连续的。

回到 ijk 版本。最内层是 k，访问的是 A[i * N + k] 和 B[k * N + j]：

• A 沿着连续的行移动，局部性很好；
• B 的地址每次跳跃 N 个元素，也就是 8 KiB 的跨度。

对于 1024×1024 的 double 矩阵，8 KiB 意味着每次访问几乎都落在不同的 cache line 上。而 B 的总大小是 8 MiB，远超 L2，甚至可能超过共享 L3 的可用部分。

于是 ijk 版本在不停地把 B 的 cache line 搬进来、赶出去，再搬进来。

ikj 版本则相反。最内层是 j，B[k * N + j] 和 C[i * N + j] 都是连续的：

• 读 B 时一整条 cache line 被用满；
• 写 C 时也沿着连续地址推进；
• A[i * N + k] 在外层循环里只读一次，几乎不影响。

数据的复用模式改变了，缓存命中率也随之改变。

局部性：性能的隐形货币

计算机体系结构里常说的两种局部性，在这段代码里同时出现：

空间局部性：访问了 B[k][j]，就很可能很快访问 B[k][j+1]。因为它们在同一个 cache line 里。

时间局部性：C[i][j] 在 k 循环中被反复读写。如果它能留在缓存里，就不必每次都去主存。

ijk 版本对 C 的使用是「算完一个格子就再也不碰」，时间局部性很弱；对 B 的使用是「每个元素都跨行访问」，空间局部性也很弱。

而 ikj 版本把 C 的一整行留在缓存里反复累加，同时顺序扫描 B 的一行。两种局部性都被照顾到了。

分块：当矩阵大到装不下

如果矩阵更大，比如 8192×8192，即使 ikj 版本也会开始吃力，因为一行 B 就有 64 KiB，已经超过 L1。

这时可以用分块（blocking / tiling）：把矩阵切成能塞进缓存的小块，在块内完成所有乘法。

#define BS 64

void matmul_blocked(const double *A, const double *B, double *C) {
    for (int ii = 0; ii < N; ii += BS) {
        for (int kk = 0; kk < N; kk += BS) {
            for (int jj = 0; jj < N; jj += BS) {
                for (int i = ii; i < ii + BS; i++) {
                    for (int k = kk; k < kk + BS; k++) {
                        const double a = A[i * N + k];
                        for (int j = jj; j < jj + BS; j++) {
                            C[i * N + j] += a * B[k * N + j];
                        }
                    }
                }
            }
        }
    }
}

块大小的选择不是玄学。一个常用的经验是让三个块同时放进 L2：

$$3 \cdot B^2 \cdot \text{sizeof(double)} \lesssim \text{L2 size}$$

如果 L2 是 512 KiB，每个 double 占 8 字节，那么 $B \approx \sqrt{512 \cdot 1024 / 24} \approx 148$。实际中还要考虑缓存关联度和预取，64 或 128 通常是安全的起点。

用 perf 看见缓存

不要凭感觉优化。Linux 的 perf 可以直接统计缓存事件：

perf stat -e cache-references,cache-misses,\
L1-dcache-load-misses,LLC-load-misses ./bench

三个数字值得关注：

• cache-references 与 cache-misses 的比值，反映整体命中率；
• L1-dcache-load-misses：L1 数据缓存未命中次数；
• LLC-load-misses：最后一级缓存未命中，通常意味着要去主存。

对比两个版本的 LLC-load-misses，你会看到数量级的差异。这比「感觉快了」有说服力得多。

别把它优化成玄学

缓存优化很容易走火入魔。几个现实的提醒：

• 先测量，再优化。 没有 profile 的优化是在猜。
• 编译器可能已经做了。 -O3、自动向量化、循环交换，现代编译器比我们想的聪明。先看生成的汇编，再决定手写什么。
• 数据规模决定策略。 小矩阵本来就在缓存里，分块反而增加开销。
• 可读性也是性能。 一个快 5% 但没人看得懂的循环，长期成本可能更高。

数组不是数学对象

在学校里，矩阵是一个抽象的二维表格；在机器里，它是一维的字节序列，被 cache line 切成一段一段。

理解这一点之后，很多「莫名其妙」的性能现象就有了形状：为什么按行遍历比按列快，为什么链表不如数组缓存友好，为什么某些看起来更复杂的算法反而更快。

算法决定你要走多少步，内存布局决定每一步有多贵。