3. 循环分块（Loop Tiling）

任务目标

本案例针对矩阵乘法 C = A × B，通过循环分块（Loop Tiling）将计算划分为适合 L1 缓存大小的子块，减少缓存行替换，提升数据局部性。你需要对比朴素三重循环与分块实现的缓存行为差异，并验证分块后执行时间下降。

环境与依赖

• 编译器：gcc（支持 -O2 优化）
• 可选：perf stat 或 time 命令用于测量
• 无第三方库依赖

核心实现

以下代码包含朴素版本与分块版本，矩阵大小 N=512，分块大小 BLOCK=32（适配典型 32KB L1 缓存，每个 double 8 字节，32×32×8×3 ≈ 24KB）。

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

#define N 512
#define BLOCK 32

void matmul_naive(double *A, 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_tiled(double *A, double *B, double *C) {
    for (int i0 = 0; i0 < N; i0 += BLOCK)
        for (int j0 = 0; j0 < N; j0 += BLOCK)
            for (int k0 = 0; k0 < N; k0 += BLOCK)
                for (int i = i0; i < i0 + BLOCK; i++)
                    for (int j = j0; j < j0 + BLOCK; j++) {
                        double sum = C[i * N + j];
                        for (int k = k0; k < k0 + BLOCK; k++)
                            sum += A[i * N + k] * B[k * N + j];
                        C[i * N + j] = sum;
                    }
}

void init(double *M) {
    for (int i = 0; i < N * N; i++)
        M[i] = (double)(rand() % 100) / 10.0;
}

double now_ms() {
    struct timespec ts;
    clock_gettime(CLOCK_MONOTONIC, &ts);
    return ts.tv_sec * 1000.0 + ts.tv_nsec / 1e6;
}

int main() {
    double *A = malloc(N * N * sizeof(double));
    double *B = malloc(N * N * sizeof(double));
    double *C1 = calloc(N * N, sizeof(double));
    double *C2 = calloc(N * N, sizeof(double));
    srand(42);
    init(A); init(B);

    double t0 = now_ms();
    matmul_naive(A, B, C1);
    double t1 = now_ms();
    matmul_tiled(A, B, C2);
    double t2 = now_ms();

    // 正确性校验
    double max_diff = 0.0;
    for (int i = 0; i < N * N; i++) {
        double d = C1[i] - C2[i];
        if (d < 0) d = -d;
        if (d > max_diff) max_diff = d;
    }

    printf("Naive: %.2f ms\n", t1 - t0);
    printf("Tiled: %.2f ms\n", t2 - t1);
    printf("Max diff: %e\n", max_diff);

    free(A); free(B); free(C1); free(C2);
    return 0;
}

编译与运行

gcc -O2 -o matmul matmul.c -lm
./matmul

验证方法

• 正确性：程序输出 Max diff 应小于 1e-10（浮点累加顺序差异导致微小误差），表明分块结果与朴素版本一致。
• 性能对比：观察 Tiled 耗时是否明显低于 Naive。在 N=512、BLOCK=32 时，分块版本通常快 2～4 倍（具体取决于 CPU 缓存大小）。
• 缓存行为（可选）：运行 perf stat -e cache-misses ./matmul，对比两次调用的 cache-misses 计数，分块版本应显著更低。

关键说明

分块的核心在于将三重循环的 k 维度也切块，使 A 的子块（32×32）和 B 的子块（32×32）能同时驻留于 L1 缓存，避免反复从内存加载。调整 BLOCK 值可观察性能变化——过小无法充分利用缓存行，过大则导致缓存溢出。若你的 CPU L1 缓存不是 32KB，可按 BLOCK = sqrt(L1_cache_size / 3 / 8) 估算。