CSR格式SpMV性能优于MKL但SpMM落后，求优化方案

一、自研实现与性能数据 #

1. 稀疏矩阵-向量乘积（SpMV）实现 #

void dcsrmv(SparseMatrixCSR *A, double *x, double *y) {
  for (int i=0; i<A->m; i++) {
    double sum = 0.;
    for (int jj=A->row_start[i]; jj<A->row_start[i+1]; jj++) {
      sum += A->val[jj] * x[A->col_idx[jj]];
    }
    y[i] = sum;
  }
}

测试结果：在SuiteSparse的ecology1矩阵上运行100次，平均耗时5.9ms，优于MKL的mkl_sparse_d_mv（8.4ms）。

2. 三种稀疏矩阵-多向量乘积（SpMM）实现 #

X、Y按列存储，当nvec=10时的性能数据：

dcsrmm0（向量循环外层）

void dcsrmm0(SparseMatrixCSR *A, double *X, double *Y, int nvec) {
  for (int vec=0; vec<nvec; vec++) {
    double *x = X + vec * A->n;
    double *y = Y + vec * A->m;
    for (int i=0; i<A->m; i++) {
      double sum = 0.;
      for (int jj=A->row_start[i]; jj<A->row_start[i+1]; jj++) {
        sum += A->val[jj] * x[A->col_idx[jj]];
      }
      y[i] = sum;
    }
  }
}

耗时：49.8ms

dcsrmm1（行循环外层，向量循环中层）

void dcsrmm1(SparseMatrixCSR *A, double *X, double *Y, int nvec) {
  for (int i=0; i<A->m; i++) {
    for (int vec=0; vec<nvec; vec++) {
      double sum = 0.;
      double *x = X + vec * A->n;
      double *y = Y + vec * A->m;
      for (int jj=A->row_start[i]; jj<A->row_start[i+1]; jj++) {
        sum += A->val[jj] * x[A->col_idx[jj]];
      }
      y[i] = sum;
    }
  }
}

耗时：38.3ms

dcsrmm2（行循环外层，非零元循环中层，向量循环内层）

void dcsrmm2(SparseMatrixCSR *A, double *X, double *Y, int nvec) {
  double *x[nvec], *y[nvec];
  for (int vec=0; vec<nvec; vec++) x[vec] = X + vec * A->n;
  for (int vec=0; vec<nvec; vec++) y[vec] = Y + vec * A->m;
  double sum[nvec];
  for (int i=0; i<A->m; i++) {
    for (int vec=0; vec<nvec; vec++) sum[vec] = 0.;
    for (int jj=A->row_start[i]; jj<A->row_start[i+1]; jj++) {
      for (int vec=0; vec<nvec; vec++) {
        sum[vec] += A->val[jj] * x[vec][A->col_idx[jj]];
      }
    }
    for (int vec=0; vec<nvec; vec++) y[vec][i] = sum[vec];
  }
}

耗时：34.8ms

MKL对比：mkl_sparse_d_mm耗时26.7ms，自研版本均落后。

二、性能差异原因分析 #

1. SpMV优于MKL的核心原因 #

• 矩阵特性适配：ecology1的稀疏结构刚好适配自研代码的行遍历逻辑，行内非零元分布让缓存命中率极高；而MKL是通用实现，要兼容所有稀疏矩阵类型，额外的分支、格式校验逻辑带来了 overhead。
• 无冗余抽象：自研代码直接操作CSR原始结构，没有MKL内部的格式转换、多层封装等额外开销。
• 编译器优化适配：极简的代码结构让编译器（如GCC/Clang）能轻松做循环展开、寄存器分配等优化，而MKL的复杂代码在特定场景下优化效果反而受限。

2. SpMM落后MKL的关键缺陷 #

MKL的mkl_sparse_d_mm针对多向量场景做了深度定制优化，自研实现存在以下致命问题：

• 内存访问模式低效：X按列存储时，x[vec][col_idx[jj]]是跨向量的随机访问，不同向量的同列元素间隔为矩阵列数A->n，远超缓存行大小，导致大量缓存 miss；dcsrmm0/1本质是重复调用SpMV，完全没利用多向量的复用性，每次都重复加载CSR的val和col_idx，浪费内存带宽。
• SIMD利用不足：dcsrmm2的内层向量循环无法被编译器有效向量化，sum[nvec]的累加没有对齐SIMD寄存器，无法发挥CPU的向量计算能力。
• 缺乏分块优化：MKL会对CSR矩阵做行分块，把非零元分布相似的行打包处理，提升缓存复用；还会对多向量做块对齐，让内存访问更连续。
• 寄存器利用率低：自研代码用栈上的sum[nvec]存储累加值，频繁读写栈内存，而MKL会尽可能把多向量累加值放在寄存器中，减少内存交互。

三、SpMM优化方案 #

1. 调整数据存储格式（优先推荐） #

如果允许修改输入格式，将X改为按行存储（同一列的nvec个向量元素连续存放），此时x[col_idx[jj]][vec]的访问是连续内存操作，缓存命中率会大幅提升。结合dcsrmm2的逻辑，能直接利用连续访问的优势。

2. 手动SIMD向量化 #

针对现有列存储X，利用编译器内置函数（如AVX/AVX512指令）手动实现向量累加，示例片段（以nvec=4、AVX2为例）：

// 假设X格式调整为：每列的nvec个元素连续存储
__m256d sum_vec = _mm256_setzero_pd();
for (int jj=A->row_start[i]; jj<A->row_start[i+1]; jj++) {
    double val = A->val[jj];
    int col = A->col_idx[jj];
    // 加载该列的4个向量元素
    __m256d x_vec = _mm256_loadu_pd(X + col * nvec);
    // val * x_vec 并累加
    sum_vec = _mm256_add_pd(sum_vec, _mm256_mul_pd(_mm256_set1_pd(val), x_vec));
}
// 将结果写入Y的第i行
_mm256_storeu_pd(Y + i * nvec, sum_vec);

3. CSR分块优化（Block CSR） #

将原CSR矩阵划分为若干行块，每个块内的非零元连续存储，对每个块的多向量做批量累加处理，提升val和col_idx的缓存复用率，减少缓存 miss。

4. 消除冗余内存操作 #

• 提前缓存x[vec]和y[vec]的指针，避免循环内重复计算；
• 用寄存器数组替代栈上的sum[nvec]，只有在每行计算完成后一次性写入Y，减少内存读写次数。

5. 编译器优化选项加持 #

编译时开启最高级优化：

• GCC/Clang：-O3 -march=native -ffast-math
• MSVC：/O2 /arch:AVX2
  这些选项会让编译器自动做循环展开、寄存器分配、SIMD向量化等优化，大幅提升性能。

内容的提问来源于stack exchange，提问作者Nicolas Venkovic