cuda-samples 深入解析：基于 cuSolverSp 低级 API 的稀疏 QR 分解与线性方程组求解（cuSolverSp_LowlevelQR）

cuda-samples 深入解析：基于 cuSolverSp 低级 API 的稀疏 QR 分解与线性方程组求解（cuSolverSp_LowlevelQR）

【免费下载链接】cuda-samplesSamples for CUDA Developers which demonstrates features in CUDA Toolkit项目地址: https://gitcode.com/GitHub_Trending/cu/cuda-samples

本文以 NVIDIA cuda-samples 仓库中的 cuSolverSp_LowlevelQR 示例 为对象，逐层剖析它如何借助 cuSolverSP 的低级（Low-level）API 对稀疏矩阵执行 QR 分解、完成A*x = b求解，并用 cuSPARSE 的 SpMV 校验残差。读完本文，你将掌握cusolverSpCreateCsrqrInfo(Host)系列 API 的完整调用链、CPU/GPU 两条求解路径的异同，以及如何用 Matrix Market（MM）格式的稀疏矩阵驱动这套流程。

示例概述：为什么需要"低级" QR 求解器

在 cuSolver 库中，求解稀疏线性系统通常有两条路线：

• 高级 API（如cusolverSpDcsrlsvqr）：一条函数调用完成"分析 → 分解 → 求解 → 残差检查"的全部工作，接口简单但对过程不可控；
• 低级 API（Low-level API）：将求解过程拆解为多个可独立调用的阶段（创建 info 结构 → 符号分析 → 工作空间查询 → setup → 数值分解 → 奇异检测 → 求解），每个阶段可复用、可调试，适合需要精细控制或反复求解同一结构矩阵的场景。

本示例正是后者，其官方定位为"A CUDA Sample that demonstrates QR factorization using cuSolverSP's low level APIs"，核心概念属于Linear Algebra（线性代数）与 CUSOLVER Library。完整实现位于 cuSolverSp_LowlevelQR.cpp，从源码结构看，它把求解过程组织成了 14 个带编号的 printf 步骤，清晰展示了低级 API 的每个阶段。

环境、架构与依赖

支持的软硬件平台

根据 README：

• 支持的 SM 架构：SM 5.0、5.2、5.3、6.0、6.1、7.0、7.2、7.5、8.0、8.6、8.7、8.9、9.0；
• 支持的操作系统：Linux、Windows；
• 支持的 CPU 架构：x86_64、armv7l；
• 依赖库：CUSOLVER（QR 分解与求解）与 CUSPARSE（CSR 矩阵描述与 SpMV 残差计算），对应仓库根目录 README.md#cusolver 与 README.md#cusparse 中的依赖说明；
• 前置条件：安装对应平台的 CUDA Toolkit。

构建配置

CMakeLists.txt 给出了可复现的构建方式：

cmake_minimum_required(VERSION 3.20) project(cuSolverSp_LowlevelQR LANGUAGES C CXX) find_package(CUDAToolkit REQUIRED) set(CMAKE_CUDA_ARCHITECTURES 75 80 86 87 89 90 100 110 120) add_executable(cuSolverSp_LowlevelQR cuSolverSp_LowlevelQR.cpp mmio.c mmio_wrapper.cpp) target_link_libraries(cuSolverSp_LowlevelQR PRIVATE CUDA::cudart CUDA::cublas CUDA::cusolver )

要点说明：

• 目标由cuSolverSp_LowlevelQR.cpp、mmio.c与mmio_wrapper.cpp三个源文件组成，其中后两者负责 Matrix Market 格式的解析与 CSR 转换；
• 默认CMAKE_CUDA_ARCHITECTURES覆盖 75（Turing）到 120（Blackwell）等主流架构；开启ENABLE_CUDA_DEBUG时追加-G以支持 cuda-gdb，否则使用-lineinfo保留行号信息；
• 通过include_directories(../../../Common)引入 helper_cuda.h 与 helper_cusolver.h 等公共头文件；
• POST_BUILD 阶段会将lap2D_5pt_n32.mtx、lap2D_5pt_n100.mtx、lap3D_7pt_n20.mtx三个数据文件复制到构建输出目录，保证运行时能定位默认输入。

输入数据：Matrix Market 格式与 CSR 转换

三个测试矩阵

仓库为示例提供了三个稀疏测试矩阵（数据文件目录）：

以 lap2D_5pt_n32.mtx 为例，其头部为：

%%MatrixMarket matrix coordinate real general % standard 5-point laplace 2D oprator with Dirichlet boundary condition 1024 1024 3008 1 1 4 2 1 -1 33 1 -1 ...

