在图像处理中,常常会需要求解超大稀疏矩阵的解,比如泊松融合,比如加权最小二乘法的保边缘平滑滤波器(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的方法