Dear all,



I believe there is an implementation error in the sparse_prod implementation:


line 73 is                                                               m (it1.index1 (), j) = temporary (j);

but (in my opinion) should be                    m (it1.index1 (), j) += temporary (j);


Similarly for line 132                                        m (i, it2.index2 ()) = temporary (i);

But (in my opinion) should be                    m (i, it2.index2 ()) += temporary (i);


Using the current implementation the following two approaches give a different result,


sparse_prod( matrix_left, matrix_right, matrix_result, false );

matrix_result += prod( matrix_left, matrix_right );


while, in fact they should give  the same result. Adapting the sparse_prod implementation as suggested above results in both lines of code giving the same result.



Kind regards,

Johannes Keustermans