第一行 Banner 声明对象为matrix、存储为coordinate（稀疏坐标格式）、数据类型为real、结构为general；第二行1024 1024 3008分别是行数、列数、非零元个数；随后每行i j value记录一个三元组。注意这些矩阵是对称的（数值上对角占优、主对角为 4、邻接为 -1），但声明为 general 后仍需完整展开。

mmio_wrapper：从 MM 文件到 CSR

核心转换逻辑在 mmio_wrapper.cpp 的模板函数loadMMSparseMatrix<T_ELEM>中，它接受元素类型字符（'d'对应 double）、是否输出 CSR 格式、是否扩展对称矩阵等参数：

1. 调用mm_read_mtx_crd（实现于 mmio.c，声明于 mmio.h）读取 COO 三元组；
2. 若矩阵声明为symmetric/hermitian/skew且开启extendSymMatrix，则把三角部分补全为完整矩阵（对角元素不复制，非对角元素按对称/反对称/共轭规则处理）——这是本示例在main中传true的原因；
3. 用qsort按行主序（CSR）或列主序（CSC）对三元组排序；
4. 自动探测基址：若存在下标 0 则为 base-0，若出现等于矩阵维度的下标则为 base-1；
5. 通过compress_index压缩出行指针csrRowPtr，生成csrColInd与csrVal；
6. 调用verify_pattern校验 nnz 一致性、基址合法性、行内列下标单调性，防止后续进入 cuSolver 时报错。

该模板针对float、double、cuComplex、cuDoubleComplex四种元素类型做了显式实例化（文件末尾），配合cuGet特化模板完成从 double 到各类型的数值转换——这也是 README 中 Driver API 部分列出cuGet、cuComplex、cuDoubleComplex的原因。

命令行参数

程序支持三个命令行选项（定义于UsageSP/parseCommandLineArguments）：

若未提供-file，程序会通过sdkFindFilePath("lap2D_5pt_n32.mtx", argv[0])（见 helper_cuda.h）在可执行文件所在目录查找默认输入文件，找不到则报错退出。运行示例：

./cuSolverSp_LowlevelQR # 使用默认 lap2D_5pt_n32.mtx ./cuSolverSp_LowlevelQR -file=lap2D_5pt_n100.mtx # 指定输入矩阵 ./cuSolverSp_LowlevelQR -file=lap3D_7pt_n20.mtx -device=0

求解主流程：14 步低级 API 调用链

主函数（cuSolverSp_LowlevelQR.cpp）先做统一的句柄与流初始化：cusolverSpCreate、cusparseCreate、cudaStreamCreate，并通过cusolverSpSetStream/cusparseSetStream绑定到同一流；再用cusparseCreateMatDescr创建 CSR 矩阵描述符descrA，根据解析出的baseA设置CUSPARSE_INDEX_BASE_ONE或CUSPARSE_INDEX_BASE_ZERO。随后分别以CPU（Host）路径（步骤 1–8）与GPU 路径（步骤 9–14）各执行一遍完整流程。

CPU（Host）路径：步骤 1–8

GPU 路径：步骤 9–14

GPU 路径使用设备端 info 结构csrqrInfo_t与设备端 CSR 数据（d_csrValA/d_csrRowPtrA/d_csrColIndA），调用序列与 CPU 路径一一对应：

cusolverSpCreateCsrqrInfo(&d_info); // step 9 cusolverSpXcsrqrAnalysis(cusolverSpH, rowsA, colsA, nnzA, descrA, d_csrRowPtrA, d_csrColIndA, d_info); // step 10 符号分析 cusolverSpDcsrqrBufferInfo(cusolverSpH, rowsA, colsA, nnzA, descrA, d_csrValA, d_csrRowPtrA, d_csrColIndA, d_info, &size_internal, &size_chol); // step 11 工作空间查询 cusolverSpDcsrqrSetup(cusolverSpH, rowsA, colsA, nnzA, descrA, d_csrValA, d_csrRowPtrA, d_csrColIndA, zero, d_info); // step 12 cusolverSpDcsrqrFactor(cusolverSpH, rowsA, colsA, nnzA, NULL, NULL, d_info, buffer_gpu); // step 12 数值分解 cusolverSpDcsrqrZeroPivot(cusolverSpH, d_info, tol, &singularity); // step 13 奇异检测 cusolverSpDcsrqrSolve(cusolverSpH, rowsA, colsA, d_b, d_x, d_info, buffer_gpu); // step 14 求解

GPU 路径的数值分解阶段实际执行在设备上，工作空间buffer_gpu由cudaMalloc分配，大小来自 step 11 查询到的size_chol（程序会打印GPU buffer size = %lld bytes）。

关键设计：分析一次、分解多次

