How does MATLAB internally store sparse arrays / Calling MKL's sparse BLAS library for sparse matrix operations

SHORT VERSION

This is a list of all the sparse formats that MKL supports. Not all of these are supported for the most up-to-date Inspector-executor Sparse BLAS routines, but MKL includes conversion functions so I can make any of these formats work. Does MATLAB use (or have the ability to convert to/from) any of these formats? If so, which ones. If not, how might I convert the format manually and/or is there a better way to call MKL's sparse BLAS routines?

.

.

MORE INFO

I would like to be able to call Intel's MKL library for sparse matrix operations, since apparently sparse matrix operations in MATLAB are still single threaded (I mean, I understand that there isnt enough time to add every new feature, but you'd think this would be a higher priority in a product called "MATrix LABratory"...).

At any rate, Im working on a project where being able to do multi-threaded sparse matrix operations is basically required, so I am trying to figure out the best way to call MKL for sparse operations.

I am in the process of trying to update the MKL library from 2017.0 to 2018.2 using this guide from Intel's website, though I suspect this won't actually work for sparse matrix operations since I'm fairly sure that MATLAB isn't using MKL for these operations (Im fairly sure that MKL 2017.0 supports some form of multi-threaded sparse operations). Even so, I had planned on doing this anyways to unlock the AVX-512 compute abilities in a computer I just finished building, so I figure I mine as well see if I can kill 2 birds with 1 stone.

Assuming updating the MKL library doesn't work, (I think) I'll have to write some sort of wrapper to manually call MKL for sparse operations. I know Matlab provides a few example c/c++ functions for this purpose (that Im hoping won't be too hard to modify, since my c/c++ is iffy at best). However, I don't know if MATLAB stores sparse arrays in a format that MKL supports (and if so which format).

[Jiaen Liu]

I'm running into the same situation. I'm wondering if there is any benchmark about the speed gain of using Intel's multithreading to calculate a sparse matrix multiplied by a dense matrix compared to the Matlab's current implementation. I know you just started this project. Thank you for initializing the effort.

Answer 1

[James Tursa]

MATLAB stores sparse matrices in the CSC 0-based format, except MATLAB does not store the pointerB and pointerE info the same way. MATLAB combines these into one array called Jc. So you would have to convert this part. But the values and row indexes are stored exactly the same in column ordered fashion.

If you are simply passing pointers to the MKL library, and the integers are the same size, you should be able to use the following:

mxGetPr(etc) for values
mxGetIr(etc) for rows
mxGetJc(etc) for pointerB
mxGetJc(etc)+1 for pointerE

If you are using complex data and R2017b memory model then you will have to copy the separate real & imaj data to interleaved format first.

[Anthony Barone]

James,

Thanks so much for the answer and the code...youve undoubetly saved me several hours of trial and error.

Ive marked the answer as "accepted", but if you had the time for a few quick followup questions Id be extremely grateful.

1) THIS is the main operation im interested in implementing, which is C = A * B * A' where B is a symmetric matrix and both A and B are sparse but C is dense. This requires the sorted CSR format. Am I correct in thinking I could use MATLAB's sparse "0-based CSC" format directly and just apply the transpose to the other A matrix? (B is symmetric, so transposing it shouldnt matter).

1a) Side question: does MATLAB store the row-indicies as sorted? MKL has a function to sort them, but I'd rather not take the time to call it if I dont need to.

2) Im using 2018a (which I would assume uses 2017b's memory model) and the arrays (Definitely B, probably A as well) are complex. However, the data are IEEE singles. Could I bit encode all the complex IEEE singles into IEEE doubles where the 1st 32 bits are the real part and the 2nd 32 bits are the imaginary part (using typecast), and then pass that to MKL? (Matlab would show meaningless values, but the bit ordering would be correct). This would actually be fantastic, since I had been struggling for a good way to implement this without casting to/from doubles (since MATLAB doesnt support sparse single arrays).

Thanks again James. you just saved me a ton of work trying to figure this out by trial and error.

EDIT: THIS routine might actually be the one i end up using.

[James Tursa]

1) You state that B is sparse, but THIS link shows that B is dense. Which is it? To be clear, you really do want the Hermetian conjugate transpose operation A', correct? I don't think you can just use the CSC format using the "transposed" operation of the operation you really want because the indexes will not be correct, but I will think about this issue some more. You may have to transpose on the MATLAB side first to obtain the equivalent of the CSR format.

