一个支持分布式(按连续行块划分)与对称存储扩展的 C++17 CSR 稀疏矩阵库,集成 MPI(可选)、OpenMP(可选)、GoogleTest 单元测试与 MUMPS 导出适配。
核心设计:模板化多标量类型支持 —— 库已完全模板化,四种标量类型(float、double、std::complex<float>、std::complex<double>)可在同一次编译中同时使用,无需切换编译开关重新构建。提供远程条目通信装配与对称矩阵的上下三角展开。
当前构建为静态库
csr4mpi+ 测试可执行csr4mpi_tests。
项目完全使用 Github Copilot 开发,用时 3 小时左右。
- 模板化多标量类型:所有核心类均为 C++ 模板,支持
float、double、std::complex<float>、std::complex<double>四种标量类型,无需重新编译即可同时使用。 - 便利类型别名:提供
cCSRMatrixF、cCSRMatrixD、cCSRMatrixCF、cCSRMatrixCD等预定义别名,简化常用类型的声明。 - 类型特征工具:
is_complex_v<T>、is_supported_scalar_v<T>、real_type_t<T>等编译期类型检测和转换工具。 - 连续行块分布:
cRowDistribution描述全局行到进程的映射。 - 本地 CSR 结构:
cCSRMatrix<Scalar>保存所属行区间的行指针、列索引与数值。 - 通信模式:
cCommPattern<Scalar>聚合跨进程目标条目 (globalRow, globalCol, value)。 - 汇总装配:
cCSRComm<Scalar>::Assemble使用 MPI Alltoallv 交换并累加远程重复条目。 - 对称存储展开:支持下三角存储在 SpMV/SpMM 中自动补全镜像贡献(本地 + 远程)。
- MUMPS 导出:
cMumpsAdapter<Scalar>输出 1-based COO 三元组。 - PETSc 原生预条件子(可选,
CSR4MPI_ENABLE_PETSC=ON):cPetscSolver<Scalar>支持 PETSc 全部预条件子类型(字符串透传,如jacobi/bjacobi/sor/ilu/lu/icc/cholesky/asm/gasm/gamg/mg/redundant等),含完整 KSP 求解与 PC-only 应用(ApplyPreconditionerOnly)。 - Benchmark:
csr4mpi_bench_spmv输出平均耗时与 GFLOPs。
推荐在独立目录:
mkdir build
cd build
cmake -G "Ninja" ..
cmake --build . --config Release参数说明:
-DCSR4MPI_ENABLE_OPENMP=ON|OFF:控制并行内核的 OpenMP 支持。- 无需指定标量类型参数 —— 四种标量类型在同一次编译中均可使用。
提供脚本 scripts/download_matrices.sh(依赖 bash + curl + tar):
bash scripts/download_matrices.sh将尝试下载:
bcsstk13.mtx(Boeing group) – 中等规模结构矩阵add32.mtx(HB group) – 较小测试矩阵
文件放置到 tests/data/。若访问失败请手动前往 https://sparse.tamu.edu/ 对应页面下载 tar.gz 后解压并将 .mtx 文件置于该目录。
也支持通过 scripts/matrices.txt 指定自定义下载清单(每行一个 .tar.gz 完整 URL;忽略以 # 开头的行):
echo "https://your-storage/MM/Boeing/bcsstk13.tar.gz" >> scripts/matrices.txt
echo "https://your-storage/MM/HB/add32.tar.gz" >> scripts/matrices.txt
bash scripts/download_matrices.sh构建后在 build:
ctest -C Release --output-on-failure多进程(示例 4 进程)针对单独的分布式测试:
mpiexec -n 4 .\csr4mpi_tests.exe --gtest_filter=MpiMatrixMarketLargeTest.DistributedLargeSpMVAccuracy
mpiexec -n 4 .\csr4mpi_tests.exe --gtest_filter=MpiMatrixMarketLargeTest.DistributedDuplicateAssembly
mpiexec -n 4 .\csr4mpi_tests.exe --gtest_filter=MpiSymmetricSpMMTest.DistributedLowerExpansionGatherMM4Proc大型矩阵测试回退逻辑:
- 首选外部下载的
bcsstk13.mtx。 - 若缺失自动尝试合成矩阵
synthetic_12k12_general.mtx(需先使用生成器创建)。 - 两者都不可用时测试标记为 skip 并输出生成提示。
生成兜底矩阵示例:
.\csr4mpi_gen_mm.exe 12000 12000 12 0 .\tests\matrices\synthetic_12k12_general.mtx构建时开启:-DCSR4MPI_ENABLE_PETSC=ON(自动经 pkg-config 或 PETSC_DIR 查找 PETSc)。
#include "PetscSolver.h"
csr4mpi::cCSRMatrix<double> A = /* 分布式 CSR */;
csr4mpi::cPetscSolver<double> solver;
csr4mpi::cPetscSolver<double>::Options opt;
opt.sKspType = "cg"; // cg / gmres / preonly / ...
opt.sPcType = "gamg"; // 任意 PETSc PC 类型字符串(透传)
opt.dRelTol = 1e-10;
opt.sPrefixOptions = "-pc_gamg_type agg"; // 附加 PETSc 选项
// Hypre(LLNL)预条件子:需 PETSc 以 --download-hypre 或系统 Hypre 构建。
// BoomerAMG 是大规模 SPD 系统的首选 AMG 之一,用法:
csr4mpi::cPetscSolver<double> hsolver;
csr4mpi::cPetscSolver<double>::Options hopt;
hopt.sKspType = "cg";
hopt.sPcType = "hypre";
hopt.sPrefixOptions = "-pc_hypre_type boomeramg"; // parasails / pilut / euclid / ams 亦可
// if (hsolver.Create(A, MPI_COMM_WORLD, hopt) && hsolver.Setup()) hsolver.Solve(x, b);
if (solver.Create(A, MPI_COMM_WORLD, opt) && solver.Setup()) {
std::vector<double> x, b /* 本地段 */;
solver.Solve(x, b); // 完整 KSP 求解
// solver.ApplyPreconditionerOnly(x, b); // 仅应用预条件子(PCApply)
}说明:
-
预条件子类型为字符串透传,任意 PETSc 注册的 PC(含外部包 Hypre/ML/MUMPS)均可直接使用;非法类型
Create干净返回false(内部在MPI_COMM_SELF上探测,多进程下不会触发 PETSc abort)。 -
多进程下
lu/cholesky等分解型 PC 需外部分解器(如-pc_factor_mat_solver_type mumps);PETSc 不提供 mpiaij 的原生并行 ILU(可改用bjacobi)。 -
矩阵约束:
MatCreateMPIAIJWithArrays要求每行的列索引按升序排列(CreateMat沿用此约束);乱序会导致装配缺对角元等隐蔽错误。 -
Hypre 预条件子:
sPcType="hypre"+sPrefixOptions="-pc_hypre_type boomeramg|parasails|pilut|euclid|ams"。需要 PETSc 编译时包含 Hypre(--download-hypre);缺少后端时Setup()返回false,测试自动跳过。BoomerAMG 与gamg同为大规模并行 AMG 首选,建议基准对比后择优。 -
Hypre 后端安全探测:
cPetscSolver::bIsHypreBackendAvailable("boomeramg")在MPI_COMM_SELF上经带前缀的临时求解器探测(不依赖 Hypre 编译期符号,MSYS2 等未编 Hypre 的构建也能链接),多进程下调用同样安全。推荐写法:auto opt = /* ... */; if (csr4mpi::cPetscSolver<double>::bIsHypreBackendAvailable("boomeramg")) { opt.sPcType = "hypre"; opt.sPrefixOptions = "-pc_hypre_type boomeramg"; } else { opt.sPcType = "gamg"; // 内置 AMG 回退 }
-
PC-only 预条件子可用性:
cPetscSolver::bIsPreconditionerTypeAvailable("hypre")可在运行时安全探测类型。 -
并行性能参考(N=262144 五点 Laplacian,4 进程,CG):
hypre-boomeramg0.22s/6 次迭代(最优);gamg0.24s/10 次迭代(近线性扩展);bjacobi/gasm3.7-4.0s;lu(MUMPS) 0.07s 直接解;jacobi/none8.8s。基准工具:csr4mpi_bench_petsc_pc(CSR4MPI_BUILD_BENCHMARKS=ON时构建)。 -
PETSc scalar 类型需与模板实参宽度一致(如 real PETSc +
double),否则Create返回false。 -
测试:
mpiexec -n 4 ./csr4mpi_tests --gtest_filter='PetscPrecond.*'。
可执行:csr4mpi_bench_spmv
mpiexec -n 4 .\csr4mpi_bench_spmv.exe bcsstk13.mtx 20输出示例字段:
avg / max / min:多次迭代(默认 10)后进程归约的平均/最大/最小单次耗时(秒)。nnz:全局非零元数。GFLOPs:采用公式2 * nnz / time / 1e9(一次乘加视为 2 浮点运算)。
非 MPI 下直接运行:
./csr4mpi_bench_spmv.exe bcsstk13.mtx 10也可使用本地合成矩阵(不依赖外部下载),通过生成器 csr4mpi_gen_mm:
.\csr4mpi_gen_mm.exe 8000 8000 12 0 .\tests\matrices\synthetic_8k12_general.mtx
.\csr4mpi_gen_mm.exe 12000 12000 10 1 .\tests\matrices\synthetic_12k10_sym_lower.mtx
mpiexec -n 4 .\csr4mpi_bench_spmv.exe synthetic_8k12_general.mtx 20
mpiexec -n 4 .\csr4mpi_bench_spmv.exe synthetic_12k10_sym_lower.mtx 20库采用 C++ 模板实现多标量类型支持,四种类型在同一次编译中均可使用,无需切换编译开关:
| 别名 | 完整类型 |
|---|---|
cCSRMatrixF |
cCSRMatrix<float> |
cCSRMatrixD |
cCSRMatrix<double> |
cCSRMatrixCF |
cCSRMatrix<std::complex<float>> |
cCSRMatrixCD |
cCSRMatrix<std::complex<double>> |
#include "CSRMatrix.h"
#include "Operations.h"
using namespace csr4mpi;
// 所有类型在同一程序中同时可用
cCSRMatrixD matDouble(0, n, n, rowPtr, colInd, valuesDouble);
cCSRMatrixCF matComplexFloat(0, n, n, rowPtr, colInd, valuesComplexFloat);
// SpMV 操作
std::vector<double> xD(n, 1.0), yD;
SpMV(matDouble, xD, yD);
std::vector<std::complex<float>> xCF(n, {1.0f, 0.0f}), yCF;
SpMV(matComplexFloat, xCF, yCF);
// 也可以直接使用模板形式
cCSRMatrix<std::complex<double>> matCD(...);Global.h 提供以下编译期类型检测与转换工具:
is_complex_v<T>:判断T是否为std::complex<>类型。is_supported_scalar_v<T>:判断T是否为支持的四种标量之一。real_type_t<T>:获取标量的实部类型(实数类型返回自身,复数返回其分量类型)。
static_assert(is_complex_v<std::complex<double>>); // true
static_assert(!is_complex_v<double>); // true (double is not complex)
static_assert(is_supported_scalar_v<float>); // true
static_assert(std::is_same_v<real_type_t<std::complex<float>>, float>); // true若矩阵以下三角形式存储(Matrix Market symmetric lower),在分布式 SpMV/SpMM 中:
- 本地展开:补全镜像列贡献。
- 远程展开:对需要的镜像条目进行分桶发送与累加,避免重复计算。
频域含损耗/复杂材料的电磁 FEM 系统矩阵通常为复对称(complex symmetric,非 Hermitian)。将本库导出的矩阵交给下游求解器时请注意:
- CG 仅适用于 Hermitian 正定矩阵;对复对称系统直接使用 CG 收敛慢且病态加剧。
- 复对称系统推荐 BiCGSTAB / GMRES 等迭代法,或 MUMPS/PARDISO 直接法(可利用对称性加速,对称分解较通用 LU 收益约 2×)。
- 与 HFSS 对标时注意:HFSS 的 symmetry 边界是几何级减半建模(非矩阵存储级对称),其官方接口(ExportNetworkData/Touchstone 等)只导出端口 S/Y/Z 网络矩阵,不提供 FEM 稀疏系统矩阵导出——对标只能做"解级"对比。
- 本库的按行块 MPI 分布 + halo 交换 SpMV 与 HFSS DDM 的分布式内存思路一致,但通信模式不同(规则 halo 交换 ≠ 子域界面耦合),不应宣称"等同 DDM"。
include/Global.h:类型定义、类型特征(is_complex_v、is_supported_scalar_v、real_type_t)与 MPI 辅助函数。include/CSRMatrix.h:模板化 CSR 本地块cCSRMatrix<Scalar>与类型别名。include/Distribution.h/src/Distribution.cpp:行分布。include/CommPattern.h:模板化远程条目通信模式cCommPattern<Scalar>。include/CSRComm.h:模板化汇总装配逻辑cCSRComm<Scalar>。include/MumpsAdapter.h:模板化 MUMPS 导出cMumpsAdapter<Scalar>。include/PardisoAdapter.h/include/PetscAdapter.h/include/TrilinosAdapter.h:PARDISO / PETSc / Trilinos 导出适配器。include/Operations.h:模板化 SpMV/SpMM 操作。include/MatrixMarketLoader.h:模板化 Matrix Market 文件加载器。include/DistributedOps.h:分布式 SpMV 操作。bench/bench_large_spmv.cpp:SpMV 性能基准。tests/:GoogleTest 单元与 MPI 分布式测试。tests/test_all_scalar_types.cpp:针对四种标量类型的全面测试(52 个测试用例)。
- 更高效装配(行内列索引二分/哈希)
- SpMM 进一步通信压缩与重用
- 更多矩阵与可复现基准统计 (CSV)
本项目基于 MIT License 开源。