3. 稀疏矩阵向量乘（SpMV）的访存优化

稀疏矩阵向量乘（SpMV）计算 y = A·x，其中 A 是 m×n 的稀疏矩阵，x 是长度为 n 的稠密向量，y 是长度为 m 的结果向量。它的浮点运算量只有 2·nnz（nnz 为非零元个数），但需要读取 nnz 个矩阵非零值、nnz 个列索引以及 x 中对应的元素。以 CSR 格式为例，一次 SpMV 的访存字节数约为 nnz×(sizeof(val)+sizeof(col_idx)) + (m+1)×sizeof(row_ptr) + n×sizeof(x) + m×sizeof(y)。当矩阵规模达到百万级非零元时，算术强度通常只有 0.1～0.3 FLOP/Byte，远低于现代处理器的机器平衡点，因此 SpMV 几乎总是访存带宽受限而非计算受限。优化的核心目标不是减少浮点运算，而是让每一次内存事务尽可能被充分利用，并让访存与计算更好地重叠。

CSR 格式的访存模式分析

CSR 用三个数组描述矩阵：row_ptr[m+1] 记录每行第一个非零元在 val 与 col_idx 中的偏移，col_idx[nnz] 存储列号，val[nnz] 存储非零值。计算第 i 行时，代码遍历 k = row_ptr[i] … row_ptr[i+1]-1，执行 y[i] += val[k] * x[col_idx[k]]。

这里存在两类截然不同的访存行为。流式访问：val、col_idx、row_ptr 都是按 k 递增顺序读取的连续数组，硬件预取器容易识别，缓存行利用率高。间接访问（gather）：x[col_idx[k]] 的地址由列索引决定，列号在行内通常无序且跨度大，一次缓存行加载可能只用到其中 4 字节，造成严重的缓存行浪费；同时它是指针追逐式的依赖加载，延迟难以被流水线掩盖。因此，SpMV 的瓶颈往往集中在 x 的随机收集上，而非矩阵数组本身。

另一个容易被忽视的问题是行间负载不均。若采用“一线程一行”的并行划分，非零元分布极不均匀的矩阵（如幂律图）会让少数线程承担绝大部分工作，其余线程空转，带宽利用率被拖低。这也是后续优化必须与负载划分一起考虑的原因。

循环展开：先解决“请求不够多”

基础实现中，每次迭代只发出一次 x 的收集请求，而下一次迭代的地址依赖 col_idx[k+1] 的加载结果，形成串行链。循环展开把一次迭代处理多个非零元，例如展开 4 次：先一次性读出 col_idx[k..k+3]，再发出 4 个独立的 x 收集，最后用 4 个乘加更新同一个 y[i]。由于 4 个收集请求之间没有依赖，处理器可以同时挂起多个未完成的内存事务，内存级并行度（MLP）显著提高。

实现时要注意尾部处理：当 row_ptr[i+1]-row_ptr[i] 不是展开因子的整数倍时，剩余元素用标量循环收尾。展开因子并非越大越好，4～8 通常是合理区间；过大时寄存器压力上升，编译器可能把累加器溢出到栈上，反而增加访存。判断展开是否生效，可以观察编译后的汇编中 gather 加载是否成组出现，以及性能计数器里 mem_inst_retired 与 cycles 的比值是否上升。

寄存器分块：让 x 的收集被多行复用

寄存器分块（register blocking）一次处理 R 行，把这 R 行的部分和保存在 R 个寄存器累加器中。对每一列位置，若多行共享相同的列索引，则只需一次 x 收集即可服务多个乘加。对于块状或带状结构的矩阵，这种复用非常有效；对于列索引完全随机的矩阵，复用机会少，但仍能通过同时推进多行使多行指针读取与 gather 交错，提升流水线占用。

一个实用的做法是采用固定 R×C 的小块划分（如 4×4），把行内连续的非零元按列号对齐成小块，块内用展开的乘加处理，块间用循环推进。这样既控制了代码体积，又让编译器能把块内累加器全部放在寄存器中。需要强调的是，分块改变了遍历顺序，必须保证每个 y[i] 的累加顺序不影响结果正确性——浮点加法不满足结合律，若对结果有严格可复现要求，应固定累加顺序或使用补偿求和。

软件预取：把延迟藏到计算后面

预取的关键是“提前量”。对 val 与 col_idx，由于是顺序访问，硬件预取通常已足够，手动预取收益有限；真正需要干预的是 x 的收集。可以在处理第 k 个元素时，读取 col_idx[k+D] 并对其指向的 x 地址发出预取，D 为预取距离，常取 8～32 个元素。D 太小无法覆盖内存延迟，太大则预取的数据在真正使用前已被挤出缓存，还会挤占带宽。

for (int i = 0; i < m; ++i) {
  double sum = 0.0;
  int start = row_ptr[i], end = row_ptr[i + 1];
  for (int k = start; k < end; ++k) {
    int kp = k + PREFETCH_DIST;
    if (kp < end) {
      __builtin_prefetch(&x[col_idx[kp]], 0, 1);
    }
    sum += val[k] * x[col_idx[k]];
  }
  y[i] = sum;
}

预取距离应通过实验扫描确定，而不是照搬常数。不同矩阵的平均行长度差异很大：行很短的矩阵，行内预取窗口很快越界，收益有限，此时更应考虑跨行的预取或改用分块策略。

测量：带宽利用率与加速比怎么看

评估优化效果需要两类指标。第一类是有效带宽：按实际搬运字节数除以耗时计算，BW = bytes_moved / time，再与平台可持续内存带宽（可用 STREAM 类基准测得）相比得到带宽利用率。第二类是加速比：以基础 CSR 实现为基准，speedup = T_base / T_opt。两者必须一起看——若加速比提升但有效带宽已接近上限，说明优化空间基本用尽；若带宽利用率仍低而加速比停滞，则瓶颈可能在负载不均或 gather 延迟未被掩盖。

下面是一组示例性的测量记录（数据为教学假设，非实测）：矩阵取 100 万行、约 800 万非零元，双精度，单线程运行于可持续带宽约 40 GB/s 的平台。基础 CSR 耗时 12.0 ms，有效带宽约 11.3 GB/s，利用率约 28%；展开 4 次后耗时 8.6 ms，加速比 1.40；叠加 4 行寄存器分块后耗时 7.1 ms，加速比 1.69；再加入距离 16 的 x 预取后耗时 6.4 ms，加速比 1.88，有效带宽约 21.2 GB/s，利用率约 53%。这组数字说明：展开主要改善 MLP，分块改善复用，预取改善延迟隐藏，三者作用点不同，收益会逐步递减而非线性叠加。

把上述手段落到工程上，建议按“先展开、再分块、后预取”的顺序推进，每步单独测量并记录有效带宽与加速比，避免多个变量同时变化导致无法归因。同时始终保留一个正确性校验：用随机稠密向量与参考实现比对 y 的最大绝对误差，确认优化没有改变计算结果。