符号分析（XcsrqrAnalysis）与数值分解（Factor）分离是低级 API 的核心价值：对结构相同、数值不同的多个矩阵，只需做一次符号分析与工作空间查询，即可反复调用Setup + Factor + Solve，这是高级 API 难以直接提供的复用能力。

奇异判定与容差

程序在分解后调用cusolverSpDcsrqrZeroPivot(Host)检测零主元：

const double tol = 1.e-14; int singularity = 0; checkCudaErrors(cusolverSpDcsrqrZeroPivotHost(cusolverSpH, h_info, tol, &singularity)); if (0 <= singularity) { fprintf(stderr, "Error: A is not invertible, singularity=%d\n", singularity); return 1; }

源码注释明确了两点语义：singularity == -1表示在容差tol下 A 可逆；tol决定了奇异的判定条件。任何singularity >= 0的值都指向首个零主元所在位置，程序立即报错退出。三个拉普拉斯测试矩阵都是严格对角占优的 M 矩阵，因而都能通过该检查。

残差验证：cuSPARSE SpMV 通用 API

步骤 8 与 GPU 路径的收尾阶段都通过 cuSPARSE 计算残差r = b - A*x，这是对求解正确性的量化验证。关键点在于它使用了 cuSPARSE 的**通用 API（generic API）**对象而非旧式描述符：

cusparseCreateCsr(&matA, rowsA, colsA, nnzA, d_csrRowPtrA, d_csrColIndA, d_csrValA, CUSPARSE_INDEX_32I, CUSPARSE_INDEX_32I, baseA ? CUSPARSE_INDEX_BASE_ONE : CUSPARSE_INDEX_BASE_ZERO, CUDA_R_64F); cusparseCreateDnVec(&vecx, colsA, d_x, CUDA_R_64F); cusparseCreateDnVec(&vecAx, rowsA, d_r, CUDA_R_64F); cusparseSpMV_bufferSize(cusparseH, CUSPARSE_OPERATION_NON_TRANSPOSE, &minus_one, matA, vecx, &one, vecAx, CUDA_R_64F, CUSPARSE_SPMV_ALG_DEFAULT, &bufferSize); cusparseSpMV(cusparseH, CUSPARSE_OPERATION_NON_TRANSPOSE, &minus_one, matA, vecx, &one, vecAx, CUDA_R_64F, CUSPARSE_SPMV_ALG_DEFAULT, buffer);

这里r = (-1)*A*x + 1*b，先查询工作空间再执行。求得的h_r回拷 CPU 后，借助 helper_cusolver.h 中的工具函数计算指标：

• vec_norminf（L49–56）：向量无穷范数|r|；
• csr_mat_norminf（L77–96）：按行累加 CSR 矩阵各元素绝对值后取最大值，即|A|；
• 最终输出相对残差|b - A*x| / (|A|*|x|)，CPU 与 GPU 路径各打印一组：

(CPU) |b - A*x| = ...E-14 (CPU) |A| = ...E+00 (CPU) |x| = ...E+00 (CPU) |b - A*x|/(|A|*|x|) = ...E-14 (GPU) |b - A*x| = ...E-14 (GPU) |b - A*x|/(|A|*|x|) = ...E-14

残差量级达到1e-14左右，与双精度求解和tol = 1e-14的奇异容差相匹配，说明求解结果在浮点精度内满足原方程。

资源清理与工程实践

主函数结尾对全部句柄与内存做了成对释放，可作为编写生产代码的清单参考：cusolverSpDestroy、cusparseDestroy、cudaStreamDestroy、cusparseDestroyMatDescr、cusolverSpDestroyCsrqrInfoHost/cusolverSpDestroyCsrqrInfo、cusparseDestroySpMat/cusparseDestroyDnVec，以及所有cudaFree/free。这些释放操作与 README 中列出的 CUDA Runtime API（cudaMemcpy、cudaStreamDestroy、cudaFree、cudaMalloc、cudaStreamCreate）一一对应。

小结

cuSolverSp_LowlevelQR 是一个理解 cuSolverSP 低级 API 的理想样本：它用同一份 CSR 数据分别跑通 Host 与 Device 两条求解链，完整覆盖"符号分析 → 工作空间查询 → setup → 数值分解 → 奇异检测 → 求解 → SpMV 残差验证"的每个环节，并配套了 MM 格式解析、基址探测与结果范数校验等工程细节。开发者可以以此为模板，将其中"分析一次、多次分解"的模式直接复用到 CFD、结构力学、电路仿真等需要反复求解同结构稀疏线性系统的场景中。

【免费下载链接】cuda-samplesSamples for CUDA Developers which demonstrates features in CUDA Toolkit项目地址: https://gitcode.com/GitHub_Trending/cu/cuda-samples