3. 循环分块提升矩阵转置的缓存局部性

问题定位：为什么朴素转置会抖动

矩阵转置 B[j][i] = A[i][j] 的访存模式是“读一行、写一列”。以行主序存储的 N×N 双精度矩阵为例，A 的读取是连续的，但 B 的写入每次跨越一整行（N×8 字节）。当 N 较大时，同一缓存行内的相邻元素在时间上相隔很远才被写入，导致写缓存行在被逐出前无法被充分利用，形成缓存抖动：L1/L2 缺失率随 N 增大而快速上升。

循环分块（blocking / tiling）的思路是把大转置拆成 BS×BS 的小块，使一个块内的读写都落在少量缓存行上，从而让每个被载入的缓存行在逐出前被充分复用。分块尺寸 BS 是关键参数：太小则循环开销与边界处理占比升高，太大则工作集超出 L1/L2 容量，缺失率重新上升。

分块转置的实现

下面给出可直接编译运行的 C 程序。它包含朴素转置与分块转置两个版本，分块版本对 BS 做参数化，并支持通过命令行扫描多个分块尺寸。矩阵按行主序存储，N 不要求是 BS 的整数倍，边界块用 min 截断处理。

// transpose_block.c
// 编译: gcc -O2 -march=native -o transpose_block transpose_block.c
// 运行: ./transpose_block 4096 8 16 32 64 128
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>

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

// 朴素转置：B[j*N+i] = A[i*N+j]
void transpose_naive(const double *A, double *B, int N) {
    for (int i = 0; i < N; ++i)
        for (int j = 0; j < N; ++j)
            B[(size_t)j * N + i] = A[(size_t)i * N + j];
}

// 分块转置：外层按 BS 分块，块内做转置
void transpose_blocked(const double *A, double *B, int N, int BS) {
    for (int ii = 0; ii < N; ii += BS) {
        int i_end = ii + BS < N ? ii + BS : N;
        for (int jj = 0; jj < N; jj += BS) {
            int j_end = jj + BS < N ? jj + BS : N;
            for (int i = ii; i < i_end; ++i) {
                for (int j = jj; j < j_end; ++j) {
                    B[(size_t)j * N + i] = A[(size_t)i * N + j];
                }
            }
        }
    }
}

// 校验：与朴素结果逐元素比较，返回不一致元素个数
long verify(const double *A, const double *B, int N) {
    long bad = 0;
    for (int i = 0; i < N; ++i)
        for (int j = 0; j < N; ++j)
            if (B[(size_t)j * N + i] != A[(size_t)i * N + j]) ++bad;
    return bad;
}

int main(int argc, char **argv) {
    if (argc < 2) {
        fprintf(stderr, "用法: %s N [BS1 BS2 ...]\n", argv[0]);
        return 1;
    }
    int N = atoi(argv[1]);
    size_t elems = (size_t)N * N;
    double *A = malloc(elems * sizeof(double));
    double *B = malloc(elems * sizeof(double));
    if (!A || !B) { fprintf(stderr, "内存分配失败\n"); return 1; }

    for (size_t k = 0; k < elems; ++k) A[k] = (double)(k % 1000);

    // 基线：朴素转置
    double t0 = now_sec();
    transpose_naive(A, B, N);
    double t_naive = now_sec() - t0;
    printf("N=%d 朴素转置: %.4f s\n", N, t_naive);

    // 分块扫描
    for (int a = 2; a < argc; ++a) {
        int BS = atoi(argv[a]);
        if (BS <= 0) continue;
        memset(B, 0, elems * sizeof(double));
        t0 = now_sec();
        transpose_blocked(A, B, N, BS);
        double t = now_sec() - t0;
        long bad = verify(A, B, N);
        printf("BS=%-4d 时间=%.4f s  加速比=%.2fx  校验不一致=%ld\n",
               BS, t, t_naive / t, bad);
    }

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

测量缓存缺失率

执行时间只能反映综合效果，要确认分块确实改善了缓存局部性，需要看 L1/L2 缺失率。用 perf stat 采集缓存事件，对每个分块尺寸分别运行一次：

# 采集 L1 数据缓存缺失与 L2 缺失（事件名随 CPU 架构不同略有差异）
perf stat -e L1-dcache-loads,L1-dcache-load-misses,\
LLC-loads,LLC-load-misses \
./transpose_block 4096 8 16 32 64 128

判读要点：L1-dcache-load-misses / L1-dcache-loads 得到 L1 缺失率，LLC-load-misses / LLC-loads 得到末级缓存缺失率。分块有效时，应看到缺失率随 BS 从很小值增大而明显下降，在某个 BS 附近达到低谷，之后因工作集超出缓存容量而回升。若缺失率几乎不随 BS 变化，通常说明编译器已自动优化，或测量的事件与目标缓存层级不匹配。

参数扫描与结果判读

以 N=4096（双精度，单矩阵约 128 MiB）为例，扫描 BS ∈ {8, 16, 32, 64, 128}。下表是示例性的预期趋势，用于说明如何判读，实际数值需在目标平台上实测：

确定较优参数的实用做法：先看缺失率曲线的低谷区间，再在该区间内用执行时间做二次筛选。若两者低谷不一致，以执行时间为准，因为缺失率只是中间指标。对 N 不是 BS 整数倍的情况，还要确认边界块处理没有引入明显额外开销——可对比 N=4096 与 N=4100 在相同 BS 下的时间差异。

可复用的扫描记录表

把每次扫描的关键数据记下来，便于横向比较和复现：

完成标准：所有分块尺寸的校验不一致数均为 0（结果正确）；缺失率曲线出现可解释的先降后升；在目标平台上给出一个较优 BS 并说明依据（缺失率低谷与时间低谷的折中）。