single float sparse matrix in mex files- using cusparse

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

[James Tursa]

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?

Answer 1

[Kevin Bronik]

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

[James Tursa]

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]

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]

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.

Answer 2

[Kevin Bronik]

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