1a) Yes, MATLAB stores the indexes sorted.

2) In 2018a, complex data is always stored as interleaved real/imag data ... even when you compile using the -R2017b memory model (the default). The difference is that the -R2017b compile flag forces a deep copy-in/copy-out of the complex data, definitely something you do NOT want to do. So you need to compile with the -R2018a compile flag for sure.

Yes, you could encode the real/imag singles into an interleaved array which could then be interpreted as a double array and passed through to your library which would reinterpret it as complex single data. How to go about doing this depends on how you have the data (and the indexes) stored in MATLAB to begin with. Maybe you could elaborate the details of this. It would be ideal if the single data were interleaved to begin with.

The THIS link mentions "... a sparse matrix in the internal data structure ..." as the input. I will have to research what this is exactly (or maybe you could tell me if you know already). A struct with pointer fields that need to be filled out? Or ...? I am guessing you will have to build this from scratch using the MATLAB mxArray data pointers, but that is just a guess at this point until I know more about what it really is.

[Anthony Barone]

The 1st link was what I was originally going to use, since the description is "Computes the symmetric product of three sparse matrices and stores the result as a dense matrix", though the actual details suggest that B must be dense, so Ill probably use the 2nd link (where output C is sparse).

Also, regarding "... a sparse matrix in the internal data structure ...": Im not 100% sure, but since this is part of the "sparse inspector-executor" routines that typically is run as a 2 step process (analyse then execute), Id guess it is reffering to the output of mkl_sparse_optimize. Also, I dont have a link, but im pretty sure i saw mkl's documentation refer to "the internal data structure" in as the array handle output by a mkl function.

Thanks for the info on the sorting and the tip about using the -2018a flag.

As far as how the data is stored in matlab and whether it is interleaved...well its complicated, and Im really not sure. Let me quickly (well...somewhat quickly) describe the algorithm as it currently is.

The raw data is seismic data. Basically, there are a number of source points and each source point has a number of (active) receiver stations. The data is a 2D array, where one dimension represents the unique combination of source and receiver points, and the other represents a 1D time series (of 32-bit floats) collected by the receiver. Data is grouped such that

{(t1,x1),(t2,x1),...,(tN,x1),(t1,x2),(t2,x2),...,(tN,x2),(t1,x3), ... etc.}

First a 1D fft is done along the dimension corresponding to time. This introduces the imaginary part. (Negative frequencies are dropped and re-added before the ifft, since the time-domain is real valued).

Next, an analysis is implemented to determine the indices of a sparse matrix such that each source position gets its own row (or column) and each receiver position gets its own column (or row). This is followed by a call of the form

sparse(columnInds, rowInds, data(curFreq,:), numCols, numRows)

Reciprocity is used in this step, meaning that anywhere there is data at {source @ (x1,y1), receiver @ (x2,y2)} but not at {source @ (x2,y2), receiver @ (x1,y1)}, the data we have is copied over into the missing data location. If there is data at both locations, these are averaged. This is equivilant to saying

data = data + data.'; data = data ./ (array composed of 1's and 2's, where points that have both a source and a receiver have a 2).

This makes the data array symmetric. It also makes transposing it trivia, since I have the actual row and column indices of each point (note: the x and y indicies are linearized, so that the 2 dimensions of the sparse array represent sources abd receivers...otherwise the array would be 4D). Currently this step also involves casting to doubles, though Id rather avoid that.

Note: the data is the B matrix from the previous comments.

The main loop loops over frequency. It does a few things to the data, but the important one for this is that it regularizes both sources and receivers to a to a regular grid. The A matrix is computed on the fly (this also computes row and column indicies and uses them in a sparse call).

