Path: news.mathworks.com!newsfeed-00.mathworks.com!fu-berlin.de!uni-berlin.de!news.dfncis.de!not-for-mail
From: Joachim Selke <j.selke@tu-bs.de>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Sparse matrix product again
Date: Mon, 12 May 2008 17:11:05 +0200
Organization: Technische =?ISO-8859-1?Q?Universit=E4t_Braunschweig?=
Lines: 64
Message-ID: <68r509F2v2e1qU1@mid.dfncis.de>
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> <fvv6ff$c92$1@fred.mathworks.com>
Mime-Version: 1.0
Content-Type: text/plain; charset=ISO-8859-1
Content-Transfer-Encoding: 7bit
X-Trace: news.dfncis.de FRvPtLFJ1FiYH5wJ8Cf14A4J1/iZ6LaYt7f802T9BmQYllZxsCqnASRkEJ
Cancel-Lock: sha1:bjHdGO1AItcW95nw/dPnEjhM4jw=
User-Agent: Mozilla/5.0 (X11; U; Linux x86_64; en-US; rv:1.8.1.14) Gecko/20080501 Fedora/2.0.0.14-1.fc8 Thunderbird/2.0.0.14 Mnenhy/0.7.5.0
In-Reply-To: <fvv6ff$c92$1@fred.mathworks.com>
Xref: news.mathworks.com comp.soft-sys.matlab:467920


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

Thanks. I also changed the loop order of the inner loop as you suggested
in your code. Testing for zero makes it slightly faster. Here is my
updated code:


-------------------------------------------------------------------

// In case you compile it with gcc:
// mex -largeArrayDims CFLAGS='$CFLAGS -std=c99 -funroll-loops'
COPTIMFLAGS='$COPTIMFLAGS -O3' sparsemult_mex.c

#include "mex.h"

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

    const mxArray *a = prhs[0];
    const mxArray *b = prhs[1];
    const mxArray *x = prhs[2];

    const mwSize m = mxGetM(x);
    const mwSize k = mxGetN(a);
    const mwSize n = mxGetN(x);

    mxArray *lhs[1];
    mxArray *rhs[1] = {mxDuplicateArray(a)};
    mexCallMATLAB(1, lhs, 1, rhs, "transpose");
    const mxArray *atr = lhs[0];

    mxArray *ab = mxDuplicateArray(x);

    mwSignedIndex *rows_ab = mxGetIr(ab);
    mwSignedIndex *cols_ab = mxGetJc(ab);

    double *data_atr = mxGetPr(atr);
    double *data_b = mxGetPr(b);
    double *data_ab = mxGetPr(ab);

    for (mwSignedIndex j = 0; j < n; j++, data_b += k) {
        /* data_b points to the first entry of B's j-th column*/
        for (mwSignedIndex rows_left = cols_ab[j + 1] - cols_ab[j];
rows_left > 0; rows_left--, data_ab++, rows_ab++) {
            /* rows_ab points to AB's current entry (i, j) */
            /* data_ab points to AB's current entry (i, j) */
            mwSignedIndex i = *rows_ab;
            double z = 0;
            for (mwSignedIndex r = k - 1; r >= 0; r--) {
                z += data_atr[i * k + r] * data_b[r];
            }
            *data_ab = z;
        }
    }

    plhs[0] = ab;
}

-------------------------------------------------------------------

Joachim