Fast sparse matrix vector product

sparse matrix vector product

Description

This is a dedicated mex-file to efficiently compute the Sparse-Matrix-Vector-Product (smvp).

x = smvp(A,b);

This is generally 4 times faster than x = A*b for large A (On my machine).

Type "help smvp" and "smvp_demo" for more information.

NOTE: You will need either windows or the MATLAB compiler to use this function. A pre-compiled version is included, but you can re-compile via:

mex -O smvp.c

You may also find my efficient LU substitution function useful:

http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=10546&objectType=File

P.S This again can be useful in the numerical solution of PDE's.

NOTE: The zip file should now work (03/04/2006) but I have removed the demo

• MATLAB release: MATLAB 6.5 (R13)

[Michal Kvasnicka]

The zip file is corrupted. The matrix.mat file can not be extracted.

[Robert Beardmore]

I had no problems using this file and it is around 5 times quicker than matlab's A*b for the problems I tried.

[Shelley Lee]

It would be excellent if it can support
complex arithmetics.

[Martin S.]

Works but offers almost very minor 8% speedup as of Matlab 2006.
I multiplied 700,000 x 700,000 matrix with 2,000,000 nonzero entries times a 700,000 vector and got only 8% faster than Matlab 2006. Not worth using unless there is more speedup in other settings.

2) The routine fails to check that the matrix allociation succeeded (may run out of memory). Need to add a NULL check to

plhs[0] = mxCreateDoubleMatrix(mxGetM(prhs[0]), 1, mxREAL);

3) Change if (nlhs!=1) to if (nlhs > 1)

Summary: works but very limited speedup.

[Tim Davis]

Regarding the comment on checking for NULL return of mxCreateDoubleMatrix: it's not necessary. If the allocation fails, MATLAB will automatically terminate the mexFunction and clean up any allocated memory. If it was malloc, then yes, you would need to check for NULL. This is a fine example of how to use mxCreateDoubleMatrix.

The only design issue that should be corrected is the "nlhs!=1" test (a minor point). If you simply do

>> smvp(A,b)

with no output arguments, then nlhs is set to zero. You can still write to plhs[0], and the result will go in the "ans" variable.

[Stephen Becker]

There's actually a very important point here. Because of how Matlab stores a sparse matrix (see Tim Davis' book, or read the help on, say, mxGetIc), and because a major bottleneck for these computations is loading data into the CPU's cache, it is MUCH faster to store a matrix by row, instead of column, when doing matrix-vector multiples.

In Matlab, this is easy to remedy: simply take the transpose.

So, this cool trick can save you time. If you want to do this repeatedly:
>> y = A*b
instead, do this:
>> At = A.'; % a one-time cost
>> y = At.'b

But, this trick only works on newer versions of Matlab. On, say, R2006b, a call like
>> y = At.'b
will actually calculate the transpose of At, so this is very slow!

On older versions of Matlab, it's likely that this SMVP code will really help, but it might not be as helpful on R2008, for example.

On linux, non-multithreaded R2006a, I found that SMVP took about 60% the time of Matlab's own multiply. So, good work!