2. 循环展开与寄存器分块对矩阵乘法的性能影响

任务边界与基线程序

本节固定矩阵规模为 N=1024，数据类型为 double，采用行主序存储，计算 C = A × B。所有版本共用同一份随机初始化数据与同一套计时框架，只改变内层循环的组织方式，从而把性能差异归因到循环展开与寄存器分块本身。

朴素三重循环按 i-j-k 顺序遍历，内层对 B 的访问步长为 N，跨行跳跃导致缓存行利用率低；同时 C[i][j] 在每次 k 迭代中被反复读写，编译器难以把它长期留在寄存器中。基线程序如下：

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define N 1024

static double A[N][N], B[N][N], C[N][N];

static void init(void) {
    for (int i = 0; i < N; i++)
        for (int j = 0; j < N; j++) {
            A[i][j] = (double)((i * 31 + j * 17) % 100) / 100.0;
            B[i][j] = (double)((i * 13 + j * 29) % 100) / 100.0;
            C[i][j] = 0.0;
        }
}

static double now_sec(void) {
    struct timespec ts;
    clock_gettime(CLOCK_MONOTONIC, &ts);
    return ts.tv_sec + ts.tv_nsec * 1e-9;
}

int main(void) {
    init();
    double t0 = now_sec();
    for (int i = 0; i < N; i++)
        for (int j = 0; j < N; j++) {
            double s = 0.0;
            for (int k = 0; k < N; k++)
                s += A[i][k] * B[k][j];
            C[i][j] = s;
        }
    double t1 = now_sec();
    printf("naive: %.3f s, checksum=%.6f\n", t1 - t0, C[0][0] + C[N-1][N-1]);
    return 0;
}

编译时统一使用 gcc -O2 -march=native -fno-tree-vectorize，关闭自动向量化以隔离循环展开与寄存器分块这两个变量的影响。校验和用于确认不同版本计算结果一致，避免优化引入错误。

循环展开：减少分支与提升指令级并行

对最内层 k 循环做 4 路展开，一次迭代处理 4 个乘加，循环控制开销降为原来的 1/4，同时给乱序执行引擎提供更多独立指令。展开后仍需保证尾部剩余元素被正确处理：

for (int i = 0; i < N; i++)
    for (int j = 0; j < N; j++) {
        double s0 = 0, s1 = 0, s2 = 0, s3 = 0;
        int k = 0;
        for (; k + 3 < N; k += 4) {
            s0 += A[i][k]     * B[k][j];
            s1 += A[i][k + 1] * B[k + 1][j];
            s2 += A[i][k + 2] * B[k + 2][j];
            s3 += A[i][k + 3] * B[k + 3][j];
        }
        for (; k < N; k++) s0 += A[i][k] * B[k][j];
        C[i][j] = s0 + s1 + s2 + s3;
    }

四个累加器 s0..s3 相互独立，打破了单一累加链上的浮点加法依赖，使乘加流水线可以重叠执行。但注意：B[k][j] 的访问模式没有改变，缓存缺失问题依旧存在，因此展开主要改善 IPC，对缓存缺失率影响有限。

寄存器分块：让 C 的多个元素常驻寄存器

寄存器分块（register blocking）在 i-j 平面上取 MR×NR 的小块，例如 4×4，把 16 个 C 累加值全部放在局部变量中，内层 k 每步只加载 4 个 A 元素和 4 个 B 元素，完成 16 次乘加，算术强度从 1 次乘加/2 次加载提升到 16 次乘加/8 次加载：

#define MR 4
#define NR 4

for (int i = 0; i < N; i += MR)
    for (int j = 0; j < N; j += NR) {
        double c[MR][NR] = {{0}};
        for (int k = 0; k < N; k++) {
            double a[MR], b[NR];
            for (int m = 0; m < MR; m++) a[m] = A[i + m][k];
            for (int n = 0; n < NR; n++) b[n] = B[k][j + n];
            for (int m = 0; m < MR; m++)
                for (int n = 0; n < NR; n++)
                    c[m][n] += a[m] * b[n];
        }
        for (int m = 0; m < MR; m++)
            for (int n = 0; n < NR; n++)
                C[i + m][j + n] = c[m][n];
    }

由于 N=1024 能被 4 整除，此处省略边界处理；实际工程中应对不能整除的规模补上尾部循环。该版本把 B 的访问从“每元素一次跨行跳跃”变为“连续 4 个元素一次加载”，显著提升缓存行利用率。

计时与 perf 统计的采集方式

每个版本单独编译为可执行文件，先运行一次预热，再连续运行 3 次取中位数作为耗时。perf 统计使用同一命令模板，只替换可执行文件名：

perf stat -e cycles,instructions,cache-references,cache-misses,\
L1-dcache-loads,L1-dcache-load-misses ./matmul_naive

IPC 由 instructions / cycles 直接算出；缓存缺失率由 cache-misses / cache-references 得到。为保证可比性，三次测量都在同一台机器、同一 CPU 频率策略下完成，并记录 perf stat 输出的 seconds time elapsed 与计时程序结果互相印证。

结果判读与归因方法

把三个版本放在一起看，指标变化应呈现清晰的分工：V2 相对 V1 的 IPC 上升而缓存缺失率基本持平，说明收益来自指令级并行与分支开销下降；V3 相对 V2 的缓存缺失率明显下降、算术强度提升，说明收益来自数据复用。若 V3 的 IPC 反而低于 V2，常见原因是 16 个累加器加 8 个临时变量超出可用寄存器数量，编译器发生寄存器溢出（spill），此时应把分块缩小到 4×2 或 2×4 再测。

判读时按以下顺序核对，避免把测量噪声当成优化效果：

• 先确认三个版本的校验和一致，排除计算错误导致的“假加速”。
• 对比 perf stat 的 seconds time elapsed 与程序内计时，两者偏差应在 5% 以内。
• 若 cache-misses 计数接近 0 或异常巨大，检查是否缺少 perf 权限或事件不被硬件支持。
• IPC 提升但耗时不变，说明瓶颈不在前端，应转向内存带宽或频率策略排查。