利用cuda的cusparse模块计算超大型稀疏矩阵方程的解

在图像处理中,常常会需要求解超大稀疏矩阵的解,比如泊松融合,比如加权最小二乘法的保边缘平滑滤波器(WLS).对于一个300*300的图像输入,其对应的稀疏矩阵大小为90000*90000,这样的矩阵别说是求解,存储在内存中都是一笔巨大的开销.好在我们有专门用于稀疏矩阵运算的工具,比如eigen矩阵运算库和cuda中的cusolver模块.

在使用cusolver的过程中,我发现网上关于此的教程是少之又少,因而用了不少时间来阅读官方文档,和理解api,为了方便中文社区的使用者们,这里写下关键步骤的教程,以供参考.

首先,不管是cuda还是eigen还是matlab,对于稀疏矩阵的存储都不是把每一个元素存储下来,最基本的稀疏矩阵存储方式是coo,也就是只存储非零元素的值以及坐标.cusolver模块中往往需要CSR格式的稀疏矩阵,CSR是基于COO的更进一步的压缩格式,它对于值的row坐标进行了压缩.这里我推荐大家去看CUSPARSE(cuda的稀疏矩阵模块)的文档的3.3节中介绍coo以及csr的部分,比较清晰.CSR难理解一点,推荐大家去搜索"csr 稀疏矩阵"这些关键词再看看csr的row到底是怎么压缩的.

第二步,大家去看安装cuda时产生的samples的例子,一般在安装cuda时会默认安装开发工具nsight,一个基于eclipse的IDE,用起来很方便,直接新建项目--导入cuda项目--导入cuda Sample就可以了.


其中NVIDIA_CUDA-8.0_Samples/7_CUDALibraries 这个目录下有一些cusolver模块的例子,下面这几个可以运行一下看一下,帮助大家理解.

cuSolverDn_LinearSolver - cuSolverDn Linear Solver 
实现cuSolverDN的LU, QR和Cholesky因式分解
cuSolverRf - cuSolverRf Refactorization 
重新因子化。
cuSolverSp_LinearSolver - cuSolverSp Linear Solver 
实现cuSolverSP的LU, QR和Cholesky因式分解
cuSolverSp_LowlevelCholesky - cuSolverSp LowlevelCholesky Solver 
使用cuSolverSP底层API实现Cholesky因式分解
cuSolverSp_LowlevelQR - cuSolverSp Lowlevel QR Solver 
使用cuSolverSP底层API实现QR因式分解

接下来,我们创建一个coo格式的稀疏矩阵,很简单,一个coo格式的矩阵只需要3个变量:非零元素的值,非零元素的横坐标,非零元素的纵坐标,我们是预先知道这些元素的个数的,所以创建然后分配大小.

int *h_cooRowIndA = NULL; //存储非零元素的纵坐标
int *h_coo_csrColIndA = NULL;  //存储非零元素的横坐标
float *h_coo_csrValA = NULL;  //存储非零元素的值
h_coo_csrValA = (float *)malloc(sizeof(float)*n*5);  //分配n*5个元素,这也是我们的稀疏矩阵的非零元素的个数
h_cooRowIndA = (int *)malloc(sizeof(int)*n*5); 
h_coo_csrColIndA = (int *)malloc(sizeof(int)*n*5); 

解释下我的数据,是带状稀疏矩阵,总共有5个带,所以我的非零元素有 n*5 个


然后,我们需要对分配的内存空间进行赋值了,当然了这个值根据你需要解决的问题来决定.但是赋值方法可以共享给大家,我这里是这样赋值的


下面就是coo转换为csr格式了,我提供一个函数给大家(当然cusparse模块有coo2csr的函数,不过那是一个global函数,也就是在显卡上操作数据的函数,需要对cudaMalloc的device数据进行操作,如果是在host分配的数据,只能自己写函数)

void DeepAnalogy::compress_index(
    const int *Ind, //coo格式的row坐标
    int nnz,  //矩阵总共有多少非零元素?
    int m, //矩阵每一行有几个元素?
    int *Ptr, //返回的csr格式的row坐标
    int base) //下标从0开始还是从1开始?
{
    int i;

    /* initialize everything to zero */
    for(i=0; i<m+1; i++){
        Ptr[i]=0;
    }
    /* count elements in every row */
    Ptr[0]=base;
    for(i=0; i<nnz; i++){
        Ptr[Ind[i]+(1-base)]++;
    }
    /* add all the values */
    for(i=0; i<m; i++){
        Ptr[i+1]+=Ptr[i];
    }
}

一个使用的范例如下,这里h_csrRowPtrA需要预先分配好空间,大小就是矩阵的行数+1