The A matrix is basically "locally windows data, weighs the data within the window, and sums it". Basically convolution with a non-stationary kernal. In general, the operation is dataOut = A1 * dataIn * A2, where A and B both apply windows to sources and receivers (these are independent operations...this is basically splitting a 4D convolution into 2x 2D convolutions that are implemented using linear indicies. Getting the indexing right was a major pain in the ass...). However, because both dataIn and dataOut have reciprocity applied (and thus are symmetric), I can do dataOut = A * dataIn * transpose(A).

Ive checked my notes and, as it is currently setup, this is definitely the transpose, NOT the hermition conjugate. The data array is symmetric in the same way (i.e., dataIn = dataIn.' ; NOT dataIn=dataIn') (The same holds true for dataOut). I realize this might be a problem, and am working on a way to make things work with hermition conjugates instead, but if it was possible to use standard transposes with complex inputs that would be ideal.

[James Tursa]

"... definitely the transpose, NOT the hermition conjugate ..."

Then you can't use those MKL BLAS routines for this directly, since they specifically prohibit using the straight transpose form with complex variables.

[Anthony Barone]

"Then you can't use those MKL BLAS routines for this directly"

That does seem to be what the documentation says, but the documentation contradicts itself in a few places (e.g., if B is sparse or dense in mkl_sparse_?_syprd), so I figure its possible (though unlikely) complezx transposes are silently supported and the documentation just hasnt been updated. Like I said Im alsop working to convert these to using actual hermition conjugates (My guy tells me this is possible, but I havent quite figured it out yet).

Again, thank you for all the help James. Youve made this much easier than it otherwise would have been. If I could ask 1 last thing....

I wrote a prototype for calling mkl_sparse_spmm. I havent actually tried using it yet (Im tempted, but its like 4:30 am my time and I should probably get a few hours of sleep...). As I said previously though my C is iffy at best...if you have a minute to quickly look it over and let me know if Im making any major errors anywhere. Once I get this working for 1 function call im confident I can build more, but (assuming I didnt somehow get it exactly right on the first try) getting that first function right is much easier if I know whether or not im at least close.

