NIVIDIA 高性能计算CUDA笔记（五） cuSOLVER库的简介与矩阵特征值分解示例

啸妹回 · 5 小时前

NIVIDIA 高性能计算CUDA笔记（五）

cuSOLVER库的简介与矩阵特征值分解示例

cuSOLVER是NVIDIA提供基于cuBLAS和cuSPARSE库的GPU加速的线性代数操作算子库，其专注于稠密和稀疏矩阵的高级线性代数运算，，包括矩阵分解和求解方程等内容。其里面包括了三个独立库：\(cuSOlverDN、cuSolver、cuSolverRF\) 三部分。其特点为：

高性能：利用GPU并行计算能力，显著加速线性代数运算；
兼容性：支持单精度（float）、双精度（double）、复数等数据类型；
易用性：提供类似LAPACK的API，便于传统CPU代码迁移到GPU；
与cuBLAS集成：底层依赖cuBLAS进行基本线性代数运算；

cuSolver的典型应用场景：（1）科学计算（如计算流体力学、结构分析）；（2）机器学习（如PCA、线性回归）；（3）信号处理；（4）计算视觉等。 cuSOLVER与cuBLAS的区别：

cuBLAS主要提供基础的线性代数运算（如矩阵乘法，向量操作）；
cuSolver专注于高级线性代数求解（如方程组求解，矩阵分解）；

1.cuSolver库主要功能API说明

cuSOLVER主要分为三个部分及部分通用函数接口：

cuSolverDN （Dense Linear Algebra）用于稠密矩阵的运算，cuSolverSP库主要设计用于求解稠密线性方程组：

\[\mathbf{Ax}=\mathbf{b} \tag{1}\]
其中系数矩阵、右侧向量和解向量：\(\mathbf{A}\in \mathbf{R}_{n\times{n}},\mathbf{b}\in\mathbf{R}_{n},\mathbf{x}\in{\mathbf{R}_{n}}\)。
\(cuSolverDN\) 库提供了\(QR\)分解和\(LU\) 分解，以处理一般的矩阵，该矩阵可能是非对称的。对于对称或Hermit矩阵提供了Cholesky分解。对于对称不定的矩阵，提供了LDL分解。
其主要功能包括：
功能关键函数（前缀为cuSolverDN）说明矩阵分解[s/d/c/z]getrfLU分解[s/d/c/z]potrfCholesky分解[s/d/c/z]getrsQR分解线性方程组求解[s/d/c/z]getrs用LU分解结果求解[s/d/c/z]ports用Cholesky分解后求解[s/d/c/z]gels最小二乘问题（超定或欠定方程组）特征值问题[s/d]syevd对称矩阵特征值分解（分治算法）[c/z]heevd埃尔米特矩阵特征值分解奇异值分解[s/d/c/z]gesvd奇异值分解（SVD）句柄操作cusolverDnCreate()创建句柄cusolverDnDestroy()销毁句柄

cuSolverSP（Sparse Linear Algebra）用于稀疏矩阵运算，为稀疏矩阵提供了一个新的工具箱，用于解决稀疏线性系统，包括\(QR\)分解，由于不是所有稀疏模式都有性能良好的分解算法模式，所以\(cuSolver\)还提供一些CPU执行的方法。
cuSolver库主要设计用于求解稀疏线性系统

\[\mathbf{A}\mathbf{x}=\mathbf{b} \tag{2}\]
以及稀疏最小二乘问题：

\[\mathbf{x}=argmin||\mathbf{A}\mathbf{x}-\mathbf{b}|| \tag{3}\]
其中\(\mathbf{A}\in{R^{m\times{n}}}\)稀疏矩阵,右侧输出向量\(b\in{R^{m}}\) 和解向量\(x\in{R^{n}}\)。其核心算法是基于稀疏矩阵的编码分解。矩阵以\(CSR\) 格式被接受。如果矩阵是对称的或Hermit矩阵，因此用户必须提供完整的矩阵，即填充缺失的下部分或上部。如果矩阵是对称正定的，用户只需要求解，则需要Cholesky分解可行，用户只需要提供的下三角部分。除了线性和最小二乘求解器，cuSolver库还提供了一个幂乘法的简单特征求解器。
功能关键函数（前缀为cusolverSP）说明稀疏线性方程求解（CPU分析）csrlsvlu使用LU分解求解（非对称矩阵）csrlsvqr使用QR分分解求解（任意矩阵）csrlsvchol使用Cholesky分解求解（对称正定）稀疏最小二乘法问题csrlsqvqr稀疏最小二乘法（QR分解）稀疏矩阵特征值分解[s/d]csreigvsi稀疏对称矩阵特征值（迭代）句柄管理cusolverSpCreate()创建句柄cusolverSpDestroy()销毁句柄
cuSolverRF(Refactorization) 用于矩阵分解，加速解决一系列具有相同稀疏模式但不同数值的线性方程组。
\(cuSolverRF\) 库旨在通过快速重分解加速线性系统集合的求解，当给定相同稀疏度模式的新系数时