h_csrRowPtrA = (int *)malloc(sizeof(int)*(n+1));
compress_index(h_cooRowIndA, NonZero_count, n, h_csrRowPtrA, CUSPARSE_INDEX_BASE_ZERO);
下面,我们将会通过chol分解方式求解稀疏矩阵A*x=B,也可以使用LU或者QR分解,请看samples里的例子,差别不大.

目前位置,A*x=B方程我们分配好了A,接下来我们为x分配空间,并且输入B的值.

float *h_bs1 = NULL, *h_xs2 = NULL;
h_bs1 = (float*)malloc(sizeof(float)*n);
.....按照自己的需要构建B的值,此部分因人而异,省略.....
h_xs1 = (float*)malloc(sizeof(float)*n);

用cpu求解很简单(没错,就是cpu,cuda提供了cpu求解的函数;在矩阵不大的情况下cpu求解反而比gpu求解更快一点),一些必要的初始化,然后一行语句求解:

cusolverSpHandle_t handle = NULL;
cusparseHandle_t cusparseHandle = NULL;
cudaStream_t stream = NULL;
checkCudaErrors(cusolverSpCreate(&handle));
checkCudaErrors(cusparseCreate(&cusparseHandle));
checkCudaErrors(cudaStreamCreate(&stream));
checkCudaErrors(cusolverSpSetStream(handle, stream));
checkCudaErrors(cusparseSetStream(cusparseHandle, stream));
checkCudaErrors(cusolverSpScsrlsvcholHost(
	handle, n, NonZero_count,
	descrA, h_coo_csrValA, h_csrRowPtrA, h_coo_csrColIndA,
	h_bs1, tol, reorder, h_xs1, &singularity));

用gpu的话,就麻烦一点,需要拷贝数据,然后求解,最后拷贝回来.

以上所说的求解是只求解一次,有时候我们有多个b值,需要求解多次A*x=b,这时候每次都要分解A就不应该了.此时我们可以先分解A

//compute A = L*L^T
checkCudaErrors(cusolverSpScsrcholFactor(handle, n, NonZero_count,
	descrA, d_csrValA, d_csrRowPtrA, d_csrColIndA, d_info, buffer_gpu));

然后对多个b进行求解

//solve A*x = b
checkCudaErrors(cusolverSpScsrcholSolve(handle, n, d_b1, d_x1, d_info, buffer_gpu));
checkCudaErrors(cusolverSpScsrcholSolve(handle, n, d_b2, d_x2, d_info, buffer_gpu));
checkCudaErrors(cusolverSpScsrcholSolve(handle, n, d_b3, d_x3, d_info, buffer_gpu));

这样可以减少不必要的计算,samples里面的cuSolverSp_LowlevelCholesky介绍了这部分代码.

展示下各个部分运行时间,单位是秒


最后附上eigen求解A*x = b的方法


  • 5
    点赞
  • 38
    收藏
    觉得还不错? 一键收藏
  • 4
    评论
CUDA中,稀疏矩阵和稠密矩阵相乘是一个非常常见的操作。稀疏矩阵是指矩阵中大部分元素为零,而稠密矩阵则是指矩阵中大部分元素都非零。 在处理稀疏矩阵和稠密矩阵相乘时,通常需要进行以下几个步骤: 1. 稠密矩阵数据传输:将稠密矩阵数据从主机内存传输到GPU设备的全局内存中,以便后续在GPU上进行计算。 2. 稀疏矩阵数据结构转换:将稀疏矩阵由传统的行压缩存储(CSR)或列压缩存储(CSC)格式转换为适合在GPU上处理的稀疏格式,例如ELL格式(ELLPACK/ITPACK)或使用了线程合并和合并访问的CSR格式(CSR-TC)。 3. 稀疏矩阵和稠密矩阵相乘的计算:通过使用CUDA的并行计算特性,在GPU上进行稀疏矩阵和稠密矩阵的并发计算,以提高计算效率。在此过程中,我们通常会使用CUDA的线程、块和网格来处理数据并进行并行计算。 4. 结果数据传输:将计算得到的结果从GPU全局内存中传输回主机内存,以供后续的数据处理或输出。 需要注意的是,稀疏矩阵和稠密矩阵相乘的计算方法可能因具体情况而异,选择适合的算法和数据结构能够提高计算性能。此外,在实际应用中,还可以采用一些优化技术,如共享内存的使用、存储器访问模式的优化等,以进一步提高计算效率。 通过使用CUDA并行计算的能力,我们可以有效地进行稀疏矩阵和稠密矩阵的相乘操作,从而提高计算效率,并在处理大规模数据时节省时间和资源。
评论 4
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值