[ublas] possible bug (probably in diagonal_matrix)

From: Vardan Akopian (vakopian_at_[hidden])
Date: 2005-06-23 03:46:50

Hello ublas users,

First post, and I'm afraid it's a bug report. Here is the simplest
code I could come up with to illustrate:

#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/banded.hpp>
#include <boost/numeric/ublas/io.hpp>

int main()
{
        using namespace boost::numeric::ublas;

        size_t sz = 3;

        diagonal_matrix<double> z(sz), x(sz);
        matrix<double> y(sz, sz);
        for (size_t i = 0; i < sz; ++i) {
                z(i,i) = i;
                x(i,i) = 2*i;
                for (size_t j = 0; j < sz; ++j) {
                        y(i,j) = i + j;
                }
        }

        // this produces a correct result
        matrix<double> x_y = x - y;
        matrix<double> A = prod(z, x_y);

        // this produces a wrong result
        matrix<double> B = prod(z, x - y);

        std::cout << "A=" << A << std::endl;
        std::cout << "B=" << B << std::endl;
        return 0;
}

A is computed correctly, while B is not (they should be equal). Am I
doing something wrong/illegal?  I've tested this with boost version
1.32.0 and the CVS head; gcc 3.4.3 and gcc 4 under Linux.
Interestingly, when NDEBUG is off, ublas's internal runtime check
catches the problem with an assertion failure:
Assertion failed in file
../../../../boost_1_32_0/boost/numeric/ublas/matrix_assign.hpp at line
797:
equals (m, cm)
terminate called after throwing an instance of
'boost::numeric::ublas::external_logic'
  what():  external logic
Aborted

I'd be very interested if someone would confirm and/or propose a fix.
Thanks.

-Vardan Akopian
Jet Propulsion Laboratory, NASA