\[\mathbf{A}_i\mathbf{x}_i=f_i \tag{4}\]
其中给出了一组系数矩阵、右侧的观测向量和解向量：\(\mathbf{A}_i\in{\mathbf{R}^{n\times{n}}},\mathbf{f}_i\in\mathbf{R}^n,\mathbf{x}_i\in{R^n}(i=1,...,k)\)
功能关键函数（前缀为cusolverRf）说明重分解求解cusolverRfCreate()创建句柄cusolverRfSetup[Host/Device]()设置矩阵cusolverRfAnalyze()分析矩阵的结构cusolverRfRefactor()数值重分解cusolverRfSolve()线性方程矩阵求解cusolverRfDestroy()销毁句柄
通用接口及辅助功能主要接口：

功能关键函数的接口说明错误处理cusolverGetStatusString()获取错误信息设备管理cusolverGetProperty()获取库版本号内存管理cusolverDn/SpSetStream()设置CUDA流cusolverDn/SpGetStream()获取CUDA流缓存管理cusolverDn/SpSetWorkspace()设置工作空间cusolverDn/Sp[Device/Host]BufferSize()计算所需缓冲区大小注意： 数据类型前缀说明：
数据类型的前缀数据类型的说明\(s\)单精度浮点(float)\(d\)双精度浮点(double)\(c\)单精度复数(cuComplex)z双精度复数(cuDoubleComplex)主要数据结构：
句柄数据结构功能描述cusolverDnHandle_tcuSolverDN句柄（上下文）cusolverSpHandle_tcuSolverSP 句柄（上下文）cusolverRfHandle_tcuSolverRF 句柄（上下文）csr[lu/qr/chol]Info_t稀疏分解信息结构体2.cusolver的特征值分解示例

cuSolver的特征值分解的流程：

创建\(cuSolverDN\)句柄；
在设备上分配内存，并将矩阵A和数据传输到设备；
准备workspace，并计算所需workspace大小；
调用syevd函数进行计算；
将结果（特征值、特征向量）传输回主机；
释放设备内存和句柄；

注意：矩阵A是列优先存储（和Fortran一样，即列主序）。在C/C++中我们通常按行主序存储，因此需要注意转换。
下面以双精度实数对称矩阵为例，使用cusolverDnDsyevd函数。
函数矩阵类型算法特点适用场景syevd实对称/复埃尔米特分治法稳定，内存占用较大中小型矩阵，需要高精度syevdx实对称/复埃尔米特二分法+逆迭代可计算特征值子集只需要部分特征值syevj实对称/复埃尔米特Jacobi法高精度，可并行性好需要高精度特征向量syevjBatched实对称/复埃尔米特Jacobi法批量处理多个小矩阵同时计算gesvd一般矩阵SVD分解计算奇异值和向量奇异值分解问题2.1实数对称矩阵的特征值分解