Thanks again James.

    /*==============================================================
     * sparseMatrixMultiply.c - Example for illustrating how to use 
     *      BLAS within a C MEX-file. sparseMatrixMultiply uses the 
     *      MKL Sparse BLAS function 'mkl_sparse_sp2m'.
     *      
     * Auxillary function calls are also made to:
     *     'mkl_sparse_c_create_csc' and 'mkl_sparse_convert_csr'
     *
     * C = sparseMatrixMultiply(A,B) computes the matrix product of A*B,
     *     where A, B, and C are sparse matrices containing real 64-bit
     *     doubles representing bit-encoded 32-bit complex singles.
     *
     * This is based off of the MATLAB mex routine 'matrixMultiply.c'
     *============================================================*/

    #include "mex.h"
    #include "matrix.h"
    #include "mkl_spblas.h"

    void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
    {
        /* DEFINE VARIABLES */

        double *A, *B, *C; /* pointers to input & output matrices*/
        size_t nRow0, nCol0, nRow1, nCol1, nRow2, rowInd0, rowInd1, rowInd2, pntrB0, pntrE0, pntrB1, pntrE1, pntrB2, pntrE2;      /* matrix dimensions and sparse indicies */

        // form of op(A) & op(B) to use in matrix multiplication
        char *opA = "";
        char *opB = "";

        // description of symmetries in A and B. For now assume general sparse matricies.
        char *descrA = "SPARSE_MATRIX_TYPE_GENERAL";
        char *descrB = "SPARSE_MATRIX_TYPE_GENERAL";

        /* EXTRACT INFO FROM INPUTS */

        // get input matrix values
        A = mxGetPr(prhs[0]);    /* first input matrix */
        B = mxGetPr(prhs[1]);    /* second input matrix */

        // dimensions of input matrices
        nRow0 = mxGetM(prhs[0]);
        nCol0 = mxGetN(prhs[0]);
        nRow1 = mxGetM(prhs[1]);
        nCol1 = mxGetN(prhs[1]);

        // get row indicies 
        rowInd0 = mxGetIr(prhs[0]);
        rowInd1 = mxGetIr(prhs[1]);

        // get pointers B and E
        pntrB0 = mxGetJc(prhs[0]);
        pntrE0 = mxGetJc(prhs[0])+1;
        pntrB1 = mxGetJc(prhs[1]);
        pntrE1 = mxGetJc(prhs[1])+1;

        // check that sizes are compatable (i.e., nCols0 == nRows1)
        if ((nCol0 != nRow1) && (nCol0 != nCol1) && (nRow0 != nRow1) && (nRow0 != nCol1) ) {
            mexErrMsgIdAndTxt("MATLAB:matrixMultiply:matchdims", "Inner dimensions of matrix multiply do not match.");
        } else if (nCol0 == nRow1) {
            opA = "SPARSE_OPERATION_NON_TRANSPOSE";
            opB = "SPARSE_OPERATION_NON_TRANSPOSE";
        } else if (nCol0 == nCol1) {
            opA = "SPARSE_OPERATION_NON_TRANSPOSE";
            opB = "SPARSE_OPERATION_TRANSPOSE";
        } else if (nRow0 == nRow1) {
            opA = "SPARSE_OPERATION_TRANSPOSE";
            opB = "SPARSE_OPERATION_NON_TRANSPOSE";
        } else {
            opA = "SPARSE_OPERATION_TRANSPOSE";
            opB = "SPARSE_OPERATION_TRANSPOSE";        
        }

        /* CALL MKL's SPARSE BLAS LIBRARY TO IMPLEMENT 2-STEP INSTECTOR-EXECUTOR SPARSE MATRIX MULTIPLY SUBROUTINE */

        // (PART 1) Get handles for use with MKL...This converts things to MKL's internal format
        // This call should (I think) interpret the 64-bit real bit-encoded doubles as 32-bit complex singles   

        // A matrix (1st input) -->  cscA
        status = mkl_sparse_c_create_csc(cscA, SPARSE_INDEX_BASE_ZERO, &nRow0, &nCol0, &pntrB0, &pntrE0, &rowInd0, &A);
        // B matrix (2nd input) -->  cscB
        status = mkl_sparse_c_create_csc(cscB, SPARSE_INDEX_BASE_ZERO, &nRow1, &nCol1, &pntrB1, &pntrE1, &rowInd1, &B);

        // (PART 2) Pass arguments to MKL and implement 2-step inspector-executor mkl_sparse BLAS operation

        // STAGE 1A - Analyse and compute row indicies of output    
        status = mkl_sparse_sp2m(op0, descrA, cscA, op1, descrB, cscB, SPARSE_STAGE_NNZ_MULT, &cscC);

        // STAGE 1B - Estimate memory requirements for output
        status = mkl_sparse_c_create_csc(cscC, SPARSE_INDEX_BASE_ZERO, &nRow2, &nCol2, &pntrB2, &pntrE2, &rowInd2, &C);
        MKL_INT nnz2 = pntrB2[rowInd2] - pntrE2[0];

        // STAGE 2 - Finish Computation (compute row start/end and values)
        status = mkl_sparse_sp2m(opA, descrA, cscA, opB, descrB, cscB, SPARSE_STAGE_FINALIZE_MULT, &cscC);

        /* EXPORT RESULTS */

        // Convert internal representation back to CSC format (by converting to CSR and transposing)
        status = mkl_sparse_convert_csr(cscC, SPARSE_OPERATION_TRANSPOSE, plhs[0]);     

    }

[James Tursa]

Let's start with the prototype of mkl_sparse_c_creat_csc from the doc:

sparse_status_t mkl_sparse_c_create_csc (sparse_matrix_t *A, sparse_index_base_t indexing, MKL_INT rows, MKL_INT cols, MKL_INT *cols_start, MKL_INT *cols_end, MKL_INT *row_indx, MKL_Complex8 *values);

The general comment here is that you are not currently using the exact variable types specified in the prototype. You might get lucky and they will match, or you might not be lucky in which case the call will probably crash MATLAB. Also some of your arguments are just plain wrong. E.g.,

