single float sparse matrix in mex files- using cusparse

2018 7 26
2 个回答
回答已采纳
更新时间：2018 8 16
2 次查看（30 天）

0 个投票

Hi,
As you know "mxCreateSparse" return a double sparse matrix instead of single matrix, and "mxCreateSparseNumericMatrix", is not available there any more (on R2018a). I did try to compile this (mxCreateSparseNumericMatrix) on R2018a, and after adding library file(libmx.lib from R2013a) I could compile the mex file including "mxCreateSparseNumericMatrix" on R2018a but I could not run it despite adding dll library files, so I would like to know is there any suggestion I could use to successful compile and run this on R2018a, or maybe you can create and send me some file that can be compile and run on R2018a, I truly appreciate your time and effort. I have implemented all the available operation on double sparse matrix using CUDA library(such as matrix matrix addition, multiplication etc )and I could use this for single sparse matrix.
Thank you very much.
Best regards, Kevin

1 个评论
显示 -1更早的评论 隐藏 -1更早的评论

James Tursa
James Tursa 2018-7-26
编辑：James Tursa 2018-7-26
Can you post a short example where mxCreateSparseNumericMatrix used to work on older versions? What were you doing with the result of this call in MATLAB?

请先登录，再进行评论。

请先登录，再回答此问题。

 采纳的回答

Kevin Bronik
Kevin Bronik 2018-7-27

0 个投票

Hi, Here the code:
// floatsparse.c
#include < math.h >
#include "mex.h"
#include < string.h >
#include "matrix.h"
EXTERN_C mxArray *mxCreateSparseNumericMatrix(mwSize m, mwSize n, mwSize nzmax, mxClassID classid, mxComplexity ComplexFlag);
void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] ) {
mwSize m,n;
mwSize i,j,k;
mwSize nzmax, nnz;
double *pri;
double *irx;
double *jcx;
float *pr;
mwIndex *ir;
mwIndex *jc;
size_t buflenp, buflenpi, buflenpj;
m = (mwSize)mxGetScalar(prhs[3]);
n = (mwSize)mxGetScalar(prhs[4]);
nnz =(mwSize)mxGetScalar(prhs[5]);
nzmax=(mwSize)mxGetScalar(prhs[6]);
pri = mxGetPr(prhs[0]);
irx = mxGetPr(prhs[1]);
jcx = mxGetPr(prhs[2]);
buflenp = nnz*sizeof(float);
pr = mxMalloc(buflenp);
buflenpi = nnz*sizeof(mwIndex);
ir = mxMalloc(buflenpi);
buflenpj = (n+1)*sizeof(mwIndex);
jc = mxMalloc(buflenpj);
for (i=0 ; i<nnz; i++){
pr[i]=(float) (pri[i]);
ir[i]=(mwIndex) (irx[i]);
}
for (j=0 ; j<n+1; j++){
jc[j]=(mwIndex) (jcx[j]);
}
plhs[0] = mxCreateSparseNumericMatrix(m, n, nzmax, mxSINGLE_CLASS, mxREAL);
memcpy((void*)mxGetPr(plhs[0]), (const void*)pr, (nnz)*sizeof(float));
memcpy((void*)mxGetIr(plhs[0]), (const void*)ir, (nnz)*sizeof(mwIndex));
memcpy((void*)mxGetJc(plhs[0]), (const void*)jc, (n+1)*sizeof(mwIndex));
}
Compiling [on a window 10 machine (MATLAB, R2013a)]:
mex -largeArrayDims floatsparse.c
Testing:
>> val=[1.0, 7.0, 5.0, 3.0, 4.0, 2.0, 6.0]
val =
1 7 5 3 4 2 6
>> ir = [0, 2, 4, 2, 3, 0, 4]
ir =
0 2 4 2 3 0 4
>> jc = [0, 3, 5, 5, 7]
jc =
0 3 5 5 7
>> f=floatsparse(val, ir, jc, 5, 4, 7 ,10)
f =
(1,1) 1
(3,1) 7
(5,1) 5
(3,2) 3
(4,2) 4
(1,4) 2
(5,4) 6
>> whos f
Name Size Bytes Class Attributes
f 5x4 160 single sparse
I really appreciate your time and help.
Kevin

3 个评论
显示 1更早的评论 隐藏 1更早的评论

