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