1st input argument: sparse_matrix_t *A

Currently I don't even see where in your code you even defines these variables, so I wouldn't expect your code to even compile. I would have expected to see something like this in your code:

sparse_matrix_t cscA, cscB, cscC;

And then your calls would be something like this:

status = mkl_sparse_c_create_csc(&cscA,...
    :
status = mkl_sparse_c_create_csc(&cscB,...
    :
status = mkl_sparse_c_create_csc(&cscC,...

This is based on the assumption that sparse_matrix_t is some type of struct in the background, and that the mkl_sparse_c_create_csc is simply attaching pointers and setting integer values in this struct.

2nd input argument: sparse_index_base_t indexing

This looks correct in your code.

3rd & 4th input arguments: MKL_INT rows, MKL_INT cols

In your code you have stuff like this: &nRow0, &nCol0

I.e., where the prototype calls for passing in integers, you are passing in pointers to integers. This is not correct. You should get rid of the & and pass in the integers directly, e.g., just nRow0, nCol0.

5th & 6th & 7th input arguments: MKL_INT *cols_start, MKL_INT *cols_end, MKL_INT *row_indx

Here is where you might be lucky or unlucky. The sparse integer indexing type in MATLAB is mwIndex, which is probably going to be a size_t on your system. As long as this is the same size integer as MKL_INT then you will probably be OK. But if it is different from MKL_INT then your code will likely crash MATLAB. You could check this by comparing sizeof(mwIndex) with sizeof(MKL_INT) to see if you are in trouble. If they don't match then you will have to deep copy these arrays (painful but necessary).

In addition, you are not passing in the proper variables. Here is what you have:

size_t ... rowInd0, ... pntrB0, pntrE0 ...;
  :
rowInd0 = mxGetIr(prhs[0]);
  :
pntrB0 = mxGetJc(prhs[0]);
pntrE0 = mxGetJc(prhs[0])+1;
  :
status = mkl_sparse_c_create_csc(... &pntrB0, &pntrE0, &rowInd0 ...);

You are supposed to be passing in the pointers to the integer data directly, but you don't even have the variables typed correctly. Here is what it should look like:

mwIndex ... *rowInd0, ... *pntrB0, *pntrE0 ...;
  :
rowInd0 = mxGetIr(prhs[0]);
  :
pntrB0 = mxGetJc(prhs[0]);
pntrE0 = mxGetJc(prhs[0])+1;
  :
status = mkl_sparse_c_create_csc(... pntrB0, pntrE0, rowInd0 ...);

But, as stated earlier, this will only work if sizeof(mwIndex)==sizeof(MKL_INT). You might even have a signed/unsigned mismatch that you will have to cast away (assuming the integer data is not too large to spill over into the signed bit). If they don't, then you will need to deep copy these integer arrays.

8th input argument: MKL_Complex8 *values

I don't even see where you have A or B defined anywhere, so this would also not compile. Plus you are not passing the correct type anyway. I would expect something like this:

MKL_Complex8 *A;
  :
A = (MKL_Complex8 *)mxGetData(prhs[0]);
  :
status = mkl_sparse_c_create_csc(... ,A);

That is, you are getting a "packed" double sparse input from the caller, where the packing is two single floating point values occupying the memory of one double value, and you are just reinterpreting the memory as interleaved complex singles.

[James Tursa]

For the mkl_sparse_sp2m calls, it is unclear to me from reading the doc what the arguments need to be, so I will have to get back to you on that. E.g., does the routine dynamically allocate memory for the output argument, or is that supposed to be pre-allocated prior to entering the routine? I am trying to find this (and other) info ...

That being said, I don't think you have the op stuff correct. You have stuff like this:

        char *opA = "";
            :
            opA = "SPARSE_OPERATION_NON_TRANSPOSE";

when it should probably be something like this instead:

        sparse_operation_t opA;
            :
            opA = SPARSE_OPERATION_NON_TRANSPOSE;

Similar comments for descrA etc.

I think what is supposed to happen is this with regards to the output variable:

1) You first allocate an empty sparse matrix of the required size so that it basically only has the Jc index array allocated.

2) You build an MKL sparse matrix with this Jc.

3) You call the multiply routine with the special request flag so that it only calculates NNZ of the result (i.e., fills in the Jc array values).

