Path: news.mathworks.com!not-for-mail
From: "Andy Robb" <ajrobb@hotmail.com>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Sparse matrix product again
Date: Thu, 8 May 2008 15:33:03 +0000 (UTC)
Organization: The MathWorks, Inc.
Lines: 36
Message-ID: <fvv6ff$c92$1@fred.mathworks.com>
References: <67rg20F2pn19mU1@mid.dfncis.de> <fvfd1n$1dc$1@fred.mathworks.com> <688agnF2rcjvjU1@mid.dfncis.de> <fvmtmc$k80$1@fred.mathworks.com> <68dnitF2soa18U1@mid.dfncis.de>
Reply-To: "Andy Robb" <ajrobb@hotmail.com>
NNTP-Posting-Host: webapp-05-blr.mathworks.com
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1210260783 12578 172.30.248.35 (8 May 2008 15:33:03 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Thu, 8 May 2008 15:33:03 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 1374918
Xref: news.mathworks.com comp.soft-sys.matlab:467422


Your code might not compile with some C compilers as you
intermix declarations and statements. I have replaced
separate declarations and assignments with 'definitions'.

I don't have a MEX compiler system handy so this is untested!

#include "mex.h"

void mexFunction
( int nlhs, mxArray *plhs[]
, int nrhs, const mxArray *prhs[]
) {
  const mxArray *atr = prhs[0], *b = prhs[1], *x = prhs[2];
  const mwSignedIndex m = mxGetM(x), n = mxGetN(x), k =
mxGetM(atr);
  mxArray *ab = mxDuplicateArray(x);
  const mwIndex *rows = mxGetIr(ab), *cols = mxGetJc(ab);
  const double *data_atr = mxGetPr(atr), data_b = mxGetPr(b);
  double * data_ab = mxGetPr(ab);

  while (--n >= 0) {
    double *A = atr + cols[n];
    double *B = b + rows[n] * k;
    double z = 0;
    mwSignedIndex i;

    for (i = 0; i < k; ++i) {
      z += *A * B[i];
      A += m;
    }

    data_ab[n] = z;
  }

  plhs[0] = ab;
}