A sparse matrix vector multiplication in C mex programming

Hi everyone,

I'm new in C mex programming and I'm really stuck. I would like to write a code that computes a sparse matrix vector multiplication (sparse_matrix*vector=result_vector). Given how sparse matrices are managed in C mex, I have difficulties getting the non zero elements of each row (to do the multiplication with the vector), because the pointers Ir=mxGetIr(sparse_matrix) and Jc=mxGetJc(sparse_matrix) give the information about the columns only and not the rows.

I would be really grateful if somebody helps me out because I can't move on without solving this problem.

Thanks in advance.

Sincerely,

J

[Matt J]

Why must you re-implement an operation that MATLAB has already implemented for you?

In any case, why must you implement sparse_matrix(i,:)*vector? Why can't you implement sum_i sparse_matrix(:,i)*vector(i)?

[James Tursa]

I'm with Matt J. Why are you doing this? And if for some reason you really do need to do this, what have you done so far? Where specifically are you stuck? The mxGetIr( ) gives the row indexing, btw. Is the vector full or sparse?

[JM]

Thanks to both of you for you answers. Matt J I'm not coding in matlab but I'm coding in C and then I compile from matlab using mex engine, and as far as I know, there's no routine already implemented that does sparse matrix vector multiplication in C.

As you know, in order to carry out the multiplication, each row of the sparse matrix is multiplied by the vector which is a full one, and then the sum of the multiplication result gives the result that is stored in the result vector. So how can I, for each row of the sparse matrix, access the non zeros elements in order to multiply them with the vector ? I can't see how to do this with Ir and Jc because they contain the information about the columns only and not the rows.

Please let me know if I'm not clear enough.

Thanks in advance.

J

Answer 1

[James Tursa]

Here is the bare bones C code for this. Most of the code consists of argument checking etc. The actual computation part of the code is only a few lines and is fairly trivial. It simply accumulates the individual element multiplies into the appropriate resulting element. The j for-loop loops over the columns of M, and the while loop does the individual element multiplies within each column.

/* Bare bones sparse matrix * full vector                               */
/* Produces a full vector result                                        */
/* Assumes all input elements are finite                                */
/* E.g., does not account for NaN or 0*inf etc.                         */
/* Programmer:  James Tursa                                             */

/* Includes ----------------------------------------------------------- */

#include "mex.h"

/* Gateway ------------------------------------------------------------ */

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    mwSize m, n, j, nrow;
    double *Mpr, *Vpr, *Cpr;
    mwIndex *Mir, *Mjc;

/* Argument checks */
    if( nlhs > 1 ) {
        mexErrMsgTxt("Too many outputs");
    }
    if( nrhs != 2 ) {
        mexErrMsgTxt("Need exactly two inputs");
    }
    if( !mxIsDouble(prhs[0]) || !mxIsSparse(prhs[0]) || mxIsComplex(prhs[0]) ) {
        mexErrMsgTxt("1st argument must be real double sparse matrix");
    }
    if( !mxIsDouble(prhs[1]) || mxIsSparse(prhs[1]) || mxIsComplex(prhs[1]) ||
        mxGetNumberOfDimensions(prhs[1])!=2 || mxGetN(prhs[1])!=1 ) {
        mexErrMsgTxt("2nd argument must be real double full column vector");
    }
    m = mxGetM(prhs[0]);
    n = mxGetN(prhs[0]);
    if( mxGetM(prhs[1])!=n ) {
        mexErrMsgTxt("Inner matrix dimensions must agree.");
    }

/* Create output */
    plhs[0] = mxCreateDoubleMatrix( m, 1, mxREAL );

/* Get data pointers */
    Mpr = mxGetPr(prhs[0]);
    Vpr = mxGetPr(prhs[1]);
    Cpr = mxGetPr(plhs[0]);
    Mir = mxGetIr(prhs[0]);
    Mjc = mxGetJc(prhs[0]);

/* Calculate result */
    for( j=0; j<n; j++ ) {
        nrow = Mjc[j+1] - Mjc[j];  /* Number of row elements for this column */
        while( nrow-- ) {
            Cpr[*Mir++] += *Mpr++ * Vpr[j];  /* Accumulate contribution of Vpr[j] */
        }
    }
}

[JM]

Thank you so much James, now I understand many things related to the handling of sparse matrices under C mex.

Cheers.

J

[Royi Avital]

Isn't there an internal function to do this efficiently (Multi Threaded)?