4) You allocate the index arrays and value arrays of the result and attach them to your sparse matrix (i.e., you rebuild it).

5) You call the multiply routine, where the output is now pre-allocated to handle a result up to NNZ elements.

You are missing some of the steps above.

[Anthony Barone]

Thanks again for all the fantastic info James. I had to temporarily put this aside to deal with something time-sensitive, but getting this working is on my to-do list for today.

As a side note - this is a good example of why yoou shouldnt write codes in languages you dont know very well when you havent slept for something like 30 hours...for some reason it didnt "click" at the time that things like 'MKL_INT' and 'sparse_operation_t' were classes that variables could be defined as (so I was mainly trying to define the variables using the nearest analogue I knew, (e.g. double, size_t, char, ...)

I just re-downloaded the newest MKL update (2018_3) and will post some (hopefully working) revised code in a bit.

I might (for the moment) just switch to the SPARSE_STAGE_FINALIZE_MULT operating mode, since that seems more straightforward. This seems like the type of thing that would be easier to add in later when I know everything else works properly.

[Anthony Barone]

James,

Sorry it took a few days longer than expected...had some computer trouble with a new computer.

I almost have it compiling, but there are a couple of warnings/errors im having trouble figuring out. If you have any thoughts / suggestions on fixing these id be very appreciative.

On Windows I am using visual studio 2017. On linux I am using gcc (4.4.x I believe...the system is running an older centOS version).

I am linking to the sparse blas libraries of MKL 2018 v2 by starting matlab using the following commands:

WINDOWS

cd "C:\Program Files (x86)\IntelSWTools\parallel_studio_xe_2018.3.054\compilers_and_libraries_2018\windows\mkl\bin" && @call mklvars.bat intel64 && set BLAS_VERSION=mkl_rt.dll && set LAPACK_VERSION=mkl_rt.dll && cd "%systemdrive%%homepath%\Documents\MATLAB" && "C:\Program Files\MATLAB\R2018a\bin\matlab.exe"

LINUX

cd "/workdisk/matlab/intel/parallel_studio_xe_2018.3.051/compilers_and_libraries_2018/linux/mkl/bin" && ./mklvars.sh intel64 && export "BLAS_VERSION=/workdisk/matlab/intel/parallel_studio_xe_2018.3.051/compilers_and_libraries_2018/linux/mkl/lib/intel64_lin/libmkl_rt.so" && export "LAPACK_VERSION=/workdisk/matlab/intel/parallel_studio_xe_2018.3.051/compilers_and_libraries_2018/linux/mkl/lib/intel64_lin/libmkl_rt.so" && cd ${HOME}/Documents/MATLAB && /opt/MATLAB/$(ls /opt/MATLAB/ | sort -g | tail -n 1)/bin/matlab -nodesktop -nosplash

MAIN FUNCTION

    /*==============================================================
     * sparseMatrixMultiply.c - Example for illustrating how to use 
     *    BLAS within a C MEX-file. sparseMatrixMultiply uses the 
     *    MKL Sparse BLAS function 'mkl_sparse_sp2m'.
     *    
     * Auxillary function calls are also made to:
     *   'mkl_sparse_c_create_csc' and 'mkl_sparse_convert_csr'
     *
     * C = sparseMatrixMultiply(A,B) computes the matrix product of A*B,
     *   where A, B, and C are sparse matrices containing real 64-bit
     *   doubles representing bit-encoded 32-bit complex singles.
     *
     * This is based off of the MATLAB mex routine 'matrixMultiply.c'
     *============================================================*/

    #include "mex.h"
    #include "matrix.h"
    #include "mkl_spblas.h"
    //#include "/workdisk/matlab/intel/parallel_studio_xe_2018.3.051/compilers_and_libraries_2018/linux/mkl/include/mkl_spblas.h"

    void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
    {
        /* DEFINE VARIABLES */

        // pointers to matricies that are input from matlab
        MKL_Complex8 *A, *B; 
        mxArray *C;

        // pointers to matricies that are used in MKL
        sparse_matrix_t cscA, cscB, cscC;

        // matrix dimensions and sparse indicies 
        mwIndex nRow0, nCol0, nRow1, nCol1, rowInd0, rowInd1, pntrB0, pntrE0, pntrB1, pntrE1;
        MKL_INT nRow0m, nCol0m, nRow1m, nCol1m, rowInd0m, rowInd1m, pntrB0m, pntrE0m, pntrB1m, pntrE1m, nnz2;   

        // form of op(A) & op(B) to use in matrix multiplication
        sparse_operation_t op_T = SPARSE_OPERATION_TRANSPOSE;
        sparse_operation_t op_NT = SPARSE_OPERATION_NON_TRANSPOSE;  
        sparse_operation_t op_CT = SPARSE_OPERATION_CONJUGATE_TRANSPOSE;
        sparse_operation_t opA, opB;

        // description of symmetries in A and B. For now assume general sparse matricies.
        //sparse_matrix_type_t descrA, descrB;
        //descrA = SPARSE_MATRIX_TYPE_GENERAL;
        //descrB = SPARSE_MATRIX_TYPE_GENERAL;
        struct matrix_descr descrA = SPARSE_MATRIX_TYPE_GENERAL;
        struct matrix_descr descrB = SPARSE_MATRIX_TYPE_GENERAL;

        // description of sparse indexing type 
        sparse_index_base_t sparseIndType = SPARSE_INDEX_BASE_ZERO;

        // description of alrogithm stage namespace
        sparse_request_t sparseStage = SPARSE_STAGE_FULL_MULT;

        // status for mkl_sparse operations
        sparse_status_t mkl_status;

        /* EXTRACT INFO FROM INPUTS */

        // get input matrix values
        A = (MKL_Complex8 *)mxGetData(prhs[0]);  /* first input matrix */
        B = (MKL_Complex8 *)mxGetData(prhs[1]);  /* second input matrix */

        // dimensions of input matrices
        nRow0 = mxGetM(prhs[0]);
        nCol0 = mxGetN(prhs[0]);
        nRow1 = mxGetM(prhs[1]);
        nCol1 = mxGetN(prhs[1]);

        // get row indicies 
        rowInd0 = *mxGetIr(prhs[0]);
        rowInd1 = *mxGetIr(prhs[1]);

        // get pointers B and E
        pntrB0 = *mxGetJc(prhs[0]);
        pntrE0 = *mxGetJc(prhs[0])+1;
        pntrB1 = *mxGetJc(prhs[1]);
        pntrE1 = *mxGetJc(prhs[1])+1;

        // check to make sure that MKL_INT and mxIndex have the same number of bits
        if ( sizeof(MKL_INT) != sizeof(mwIndex) ) {
            mexErrMsgIdAndTxt("", "mxIndex and MKL_INT do not have the same size. A fix for this has not yet been implemented (yet).");
        } else {
            nRow0m = nRow0;
            nRow1m = nRow1;
            nCol0m = nCol0;
            nCol1m = nCol1;
            pntrB0m = pntrB0;
            pntrB1m = pntrB1;
            pntrE0m = pntrE0;
            pntrE1m = pntrE1;
            rowInd0m = rowInd0;
            rowInd1m = rowInd1;
        }

        // check that array sizes are compatable (i.e., nCols0 == nRows1)
        // should the size be ambiguous, prioritize operations with fewer transposes
        if ((nCol0 != nRow1) && (nCol0 != nCol1) && (nRow0 != nRow1) && (nRow0 != nCol1) ) {
            mexErrMsgIdAndTxt("MATLAB:matrixMultiply:matchdims", "Inner dimensions of matrix multiply do not match.");
        } else if (nCol0 == nRow1) {
            opA = op_NT;
            opB = op_NT;
        } else if (nCol0 == nCol1) {
            opA = op_NT;
            opB = op_T;
        } else if (nRow0 == nRow1) {
            opA = op_T;
            opB = op_NT;
        } else {
            opA = op_T;
            opB = op_T;    
        }

        /* CALL MKL's SPARSE BLAS LIBRARY TO IMPLEMENT 2-STEP INSTECTOR-EXECUTOR SPARSE MATRIX MULTIPLY SUBROUTINE */

        // (PART 1) Get handles for use with MKL...This converts things to MKL's internal format
        // This call should (I think) interpret the 64-bit real bit-encoded doubles as 32-bit complex singles   

        // A matrix (1st input) -->  cscA
        mkl_status = mkl_sparse_c_create_csc(&cscA, sparseIndType, nRow0m, nCol0m, &pntrB0m, &pntrE0m, &rowInd0m, A);

        // B matrix (2nd input) -->  cscB

        mkl_status = mkl_sparse_c_create_csc(&cscB, sparseIndType, nRow1m, nCol1m, &pntrB1m, &pntrE1m, &rowInd1m, B);

        // (PART 2) Pass arguments to MKL and implement 2-step inspector-executor mkl_sparse BLAS operation

        // USE COMBINED METHOD "SPARSE_STAGE_FULL_MULT"  -->  cscC
        mkl_status = mkl_sparse_sp2m(opA, descrA, cscA, opB, descrB, cscB, sparseStage, &cscC);  

        /* EXPORT RESULTS */

        // Convert internal representation back to CSC format (by converting to CSR and transposing)  -->  plhs[0]
        mkl_status = mkl_sparse_convert_csr(cscC, op_T, &cscC);   

        // output to be sent to matlab as a mxArray
        plhs = mxCreateSparse(nRow0, nCol1, 0, mxREAL);

        // Convert CSC matrix representation back to sparse mxArray  
        C = &cscC;
        mxSetData(plhs[0], C);  

    }

I get the following warnings:

Lines 85-94: these are warning from where the row / column iundicies get convereted to MKL_INT's

warning C4267: '=': conversion from 'size_t' to 'int', possible loss of data.

Lines 139 and 142: these are warnings from where I convert back to an mxArray (the 2 commands preceding the final `mxSetData` call)

(139): warning C4047: '=': 'mxArray **' differs in levels of indirection from 'mxArray *'
(142): warning C4047: '=': 'mxArray *' differs in levels of indirection from 'sparse_matrix_t *'

ERRORS from lines 46-47: These are the `struct matrix_descr ...`

 (46): error C2440: 'initializing': cannot convert from 'int' to 'matrix_descr'
   (47): error C2440: 'initializing': cannot convert from 'int' to 'matrix_descr'

Ive tried several variations (including some using the sparse_matrix_type_t declaration, since this is also mentioned in the documentation), but I cant seem to figure out how to correctly declare these.

Thanks again James. You've been a HUGE help in this.

[Anthony Barone]

So, I actually think ive figured out the main error from lines 46-47. I was initializing the structure all wrong. This seems to work.

 typedef struct {
        sparse_matrix_type_t type;
        sparse_fill_mode_t mode;
        sparse_diag_type_t diag;
    } matrix_descr;

    struct matrix_descr descrA;
    struct matrix_descr descrB;

    descrA.type = SPARSE_MATRIX_TYPE_GENERAL;
    descrB.type = SPARSE_MATRIX_TYPE_GENERAL;

I'm still getting a warning where I re-assign the output to a mxArray, but I imagine I can fix this...I just need to take a closer look at the functions for doing this and the examples that matlab provides.

First though, i need to update gcc (once I fixed the struct definition error i got a notice that mex officially supports gcc 6.3.x, and 4.4.7-17 is unsupported). Unfortunately I dont have root access on this machine, and I think to get a local copy of gcc I need to build it from source (which isnt THAT hard, but is somewhat time consuming).

[James Tursa]

For converting the output back to mxArray, I thought our strategy was to build the output mxArray first, then build the output MKL sparse matrix from that. That way the MKL sparse output would already be part of the output mxArray and there would be no conversion to do. (All subject to mwIndex being the same size as MKL_INT of course).

But in any event, I will wait until your next post giving me your status and we can work from there.