科学计算中大量问题——有限元、图分析、流体、电磁——最终都归结为求解大规模稀疏线性系统 Ax = b。稀疏矩阵的绝大多数元素为零，存储与运算必须跳过这些零，否则内存和算力会被白白浪费。稀疏线性代数的核心难点不在于数学，而在于数据布局与访存效率：同一个矩阵用不同格式存储，SpMV 的性能可以差好几倍。本文先讲存储格式的内存语义，再剖析 SpMV 为何访存受限，最后落到优化手段与求解器选型。
稀疏矩阵存储格式
COO：坐标格式
坐标格式（Coordinate, COO）用三个数组描述非零元：行索引、列索引、数值。
typedef struct {
int *row; /* 非零元的行号，长度 nnz */
int *col; /* 非零元的列号，长度 nnz */
double *val; /* 非零元数值，长度 nnz */
int nnz; /* 非零元个数 */
int nrows, ncols;
} COO;
COO 结构简单、易于构造与合并，是组装阶段（从单元刚度矩阵拼装全局矩阵）的首选，但随机访存和无法直接按行遍历使其不适合作为计算格式。
CSR：压缩稀疏行
压缩稀疏行（Compressed Sparse Row, CSR）是最主流的计算格式，用行指针消除行索引数组：
typedef struct {
double *val; /* 非零元数值，按行主序，长度 nnz */
int *col; /* 每个非零元的列号，长度 nnz */
int *rowptr; /* 行起始偏移，长度 nrows+1 */
int nnz, nrows, ncols;
} CSR;
第 i 行的非零元位于 val[rowptr[i] .. rowptr[i+1]-1]。行 i 与列 j 的取值：
double csr_get(const CSR *A, int i, int j) {
for (int k = A->rowptr[i]; k < A->rowptr[i + 1]; k++)
if (A->col[k] == j) return A->val[k];
return 0.0;
}
| 格式 | 存储开销 | 适用场景 |
|---|---|---|
| COO | 3 数组 | 组装、格式转换 |
| CSR | val + col + rowptr | 通用、行操作、SpMV |
| CSC | val + row + colptr | 列操作、转置乘 |
| ELL | 等长行存储 | 行长度均匀、GPU |
| DIA | 对角带 | 带状矩阵 |
| BSR | 分块 CSR | 块状结构、多右端 |
ELL 与 DIA：规整化换取向量化
ELL 把每行补齐到最长行的长度，形成二维规则数组，天然适合向量化与 GPU 线程映射：
CSR 行长度: [3, 1, 4, 2] ELL (padded to 4):
val: [a b c | d | e f g h | i j]
col: [0 2 5 | 1 | 0 3 4 7 | 2 6]
val[4][4]: col[4][4]:
a b c * 0 2 5 -
d - - - 1 - - -
e f g h 0 3 4 7
i j - - 2 6 - -
补齐用 -（无效槽）填充，列索引可用特殊值标记。ELL 的代价是填充率：若各行长度差异大，填充会浪费大量带宽。于是有 ELL+COO 混合格式——把超长行剥离到 COO，其余进 ELL，兼顾规整与稀疏。
DIA 则把矩阵按对角线存储，对带状矩阵（如有限差分）零填充最少，但要求非零元集中在对角线附近。
BSR：分块稀疏行
当矩阵具有块结构（如每个网格节点带多个自由度），BSR 把 r×c 的小稠密块作为单位存储：
typedef struct {
double *val; /* 块值，长度 nnzb * r * c */
int *col; /* 块列号 */
int *rowptr; /* 块行指针 */
int r, c, nnzb, nrows;
} BSR;
块内是稠密小矩阵，可用 SIMD 一次处理多个元素，是结构力学、多物理场问题的标准格式。
格式转换：COO → CSR
从组装（COO）到计算（CSR）的转换是每个稀疏求解流程的第一步，标准做法是计数排序（counting sort），复杂度 O(nnz + nrows)：
void coo_to_csr(const COO *A, CSR *B) {
int n = A->nrows;
B->rowptr = calloc(n + 1, sizeof(int));
/* 第一步：统计每行非零元个数 */
for (int k = 0; k < A->nnz; k++)
B->rowptr[A->row[k] + 1]++;
/* 第二步：前缀和得到行起始偏移 */
for (int i = 0; i < n; i++)
B->rowptr[i + 1] += B->rowptr[i];
/* 第三步：按行散射非零元 */
B->val = malloc(A->nnz * sizeof(double));
B->col = malloc(A->nnz * sizeof(int));
int *next = malloc(n * sizeof(int));
memcpy(next, B->rowptr, n * sizeof(int));
for (int k = 0; k < A->nnz; k++) {
int i = A->row[k];
int pos = next[i]++;
B->val[pos] = A->val[k];
B->col[pos] = A->col[k];
}
free(next);
B->nnz = A->nnz; B->nrows = n; B->ncols = A->ncols;
}
若 COO 未排序，同一行内列号可能乱序，转换后需对每行做一次局部排序，或在 SpMV 前接受乱序（数值结果不变，但访存局部性下降）。
去重与压缩
有限元组装会产生重复的非零元（同一 (i,j) 被多个单元写入）。CSR 转换后需归并相同列号：
组装结果（含重复）: (0,1)=a (0,1)=b (0,3)=c
归并后: (0,1)=a+b (0,3)=c
MPI_Allreduce 或哈希表都能完成归并；并行组装时，通常先在本地归并再跨进程通信，减少消息量。
SpMV 的计算特征
稀疏矩阵向量乘 y = A·x 是几乎所有 Krylov 迭代法每一步的核心，其性能直接决定求解器速度。
/* CSR 上的串行 SpMV */
void spmv_csr(const CSR *A, const double *x, double *y) {
for (int i = 0; i < A->nrows; i++) {
double sum = 0.0;
for (int k = A->rowptr[i]; k < A->rowptr[i + 1]; k++)
sum += A->val[k] * x[A->col[k]];
y[i] = sum;
}
}
算术强度极低
设每行平均非零元数为 nnzr。每读一个非零元，需读取 val（8 字节）与 col（4 字节），并间接读取 x[col]（8 字节），但只做一次乘加（2 flops）。算术强度约为：
AI ≈ 2 flops / (8 + 4 + 8) bytes ≈ 0.1 flop/byte
远低于现代 CPU 的机器平衡点（通常 5~10 flop/byte）。这意味着 SpMV 是**访存受限（memory-bound）**内核：性能天花板由内存带宽决定，而非浮点算力。用 Roofline 模型 分析，SpMV 稳稳落在斜率为带宽的那条线上。
间接访存与不规则性
x[A->col[k]] 是间接访存（gather）：列索引不连续，缓存行利用率低。若矩阵带宽窄、列号聚集，x 的访问有局部性；若是随机稀疏结构（如无标度图），每次 gather 都可能触发缓存缺失。这一特性与 内存层次结构
中的 cache/TLB 行为紧密相关。
SpMV 经典优化
寄存器分块
把 CSR 按 r×c 小方块重排（register blocking），一个 x 分量被复用 r 次，减少 x 的重复加载：
/* 2x2 寄存器分块：一次处理两行两列 */
for (int i = 0; i < n; i += 2) {
double y0 = 0, y1 = 0;
for (int k = 0; k < rowblocks; k++) {
double x0 = x[col[k]], x1 = x[col[k] + 1];
y0 += val[k][0] * x0 + val[k][1] * x1;
y1 += val[k][2] * x0 + val[k][3] * x1;
}
y[i] = y0; y[i + 1] = y1;
}
分块把「间接访存比例」降低，但引入了显式零填充，需在带宽节省与填充浪费之间权衡。
SELL-C-σ：按行长度排序
SELL-C-σ 按行长度排序后切分成固定高度的切片（slice），每个切片内部用 ELL 布局，切片长度自适应：
按 nnzr 排序后:
slice 0 (短行): 宽度 2
slice 1 (中行): 宽度 4
slice 2 (长行): 宽度 7
排序重排了行顺序，需同步置换 x 与 y，但换来极高的向量化效率与低填充率，是 CPU/GPU 通用性最好的格式之一。
Merge-based（CSR-Stream）
merge-based SpMV 把「行划分」改成「非零元划分」：把所有非零元均分成 P 段，每段独立处理，避免行长不均导致的负载失衡。结合 merge-path 算法在线程间分配工作：
总非零元 nnz 均匀切成 P 段，每段 (nnz/P) 个非零元
线程 t 处理 [t·nnz/P, (t+1)·nnz/P)，跨行边界由归并定位
这对 GPU 尤其重要——避免「一行一个线程」在长行上退化为串行。
向量化与对齐
CSR 的 col 用 32 位整数、val 用 64 位双精度，混排会破坏 SIMD 对齐。可拆分数组（SoA）并对齐到 32/64 字节边界：
double *val __attribute__((aligned(64)));
int *col __attribute__((aligned(64)));
配合 #pragma omp simd 让编译器对 sum += val[k] * x[col[k]] 向量化——注意 col[k] 的 gather 在多数架构上无法直接 SIMD 化，需依赖硬件 gather 指令（AVX2/AVX-512）或改写为分块稠密形式。
多核与 GPU 实现
| 平台 | 推荐格式 | 并行策略 |
|---|---|---|
| 多核 CPU | CSR + SELL-C-σ | OpenMP 按行/按非零元划分 |
| NVIDIA GPU | ELL / hybrid ELL+COO | 一 warp 一行或 merge-based |
| AMD GPU | CSR / BSR | wavefront 映射 |
| 多右端 | BSR / dense block | 块内 GEMM |
GPU 上的 SpMV 关键是避免线程发散：ELL 让所有线程循环次数一致，但长行差异大时填充浪费；merge-based 让所有线程处理相同数量的非零元，负载完全均衡，是现代库（如 cuSPARSE 的 SpMV 新算法、Ginkgo）的主流方案。多右端场景（Y = A·X）本质是稀疏-稠密矩阵乘，可按 BSR 块用 GEMM 处理，算术强度显著提升。
稀疏矩阵的其他核心运算
SpMV 之外，稀疏线性代数还有几个高频内核，它们共享「间接访存 + 不规则」的挑战。
SpGEMM：稀疏矩阵乘稀疏矩阵
稀疏-稀疏矩阵乘 C = A·B 的结果规模事先未知，需要符号阶段预测 C 的非零结构、数值阶段填充：
符号阶段: 对每行 i, 合并 A 行 i 的非零列对应的 B 行结构 → C 行 i 的模式
数值阶段: 按模式累加 a_ik * b_kj
SpGEMM 用于 AMG 的粗化算子构造、图算法的多步传播，内存分配是主要难点——常用「两遍法」先算大小再分配。
SpTRSV：稀疏三角求解
L y = b（下三角）是 ILU 预条件的核心。第 i 个分量依赖前面分量，天然串行：
for (int i = 0; i < n; i++) {
double sum = b[i];
for (int k = rowptr[i]; k < rowptr[i + 1]; k++)
if (col[k] < i) sum -= val[k] * y[col[k]];
y[i] = sum / diag[i];
}
并行化靠层次划分（level scheduling）：按依赖关系把行分层，同层内无依赖可并行。多色排序（multicoloring）是生成层次的常用方法。
转置与稀疏结构变换
CSC 与 CSR 互为转置，转置等价于一次 COO→CSR 重排。转置乘 Aᵀ·x 在 BiCGSTAB 等算法中需要，可在 CSR 上直接用「散射」形式实现而无需显式转置：
/* y = Aᵀ x，CSR 上直接散射 */
memset(y, 0, ncols * sizeof(double));
for (int i = 0; i < nrows; i++)
for (int k = rowptr[i]; k < rowptr[i + 1]; k++)
y[col[k]] += val[k] * x[i];
散射写 y 存在写冲突，多线程下需原子操作或按列划分；这也是转置乘通常比正向 SpMV 更难并行化的原因。
Krylov 求解器与预条件
SpMV 本身不求解方程，它嵌在 Krylov 迭代法中。理解求解器才能理解 SpMV 优化的收益边界。
| 方法 | 适用矩阵 | 每步 SpMV 数 | 内存 |
|---|---|---|---|
| CG | 对称正定 | 1 | 少量向量 |
| BiCGSTAB | 非对称 | 2 | 少量向量 |
| GMRES(m) | 非对称 | 1 + 正交化 | 重启向量组 |
| MINRES | 对称不定 | 1 | 少量向量 |
预条件子（Preconditioner）把病态系统转化为良态系统，减少迭代次数，但每步可能引入额外的稀疏三角求解（前向/后向替换），这部分是串行瓶颈：
Jacobi (对角): 完全并行, 效果弱
ILU(0): 串行/层次并行, 效果中
AMG (代数多重网格): 复杂, 效果强, 近似 O(n) 收敛
ILU 的前向替换 L y = r 依赖前面行的结果，天然串行；多色排序（graph coloring）可解除部分依赖，实现并行三角求解。对于大规模问题，代数多重网格（AMG） 常作为最优预条件，收敛步数与问题规模近似无关。完整的求解器库选型与配置见 PETSc 求解器
一文。
性能评估与选型
评估 SpMV 性能应看有效带宽而非 flops：
# 用 STREAM 测机器带宽上限
stream
# 计算 SpMV 有效带宽 = (nnz*(4+8) + nrows*8 + ncols*8) / time
| 优化手段 | 典型收益 | 代价 |
|---|---|---|
| 寄存器分块 | 1.2~2× | 零填充、格式转换 |
| SELL-C-σ | 1.5~2.5× | 行重排、置换开销 |
| merge-based | 1.3~2× (GPU) | 实现复杂 |
| 混合精度（val 存 FP32） | 1.5× | 精度损失 |
| 多右端批处理 | 2~5× | 内存占用 |
经验法则：先测量矩阵的行长度分布与带宽分布——行长度方差大选 SELL-C-σ，带状选 DIA，块状选 BSR，通用选 CSR。格式与算法的联合设计（如为 AMG 层级选择不同格式）比单点优化收益更大，这与 并行算法设计 中「负载均衡优先于局部优化」的原则一致。
小结
稀疏线性代数的性能之争，本质是内存布局之争。CSR 提供通用性，ELL/SELL-C-σ 用规整化换取向量化，BSR 用块结构换取 SIMD，merge-based 用非零元划分换取负载均衡。没有万能格式，只有与矩阵结构和硬件匹配的格式。掌握了存储格式的内存语义与 SpMV 的访存受限本质，你就能在求解器选型与内核优化之间做出有依据的权衡，而不是盲目套用某个库的默认配置。
继续阅读
探索更多技术文章
浏览归档，发现更多关于系统设计、工具链和工程实践的内容。