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

Question

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

• 对于大方程系统(实际上从只有100个元素开始)是非常不稳定的，
• 它比Matlab CPU替代品慢几十倍或数百倍。

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

硬件规格:Intel i7 3630QM, GT640M 2G, 8gb 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 FilesMATLABR2012bexterninclude" -I "C:Program Files (x86)Microsoft Visual Studio 10.0VCinclude
>> mex ex.obj -L"C:Program FilesNVIDIA GPU Computing ToolkitCUDAv5.0libx64" -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

where times -单独的执行时间，where last element -用于解决方案的迭代计数(bicstab监视器)，result - K*x= f的解决方案

结果(http://www.failai.lt/rupaliln7kfb/results.7z.htm):

• K_int_6, F_int_6 - ok
• K_11, F_11 - x(1)错误，其他正常
• K_100000, F_100000 - x(1)错误，其他从一开始是好的，但后来与正确结果相比减少。
• K_100000, F_100000 -执行时间在GPU (MEX)上为0.6 s，在CPU (tic;xcpu=KF;toc;)上为0.014 s。

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

Answer 1

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

稀疏矩阵向量乘法重载CUSP