More details about a possible bug in MTL

Hello again.

As posted a few days ago, we have an issue with simple matrix  
multiplication over a dense symmetric matrix with MTL. We have been  
trying to understend it, but debugging MTL is very complex unless one  
knows the package in its deepest details.

The problem can be shown by the following simple example:

#include <iostream>
#include <mtl/matrix.h>
#include <mtl/mtl.h>
#include <mtl/utils.h>
#include <mtl/linalg_vec.h>

using namespace mtl;
using namespace std;

typedef matrix< double, symmetric<upper>, dense<>, row_major>::type  
Matrix;

typedef dense1D<double> Vec;

int main()
{
     const Matrix::size_type matrix_size=5;

     Matrix A(matrix_size,matrix_size);
     Vec y(matrix_size,0), x(matrix_size, 1), z(matrix_size);

     int r;
     for (r = 0;r < matrix_size;++r) {
         (A)(r,(r)) = 1;
         if (r + 1 < matrix_size)
             (A)(r,(r+1)) = 1;
     }
     mult(A,x,y,z);

     return 0;
}

The bug occurs inside

template <class Matrix, class VecX, class VecZ>
inline void
mult_symm__(const Matrix& A, const VecX& x, VecZ& z, row_tag)
{
   typedef typename matrix_traits<Matrix>::value_type T;
   typename Matrix::const_iterator i;
   typename Matrix::OneD::const_iterator j, jend;

   for (i = A.begin(); i != A.end(); ++i) {
     T tmp = z[i.index()];
     j = (*i).begin();
     jend = (*i).end();
     if (A.is_upper()) {
       tmp += *j * x[j.column()];
       ++j;
     } else
       --jend;
     for (; j != jend; ++j) {
       /* normal side */
       tmp += *j * x[j.column()];
       /* symmetric side */
       z[j.column()] += *j * x[j.row()];
     }
     if (A.is_lower())
       tmp += *j * x[j.column()];
     z[i.index()] = tmp;
   }
}

and can by easily spotted by having jend.column() printed at every  
iteration of the for() loop. It is a 5-by-5 dense matrix, so one  
would assume jend.column() to be always equal to 5. However, this is  
not the case: jend.column() is actually 5 at the first iteration, but  
then it becomes 6, 7, 8, and 9.  This results in wrong data to be  
read from x[], and, even worse, data to be written outside the bounds  
of z[].

As I said, we tried to follow the sequence of who is jend and why the  
value is wrong, but we after the 5th consecutive template  
instantiation we sadly lost our way.

By changing to

typedef matrix< double, symmetric<lower>, dense<>, row_major>::type  
Matrix;
                                   ^^^^^

the problem "sort of disappears": the iterators do not go out of  
bounds. However, each computation appears to be performed twice, and  
one of the two copies gets overwritten.

A symmetric pair of occurrences happen with

typedef matrix< double, symmetric<***er>, dense<>,  
column_major>::type Matrix;
                                   ^^^^^            ^^^^^^^^^^^^

When the matrix is <upper> the jend.row() goes out of bounds, when  
the matrix is <lower> the result is correct but the work is done twice.

It actually looks a pretty funny bug. Any help in this matter will be  
greatly appreciated.

			Best Regards

			  Antonio


_______________________________________________
This list is archived at http://www.osl.iu.edu/MailArchives/mtl-devel/