uSOLVER库的简介与矩阵特征值分解示例​ 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库主要设计用于求解稠密线性方程组Axb(1)(1)其中系数矩阵、右侧向量和解向量A∈Rn×n,b∈Rn,x∈Rn∈×,∈,∈。​ 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()销毁句柄cuSolverSPSparse Linear Algebra用于稀疏矩阵运算为稀疏矩阵提供了一个新的工具箱用于解决稀疏线性系统包括QR分解由于不是所有稀疏模式都有性能良好的分解算法模式所以cuSolver还提供一些CPU执行的方法。cuSolver库主要设计用于求解稀疏线性系统Axb(2)(2)以及稀疏最小二乘问题xargmin||Ax−b||(3)(3)||−||其中A∈Rm×n∈×稀疏矩阵,右侧输出向量b∈Rm∈ 和解向量x∈Rn∈。其核心算法是基于稀疏矩阵的编码分解。矩阵以CSR 格式被接受。如果矩阵是对称的或Hermit矩阵因此用户必须提供完整的矩阵即填充缺失的下部分或上部。如果矩阵是对称正定的用户只需要求解则需要Cholesky分解可行用户只需要提供的下三角部分。除了线性和最小二乘求解器cuSolver库还提供了一个幂乘法的简单特征求解器。功能关键函数前缀为cusolverSP说明稀疏线性方程求解CPU分析csrlsvlu使用LU分解求解非对称矩阵csrlsvqr使用QR分分解求解任意矩阵csrlsvchol使用Cholesky分解求解对称正定稀疏最小二乘法问题csrlsqvqr稀疏最小二乘法QR分解稀疏矩阵特征值分解[s/d]csreigvsi稀疏对称矩阵特征值迭代句柄管理cusolverSpCreate()创建句柄cusolverSpDestroy()销毁句柄cuSolverRF(Refactorization)用于矩阵分解加速解决一系列具有相同稀疏模式但不同数值的线性方程组。cuSolverRF 库旨在通过快速重分解加速线性系统集合的求解当给定相同稀疏度模式的新系数时Aixifi(4)(4)其中给出了一组系数矩阵、右侧的观测向量和解向量Ai∈Rn×n,fi∈Rn,xi∈Rn(i1,...,k)∈×,∈,∈(1,...,)功能关键函数前缀为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风格矩阵对称性确保输入矩阵满足对称性要求工作空间必须先查询正确的工作空间大小