改进 CUDA GPU 上 A*x=B 的 Matlab + CUSP MEX 解决方案

Question

Matlab 仍然无法在 CUDA GPU 上计算稀疏矩阵。也没有这样的工具箱（夹克已停产）。这就是为什么我使用通过 MEX 文件集成到 Matlab 的 CUSP。但是，我开发的工具有两个问题：

• 对于大型方程组（实际上仅从 100 个元素开始），它非常不稳定，
• 它比 Matlab CPU 替代品慢几十或几百倍。

我正在求解 A*x=b，其中 A 是一个稀疏的对称矩阵，b 是一个向量。

硬件规格：Intel i7 3630QM、GT640M 2G、8 GB DDR3。 软件：Windows 8 64 位、Matlab R2012b 64 位、CUDA 5.0 64 位、CUSP 0.3.1、Windows SDK v7.0、VS2010 编译器。

墨西哥代码：

#include<cusp/csr_matrix.h>
#include <cusp/krylov/bicgstab.h>
#include <matrix.h>
#include <mex.h> 
#include <time.h>

void mexFunction(int nlhs,mxArray *plhs[],int nrhs,const mxArray *prhs[])
{
        double t1 =  clock();
          // data from Matlab       
        double *b = mxGetPr(prhs[1]);
        double *A = mxGetPr(prhs[0]);
        int n = mxGetM(prhs[0]);
        mwIndex *ir = mxGetIr(prhs[0]);
        mwIndex *jc = mxGetJc(prhs[0]);
        int N = jc[n];
        t1 = clock() - t1;

        double t2 =  clock();
          // initialization of matrix A in CSR format (jc and ir are exchanged, because Matlab uses CSC format
        cusp::csr_matrix<int,float,cusp::device_memory> Ag(n,n,3*n-2);
        thrust::copy(jc, jc + n + 1, Ag.row_offsets.begin());
        thrust::copy(ir, ir + N,     Ag.column_indices.begin());
        thrust::copy(A,  A  + N,     Ag.values.begin()); 
          // initialization of vector b
        cusp::array1d<float, cusp::device_memory> bg (b, b+n);
        cusp::array1d<float, cusp::device_memory> xg (n, 0);
        t2 = clock() - t2;

        double t3 =  clock();
          // bicgstab algorithm solution for vector x, when using 0.001 accuracy and precondition M
          // this is the slowest part, much slower than others
        cusp::verbose_monitor<float> monitor(bg, 5000, 1e-3);
        cusp::identity_operator<float, cusp::device_memory> M(n, n);
        cusp::krylov::bicgstab(Ag, xg, bg, monitor, M);        
        t3 = clock() - t3;

        double t4 =  clock();     
          // gathering solution vector bact on host to Matlab array T
        mxArray *T = mxCreateDoubleMatrix(n, 1, mxREAL);
        double *x  = mxGetPr(T);
        thrust::copy(xg.begin(), xg.end(), x);
        t4 = clock() - t4;

          // gathering execution times to Matlab array times
        mxArray *times=mxCreateDoubleMatrix(5, 1, mxREAL);
        double *timesb=mxGetPr(times);
        timesb[0]=t1; timesb[1]=t2; timesb[2]=t3; timesb[3]=t4; timesb[4]=monitor.iteration_count();

          // sending data back to Matlab
        plhs[0] = times; 
        plhs[1] = T;
} 

使用这些命令在 Matlab 上的 MEX 文件 (ex.cu) 中编译此代码（如有必要，将第二个命令更改为 32 位）：

>> !nvcc -c -arch sm_20 ex.cu -Xcompiler -fPIC -I "C:\Program Files\MATLAB\R2012b\extern\include" -I "C:\Program Files (x86)\Microsoft Visual Studio 10.0\VC\include
>> mex ex.obj -L"C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v5.0\lib\x64" -lcudart

示例矩阵、向量和编译的 64 位 MEX 函数：http://www.failai.lt/3fqkhvoslxyt/sampleData.7z.htm

用途：

tic; [times,x]=ex(K',F); toc;   %K has to be transposed for CSR

其中时间 - 单独的执行时间，其中最后一个元素 - 用于解决方案的迭代计数（bicgstab 监视器），结果 - K*x=F 的解决方案。

结果（http://www.failai.lt/rupaliln7kfb/results.7z.htm）：

• K_int_6，F_int_6 - 好的
• K_11, F_11 - x(1) 错误，其他正常
• K_100000, F_100000 - x(1) 错误，其他从一开始还可以，但后来与正确结果相比正在减少。
• K_100000, F_100000 - 在 GPU (MEX) 上执行持续 0.6 秒，而在 CPU 上执行持续 0.014 秒 (tic;xcpu=K\F;toc;)。

你能看看那个代码，也许试试 MEX 函数，报告你的结果，建议如何改进这个函数？也许您知道任何可以在 GPU 上进行稀疏计算的替代方案？我希望，在 Matlab 发布对 GPU 上稀疏矩阵的兼容性之前，它对每个人都有用 :)

Answer 1

看看 Matlab 文件交换，gpus 的尖点稀疏类，支持单精度，实数/复数： http://www.mathworks.com/matlabcentral/fileexchange/44423-gpu-sparse-accumarray-non-uniform-grid

稀疏矩阵向量乘法被 CUSP 重载。