Re: [ublas] Bad prod efficiency of sparse/symmetry matrix ?

From: Gunter Winkler (guwi17_at_[hidden])
Date: 2006-12-08 05:58:47

On Thursday 07 December 2006 16:48, Seabook wrote:
> (code)

1) there is a bug in axpy_prod:

other() is a static function and should be used as such.

Index: operation.hpp
===================================================================
RCS file: /cvsroot/boost/boost/boost/numeric/ublas/operation.hpp,v
retrieving revision 1.23
diff -u -p -r1.23 operation.hpp
--- operation.hpp	16 Jul 2006 13:13:01 -0000	1.23
+++ operation.hpp	8 Dec 2006 09:02:41 -0000
@@ -642,7 +642,7 @@ namespace boost { namespace numeric { na
                 typename matrix_column<expression1_type>::const_iterator itc (mc.begin ());
                 typename matrix_column<expression1_type>::const_iterator itc_end (mc.end ());
                 while (itc != itc_end) {
-                    if (triangular_restriction::functor_type ().other (itc.index (), it2.index2 ()))
+                    if (triangular_restriction::other (itc.index (), it2.index2 ()))
                         m (itc.index (), it2.index2 ()) += *it1 * *itc;
                     ++ itc;
                 }

2) my times (g++4.1 -DNDEBUG -O2) are
prod()
sp*sp, 160, 1620, 8520, 65900
sp*dc, 180, 1660, 9130, 70300
sp*dr, 170, 1740, 9140, 70500
sp*sy, 200, 1640, 8900, 70200
sp*sy, 190, 1710, 8830, 69600
dc*dr, 60, 780, 9430, 135100
axpy_prod()
0, 0, 10, 10
70, 240, 940, 3900
70, 260, 930, 3900
50, 200, 860, 3100
50, 210, 830, 3400
70, 400, 3520, 26600

the last row is the ublas dense*dense time. the first column is size=100x100, 
the second column is size=200x200, third is size=400x400, fourth is 
size=800x800. You see that the time of the sparse algorithm grows much 
slower than the time of the dense version. However I have to state that 
I modified you initialize code: The assign sp(i,j)=0 adds a (logical) 
nonzero to the sparse matrix, because there is no check whether you assign 0
or something else. Thus your sparse matrices were completely filled 
with zeros and that's why they computed extremely slow. Note also, 
that the time of axpy_prod(compressed, compressed, compressed) is 0, 
because the result is nearly empty. Another remark is, that 
sparse = prod(sparse, dense) is a bad design, because the result is dense.

Summary: 
* prod() works best for (small to medium) dense and structured matrices
* axpy_prod() is much more optimized for sparse matrices 
* the advantage of sparse matrices becomes larger with increasing matrix size
* The optimal workload should be O(nnz^2) for sparse O(nnz*n^2) for mixed and
  O(n^3) for dense products.

mfg
Gunter