Re: [ublas] More observations on ublas sparse matrix multiplication performance

From: Michael Stevens (mail_at_[hidden])
Date: 2005-09-22 14:26:35

On Donnerstag 22 September 2005 01:53, Frank Astier wrote:
> Ian -
>
> Thanks for the correction. It got rid of the bizarre effect that
> ublas seemed to always to the same amount of work, regardless of
> sparsity :-)

It is not all doom and gloom! In "operation.hpp" there axpy_prod which help 
considerably. Using in the test loop
       axpy_prod (m, v, r1);

I get.
0.1 0.32 0.42=
0.2 0.35 0.42=
0.3 0.46 0.43=
0.4 0.56 0.42=
0.5 0.56 0.42=
0.6 0.79 0.43=
0.7 0.95 0.42=
0.8 1.07 0.42=
0.9 1.17 0.42=

and looking in more detail
0 0.02 0.49=
0.02 0.12 0.42=
0.04 0.12 0.42=
0.06 0.13 0.41=
0.08 0.12 0.42=
0.1 0.24 0.42=
0.12 0.24 0.42=
0.14 0.23 0.42=
0.16 0.23 0.42=
0.18 0.24 0.42=
0.2 0.24 0.42=
0.22 0.35 0.42=
0.24 0.34 0.42=
0.26 0.35 0.43=
0.28 0.35 0.42=

So even at 20% sparsity we are doing better. I think this is pritty good if 
not glittering.

> I implemented my own sparse matrix too, optimized for multiplication
> with dense vectors (I didn't even write random access to elements, I
> don't need it!). Below 10% non zeros, I get at least 3 times better
> than the naive dense multiplication, and up to 40 times better than
> ublas on small matrix sizes (around 256x256). The matrices I handle
> don't have any structure I can take advantage of, unfortunately.
>
> I wish ublas documentation ad enough details about what to expect.

I think one could write several PhD on matrix multiplications involving sparse 
matrices.

That said the normal ublas::prod is unlikely to be a good choice where sparse 
matrics are involved. The are many ways the basic element access loops can be 
ordered to achive a multiplication. axpy_prod use one ordering that iterates 
over the matrix.
For prod to work in general expressions it uses an ordering which iterates of 
the elements in the result. For sparse matrices this is bad news as we have 
to do a great deal of work to find where we are in the sparse matrix.

That said I would love to see this documented further. I have some material of 
my own with diagrams showing the element access paterns I could contribute.

Michael

___________________________________
Michael Stevens Systems Engineering

34128 Kassel, Germany
Phone/Fax: [phone]

Navigation Systems, Estimation  and
                 Bayesian Filtering
    http://bayesclasses.sf.net
___________________________________