利用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的方法