3. CSR 矩阵向量乘法的缓存友好实现

CSR 的存储结构与 SpMV 的访问模式

CSR（Compressed Sparse Row）用三个数组描述矩阵：values 按行连续存放非零元数值，col_idx 存放对应的列号，row_ptr 长度为 m+1，row_ptr[i] 到 row_ptr[i+1]-1 就是第 i 行的非零元区间。SpMV 计算 y = A·x 时，第 i 行的结果是 y[i] = Σ values[k] * x[col_idx[k]]，k 遍历该行的非零元。

性能瓶颈不在乘法本身，而在 x[col_idx[k]] 的间接访存：col_idx 和 values 是顺序读取的（缓存友好），但 x 的访问位置由列号决定，是随机的。因此优化的核心是让 x 的访问尽量命中缓存，并让内层循环保持简单、可向量化。

缓存友好的实现要点

下面给出一个可直接编译运行的 C 实现，包含 CSR 构建、两种循环顺序的 SpMV、以及用稠密矩阵对照的正确性验证。示例矩阵为 5×10、共 11 个非零元，与上图数据一致。

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

typedef struct {
    int m, n, nnz;
    int    *row_ptr;   /* 长度 m+1 */
    int    *col_idx;   /* 长度 nnz */
    double *values;    /* 长度 nnz */
} CSR;

/* 由稠密矩阵构建 CSR（行优先扫描，天然有序） */
CSR csr_from_dense(const double *A, int m, int n) {
    CSR a; a.m = m; a.n = n;
    a.row_ptr = malloc((m + 1) * sizeof(int));
    a.col_idx = malloc((size_t)m * n * sizeof(int));
    a.values  = malloc((size_t)m * n * sizeof(double));
    int k = 0;
    for (int i = 0; i < m; ++i) {
        a.row_ptr[i] = k;
        for (int j = 0; j < n; ++j) {
            double v = A[(size_t)i * n + j];
            if (v != 0.0) { a.col_idx[k] = j; a.values[k] = v; ++k; }
        }
    }
    a.row_ptr[m] = k;
    a.nnz = k;
    return a;
}

/* 行优先：外层行、内层非零元，y 只写一次 */
void spmv_row_major(const CSR *a, const double *x, double *y) {
    for (int i = 0; i < a->m; ++i) {
        double sum = 0.0;
        int end = a->row_ptr[i + 1];
        for (int k = a->row_ptr[i]; k < end; ++k)
            sum += a->values[k] * x[a->col_idx[k]];
        y[i] = sum;
    }
}

/* 列优先：按非零元顺序遍历，y 随机累加（单线程下结果相同） */
void spmv_col_major(const CSR *a, const double *x, double *y) {
    memset(y, 0, a->m * sizeof(double));
    for (int i = 0; i < a->m; ++i)
        for (int k = a->row_ptr[i]; k < a->row_ptr[i + 1]; ++k)
            y[i] += a->values[k] * x[a->col_idx[k]];
}

int main(void) {
    const int m = 5, n = 10;
    double A[5][10] = {
        {2.0, 0,   0,  -1.0, 0,   0,   0,   0,   0,   0  },
        {0,   4.0, 0,   0,   1.5, 0,   0,  -2.0, 0,   0  },
        {0,   0,   3.0, 0,   0,   0.5, 1.0, 0,   0,   0  },
        {0,   0,   0,   0,   0,   0,   0,   0,   0,   0  },
        {0,  -1.0, 0,   0,   0,   0,   0,   0,   2.0, 1.0}
    };
    double x[10] = {1,2,3,4,5,6,7,8,9,10};
    double y1[5], y2[5], yd[5];

    CSR a = csr_from_dense(&A[0][0], m, n);
    printf("nnz = %d\n", a.nnz);

    spmv_row_major(&a, x, y1);
    spmv_col_major(&a, x, y2);

    /* 稠密对照 */
    for (int i = 0; i < m; ++i) {
        yd[i] = 0;
        for (int j = 0; j < n; ++j) yd[i] += A[i][j] * x[j];
    }

    int ok = 1;
    for (int i = 0; i < m; ++i) {
        printf("y[%d] row=%.3f col=%.3f dense=%.3f\n", i, y1[i], y2[i], yd[i]);
        if (y1[i] != yd[i] || y2[i] != yd[i]) ok = 0;
    }
    printf("verify: %s\n", ok ? "PASS" : "FAIL");

    free(a.row_ptr); free(a.col_idx); free(a.values);
    return 0;
}

预期输出（由给定输入直接推出）：nnz = 11，y[0] = 2·1 + (−1)·4 = −2，y[1] = 4·2 + 1.5·5 + (−2)·8 = −0.5，y[2] = 3·3 + 0.5·6 + 1·7 = 19，y[3] = 0，y[4] = (−1)·2 + 2·9 + 1·10 = 26，最后打印 verify: PASS。

循环顺序的性能差异与判读

两种写法在单线程下数值结果一致，但访存行为不同：行优先版本对 y 只写一次，内层循环连续读 values 与 col_idx，只有 x 是间接访问；列优先版本让 x 顺序读，但 y 变成随机累加，且在多线程下同一 y[i] 会被多个线程写，产生数据竞争。因此实际工程中通常采用行优先，并在此基础上做行分块、按行切分线程、对 x 做分块预取。

要量化差异，可在 main 中用 clock_gettime(CLOCK_MONOTONIC, ...) 包裹两种 SpMV，各重复 100 次取平均。判读原则：若矩阵每行非零元很少且列号分布分散，行优先的 x 间接访问会成为主要开销；若矩阵带宽很窄（非零元集中在少数列），x 的工作集小、缓存命中高，两种顺序差距会明显缩小。不要预设某一种一定更快，应以实测为准。