#include <stdio.h>
#include <stdlib.h>
#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <cusolverDn.h>
int main()
{
cusolverDnHandle_t cusolverH = NULL;
const int m = 3;
const int lda = m;
float *A = (float*)malloc(lda*m * sizeof(float));
A[0] = 3.5;
A[1] = 0.5;
A[2] = 0;
A[3] = 0.5;
A[4] = 3.5;
A[5] = 0;
A[6] = 0;
A[7] = 0;
A[8] = 2;
float W[m]; // eigenvalues最终保存结果
int info_gpu = 0;//计算状态保存
// 步骤1：声明句柄
cusolverDnCreate(&cusolverH);
// 步骤2：分配显存空间
float *d_A = NULL; cudaMalloc((void**)&d_A, sizeof(float) * lda * m);//声明Hermite矩阵（与计算后的特征向量为同一空间）
float *d_W = NULL; cudaMalloc((void**)&d_W, sizeof(float) * m);//声明特征值存储空间
int *devInfo = NULL; cudaMalloc((void**)&devInfo, sizeof(int));//声明计算结果状态空间
cudaMemcpy(d_A, A, sizeof(float) * lda * m, cudaMemcpyHostToDevice);//数据拷贝
// 步骤3：申请计算缓存空间，并在显存中申请该空间
float *d_work = NULL;
int lwork = 0;
cusolverEigMode_t jobz = CUSOLVER_EIG_MODE_VECTOR; // compute eigenvalues and eigenvectors.
cublasFillMode_t uplo = CUBLAS_FILL_MODE_LOWER;
cusolverDnSsyevd_bufferSize(cusolverH, jobz, uplo, m, d_A, lda, d_W, &lwork);//计算evd计算所需存储空间,保存到lwork中
cudaMalloc((void**)&d_work, sizeof(float)*lwork);
// 步骤4：特征分解
cusolverDnSsyevd(cusolverH, jobz, uplo, m, d_A, lda, d_W, d_work, lwork, devInfo);
cudaDeviceSynchronize();
//步骤5：数据读回
cudaMemcpy(A, d_A, sizeof(float)*lda*m, cudaMemcpyDeviceToHost);
cudaMemcpy(W, d_W, sizeof(float)*m, cudaMemcpyDeviceToHost);
cudaMemcpy(&info_gpu, devInfo, sizeof(int), cudaMemcpyDeviceToHost);
printf("%d\n", info_gpu);
printf("eigenvalue = (matlab base-1), ascending order\n");
for (int i = 0; i < m; i++) {
printf("W[%d] = %E\n", i + 1, W[i]);
}
for (size_t i = 0; i < m; i++)
{
for (size_t j = 0; j < m; j++)
{
printf("%.4f\n", A[i + j*m]);
}
}
return 0;
}

复制代码

2.2 复数Hermit矩阵的双精度特征分解的案例