James Tursa
James Tursa 2018-8-3
编辑：James Tursa 2018-8-3
Preliminary investigation results:
mxCreateSparseNumericMatrix existed as a C routine in the libmx library up through R2013b. Starting in R2014a it appears that this routine only exists as a C++ routine. So, presumably if you had the correct prototype defined you might be able to link with it and use it if you use a C++ mex routine.
That being said, to my knowledge there is no supporting functionality for single sparse matrices in MATLAB. So although you could create such a variable with your mex routine, you couldn't use it at the m-code level for any calculations. E.g.,
>> mex -largeArrayDims floatsparse.c
>> val=[1.0, 7.0, 5.0, 3.0, 4.0, 2.0, 6.0];
>> ir = [0, 2, 4, 2, 3, 0, 4];
>> jc = [0, 3, 5, 5, 7];
>> f=floatsparse(val, ir, jc, 5, 4, 7 ,10)
f =
(1,1) 1
(3,1) 7
(5,1) 5
(3,2) 3
(4,2) 4
(1,4) 2
(5,4) 6
>> whos f
Name Size Bytes Class Attributes
f 5x4 100 single sparse
>> f+5
Undefined function 'plus' for input arguments of type 'single' and attributes 'sparse 2d
real'.
>> f>0
Undefined function 'gt' for input arguments of type 'single' and attributes 'sparse 2d
real'.
Are you just creating these single sparse variables and then passing them off to your CUDA code for computation? If so, you might try the C++ route. If that doesn't work, then it would probably not be too difficult to write a simple replacement mxCreateSparseSingle routine from scratch as long as you didn't need to do any index sorting. (The index sorting could be done, of course, but the routine is no longer "simple" to write). E.g., create a double sparse matrix first and then hack into it to turn it into a single class variable by manipulating the classid etc.
Kevin Bronik
Kevin Bronik 2018-8-13
Hi, Yes I am creating these single sparse variables and then passing them off (inside) to my CUDA code [I have implemented this using double sparse (all the functions such as 'plus' etc. using cusparse and cusolver library) and would like to implement all using single sparse]. I did try all your suggestions and it didn't work so I would like to ask you if possible! you create for me a simple sparse single (a prototype using mex routine). I will definitely mention it (a collaborative work with MathWorks team) in any my future implementation related to single sparse. [After my knowledge, I am not the only one and there is a lot of student/researcher who could be potentially interested in trying single sparse applications, specially on GPU], thanks.
James Tursa
James Tursa 2018-8-13
编辑：James Tursa 2018-8-13
I assume you want only the R2017b memory model, with real & imag parts separate? What version of MATLAB are you using? Complex mxArray variables only exist as iterleaved real/imag data in R2018a and later. There is no separate real/imag mxArray version available. So if you will be working with complex sparse single matrices that need separate real/imag parts I don't even know how to create such a variable in R2018a. The internal mxArray format simply doesn't support it. Will you be working with real variables only? Or also complex? How do you pass this variable to CUDA? Do you use mxGetPr, mxGetM, mxGetN, mxGetIr, mxGetJc to get pointers etc and then pass these off to CUDA individually? I.e., your CUDA code doesn't know anything about mxArray stuff but simply gets sizes and pointers to the data? Or ...?
If you want to use complex single sparse matrices in R2018a or later but still use the separate real/imag data model, you would have to bundle the real/imag data into one data block, put it in the real spot, and then write your own mxGetImagData function to pick off the imag data pointer to pass to your CUDA routines.

请先登录，再进行评论。

更多回答（1 个）

Kevin Bronik
Kevin Bronik 2018-8-16

0 个投票

I am using MATLAB-R2018a, I would work with both real and complex [actually I did implement both (real and complex double precision) using cusparse and cusolver library] For Example the generic CUDA-OMP-MEX model I use:
......
size_t pivot_dimensions[2] = {numARows, numAColumns};
mxGPUArray *INP = mxGPUCreateGPUArray(2, (mwSize*) pivot_dimensions, mxDOUBLE_CLASS,
mxCOMPLEX, MX_GPU_DO_NOT_INITIALIZE);
cuDoubleComplex *d_A_ = (cuDoubleComplex *)mxGPUGetData(INP);
......
mxArray * ROW =mxCreateNumericMatrix(nR, 1, mxINT32_CLASS, mxREAL);
int *RowIndices = (int *)mxGetInt32s(ROW);
mxArray * COLUMN =mxCreateNumericMatrix(nC, 1, mxINT32_CLASS, mxREAL);
int *ColIndices = (int *)mxGetInt32s(COLUMN);
mxArray * VALUE =mxCreateNumericMatrix(nnz, 1, mxDOUBLE_CLASS, mxREAL);
double *Val = (double *)mxGetDoubles(VALUE);
......
mwIndex *irs,*jcs;
irs = static_cast<mwIndex *> (mxMalloc (nR * sizeof(mwIndex)));
#pragma omp parallel for shared(nR) private(i)
for ( i = 0; i < nR; ++i) {
irs[i] = static_cast<mwIndex> (Cu_Row[i])-1;
}
jcs = static_cast<mwIndex *> (mxMalloc ((NcolsA+1) * sizeof(mwIndex)));
int nc1= NcolsA+1;
#pragma omp parallel for shared(nc1) private(i)
for (i = 0; i < nc1; ++i) {
jcs[i] = static_cast<mwIndex> (Cu_Col[i])-1;
}
OUTPUTMATRIX = mxCreateSparse(NrowsA,NcolsA,nnzm,mxREAL);
......
mxFree (mxGetJc (OUTPUTMATRIX)) ;
mxFree (mxGetIr (OUTPUTMATRIX)) ;
mxFree (mxGetDoubles (OUTPUTMATRIX)) ;
mxSetNzmax(OUTPUTMATRIX, (static_cast<mwSize>(NNZMAXA)));
mxSetIr(OUTPUTMATRIX, (mwIndex *)irs);
mxSetJc(OUTPUTMATRIX, (mwIndex *)jcs);
int s = mxSetDoubles(OUTPUTMATRIX, (mxDouble *)PRS);
......
Yes, my CUDA code doesn't know anything about mxArray stuff and simply after creating mxArray/GPUArray using MATLAB C-API, it uses again exiting MATLAB function/pointers to access/manipulate the data (set/get the data etc). thanks

类别

在 帮助中心 和 File Exchange 中查找有关 GPU Computing 的更多信息

产品

版本

R2018a

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by