MATLAB: A sparse matrix vector multiplication in C mex programming

cmexsparse matrix

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,

Best Answer

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] */
        }
    }
}

Related Solutions

MATLAB: Building Mex file with Fortan

Your main problem is that you need to create the output array PLHS(1) before you access it. Also your NLHS and NRHS types are wrong. Try this:

#include <fintrf.h>
C
C Gateway routine
C
    SUBROUTINE MEXFUNCTION(NLHS, PLHS, NRHS, PRHS)
    implicit none
C mexFunction arguments:
    MWPOINTER PLHS(*), PRHS(*)
    INTEGER NLHS, NRHS
    mwPointer, external :: mxGetPr, mxGetPi
    mwPointer, external :: mxCreateDoubleMatrix, mxCreateNumericArray
    mwSize, external :: mxGetM, mxGetN
    MWPOINTER YPP 
    mwSize M, N
    integer*4 ComplexFlag
    REAL *8 RYPP(1)
    RYPP(1) = 10
    M = 1
    N = 1
    ComplexFlag = 0
    PLHS(1) = mxCreateDoubleMatrix( M, N, ComplexFlag )
    YPP = MXGETPR(PLHS(1))
    CALL MXCOPYREAL8TOPTR(RYPP, YPP, M*N)
    RETURN
    END

MATLAB: Preserving numerical symmetry in large nxn matrix

Here is a mex routine to do this calculation. It relies on inputting the diagonal matrix as a full vector of the diagonal elements. It does not check for underflow to 0 for the calculations. A robust production version of this code would check for this and clean the sparse result of 0 entries, but I did not include that code here. It also does not check for inf or NaN entries. This could be made faster with parallel code such as OpenMP, but I didn't do that either.

/* File:  spdimd.c                                                                 */
/* Compile:  mex spdimd.c                                                          */
/* Syntax  C = spdimd(M,D)                                                         */
/* Does  C = D^-1 * M * D                                                          */
/* where M = double real sparse NxN matrix                                         */
/*       D = double real N element full vector representing diagonal NxN matrix    */
/*       C = double real sparse NxN matrix                                         */
/* Programmer:  James Tursa                                                        */
/* Date: 2/24/2021                                                                 */
/* Includes ----------------------------------------------------------- */
#include "mex.h"
#include <string.h>
/* Gateway ------------------------------------------------------------ */
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    mwSize m, n, j, nrow;
    double *Mpr, *Dpr, *Cpr;
    mwIndex *Mir, *Mjc, *Cir, *Cjc;
/* 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 || (mxGetM(prhs[1]) != 1 && mxGetN(prhs[1]) != 1)) {
        mexErrMsgTxt("2nd argument must be real double full vector");
    }
    m = mxGetM(prhs[0]);
    n = mxGetN(prhs[0]);
    if( m!=n || mxGetNumberOfElements(prhs[1])!=n ) {
        mexErrMsgTxt("Matrix dimensions must agree.");
    }
/* Sparse info */
    Mir = mxGetIr(prhs[0]);
    Mjc = mxGetJc(prhs[0]);
/* Create output */
    plhs[0] = mxCreateSparse( m, n, Mjc[n], mxREAL);
    
/* Get data pointers */
    Mpr = (double *) mxGetData(prhs[0]);
    Dpr = (double *) mxGetData(prhs[1]);
    Cpr = (double *) mxGetData(plhs[0]);
    Cir = mxGetIr(plhs[0]);
    Cjc = mxGetJc(plhs[0]);
/* Fill in sparse indexing */
    memcpy(Cjc, Mjc, (n+1) * sizeof(mwIndex));
    memcpy(Cir, Mir, Cjc[n] * sizeof(mwIndex));
/* Calculate result */
    for( j=0; j<n; j++ ) {
        nrow = Mjc[j+1] - Mjc[j];  /* Number of row elements for this column */
        while( nrow-- ) {
            *Cpr++ = *Mpr++ * (Dpr[j] / Dpr[*Cir++]);  
	}
    }
}

Best Answer

Related Solutions

MATLAB: Building Mex file with Fortan

MATLAB: Preserving numerical symmetry in large nxn matrix

Related Question