#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <stdio.h>
#include <stdlib.h>
#include <cusolverDn.h>
int main()
{
cusolverDnHandle_t cusolverH = NULL;
const int m = 4;
const int lda = m;
cuDoubleComplex *A = (cuDoubleComplex*)malloc(lda*m*sizeof(cuDoubleComplex));
A[0] = make_cuDoubleComplex(1.9501e2,0);
A[1] = make_cuDoubleComplex(0.2049e2, 0.1811e2);
A[2] = make_cuDoubleComplex(0.5217e2, 0.3123e2);
A[3] = make_cuDoubleComplex(0.2681e2, 0.3998e2);
A[4] = make_cuDoubleComplex(0.2049e2, -0.1811e2);
A[5] = make_cuDoubleComplex(1.8272e2, 0);
A[6] = make_cuDoubleComplex(0.5115e2, -0.0987e2);
A[7] = make_cuDoubleComplex(0.4155e2, -0.0435e2);
A[8] = make_cuDoubleComplex(0.5217e2, -0.3123e2);
A[9] = make_cuDoubleComplex(0.5115e2, 0.0987e2);
A[10] = make_cuDoubleComplex(2.3984e2, 0);
A[11] = make_cuDoubleComplex(0.4549e2, -0.0510e2);
A[12] = make_cuDoubleComplex(0.2681e2, -0.3998e2);
A[13] = make_cuDoubleComplex(0.4155e2, 0.0435e2);
A[14] = make_cuDoubleComplex(0.4549e2, 0.0510e2);
A[15] = make_cuDoubleComplex(2.2332e2, 0);
//1.9501 + 0.0000i 0.2049 - 0.1811i 0.5217 - 0.3123i 0.2681 - 0.3998i
// 0.2049 + 0.1811i 1.8272 + 0.0000i 0.5115 + 0.0987i 0.4155 + 0.0435i
// 0.5217 + 0.3123i 0.5115 - 0.0987i 2.3984 + 0.0000i 0.4549 + 0.0510i
// 0.2681 + 0.3998i 0.4155 - 0.0435i 0.4549 - 0.0510i 2.2332 + 0.0000i
double W[m]; // eigenvalues最终保存结果
int info_gpu = 0;//计算状态保存
// step 1: create cusolver/cublas handle
cusolverDnCreate(&cusolverH);
// step 2: copy A and B to device
cuDoubleComplex *d_A = NULL; cudaMalloc((void**)&d_A, sizeof(cuDoubleComplex) * lda * m);//声明Hermite矩阵（与计算后的特征向量为同一空间）
double *d_W = NULL; cudaMalloc((void**)&d_W, sizeof(double) * m);//声明特征值存储空间
int *devInfo = NULL;cudaMalloc((void**)&devInfo, sizeof(int));//声明计算结果状态空间
cudaMemcpy(d_A, A, sizeof(cuDoubleComplex) * lda * m, cudaMemcpyHostToDevice);//数据拷贝
// step 3: query working space of syevd
cuDoubleComplex *d_work = NULL;
int lwork = 0;
cusolverEigMode_t jobz = CUSOLVER_EIG_MODE_VECTOR; // compute eigenvalues and eigenvectors.
cublasFillMode_t uplo = CUBLAS_FILL_MODE_LOWER;
cusolverDnZheevd_bufferSize(cusolverH, jobz, uplo, m, d_A, lda, d_W, &lwork);//计算evd计算所需存储空间,保存到lwork中
cudaMalloc((void**)&d_work, sizeof(cuDoubleComplex)*lwork);
// step 4: compute spectrum
cusolverDnZheevd(cusolverH, jobz, uplo, m, d_A, lda, d_W, d_work, lwork, devInfo);
cudaDeviceSynchronize();
cudaMemcpy(A, d_A, sizeof(cuDoubleComplex)*lda*m,cudaMemcpyDeviceToHost);
cudaMemcpy(W, d_W, sizeof(double)*m, cudaMemcpyDeviceToHost);
cudaMemcpy(&info_gpu, devInfo, sizeof(int), cudaMemcpyDeviceToHost);
printf("%d\n", info_gpu);
printf("eigenvalue = (matlab base-1), ascending order\n");
for (int i = 0; i < m; i++) {
printf("W[%d] = %E\n", i + 1, W[i]);
}
for (size_t i = 0; i < 4; i++)
{
for (size_t j = 0; j < 4; j++)
{
printf("%.4f + %.4f j\n", A[i * 4 + j].x, A[i * 4 + j].y);
}
}
cudaFree(d_A);
cudaFree(d_W);
cudaFree(devInfo);
cudaFree(d_work);
cusolverDnDestroy(cusolverH);
return 0;
}

复制代码

注意事项：

存储顺序：cuSolver 使用列优先存储（Fortran风格）
矩阵对称性：确保输入矩阵满足对称性要求
工作空间：必须先查询正确的工作空间大小
错误检查：总是检查 d_info 和函数返回值
异步执行：大多数 cuSolver 函数是同步的，会阻塞直到完成
流管理：可以使用 cusolverDnSetStream() 将计算绑定到特定 CUDA 流

3.参考资料

官网文档：1. 引言 — cuSOLVER 13.1 文档
参考博客：CUDA加速的线性代数求解器库cuSOLVER-CSDN博客
参考博客：cuSOLVER - 上海交大超算平台用户手册
参考博客：CUDA进阶：cuSOLVER库实战指南 - Dawoai
参考博客：CUDA求解特征值和特征向量使用cusolver库 - 灰信网（软件开发博客聚合）

来源：程序园用户自行投稿发布，如果侵权，请联系站长删除
免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！

账号		自动登录	找回密码
密码			立即注册

NIVIDIA 高性能计算CUDA笔记（五） cuSOLVER库的简介与矩阵特征值分解示例

相关帖子

签约作者

NIVIDIA 高性能计算CUDA笔记 （五） cuSOLVER库的简介与矩阵特征值分解示例

相关帖子

相关推荐

Zsh 配置笔记

从位图结构开始学内存的计算

NIVIDIA高性能计算CUDA笔记（四）cuDNN深度神经网络计算库简介及卷积操作示例

FFmpeg开发笔记（九十七）国产的开源视频剪辑工具AndroidVideoEditor

FFmpeg开发笔记（九十八）基于FFmpeg的跨平台图形用户界面LosslessCut

一个高性能的 .NET MQTT 客户端与服务器库

签约作者

NIVIDIA 高性能计算CUDA笔记（五） cuSOLVER库的简介与矩阵特征值分解示例