1

Matlab はまだ CUDA GPU でスパース行列を計算できません。そのためのツールボックスもありません (Jacket は廃止されました)。そのため、MEX ファイルを介して Matlab に統合された CUSP を使用しています。ただし、開発したツールには 2 つの問題があります。

  • 大きな方程式系 (実際には 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 コンパイラ。

MEX コード:

#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 ビットの 2 番目のコマンドを変更します)。

>> !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

ここで、times - 個別の実行時間、last element - 解に使用される反復回数 (bicgstab モニター)、result - 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 上の疎行列の互換性をリリースするまで、すべての人に役立つことを願っています :)

4

1 に答える 1

0

Matlab ファイル交換、gpus の cusp スパース クラス、単精度のサポート、実数/複素数を見てみましょう: http://www.mathworks.com/matlabcentral/fileexchange/44423-gpu-sparse-accumarray-non-uniform-grid

疎行列ベクトル乗算は CUSP でオーバーロードされます。

于 2013-12-02T19:08:04.303 に答える