[ublas] Sparse/Packed matrix assignment needs type conversion

Subject: [ublas] Sparse/Packed matrix assignment needs type conversion
From: Marco Guazzone (marco.guazzone_at_[hidden])
Date: 2010-07-08 11:05:34

Hi,

I've found a possible problem in matrix assignment when the target
matrix has a packed storage.

For instance, try this:

--- [code] ---
#include <complex>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/symmetric.hpp>

namespace ublas = boost::numeric::ublas;

int main()
{
    typedef double in_value_type;
    typedef std::complex<double> out_value_type;

    const std::size_t n = 4;

    ublas::matrix<in_value_type> IN(n,n);
    ublas::matrix<out_value_type> OUT(IN); // COMPILE
    OUT = IN; // COMPILE

    ublas::symmetric_matrix<in_value_type,ublas::lower> sym_IN(n,n);
    ublas::symmetric_matrix<out_value_type,ublas::lower>
sym_OUT(sym_IN); // NOT COMPILE
    OUT = IN; // NOT COMPILE
}
--- [/code] ---

Looking in the uBLAS code, I've realized that the problem is in
matrix_assign and also interests sparse matrices.

The problem.

Let:
  typedef typename M::value_type value_type;
  typedef F<typename M::..., typename E::value_type> functor_type;

Then:
1.
  if (v != value_type/*zero*/()) ... // where v is of type E::value_type;
  [ e.g. E::value_type is std::complex<float> and M::value_type is
std::complex<double> ]

2. functor_type::apply(*it, value_type/*zero*/()); // we are assigning
an M::value_type to an E::value_type
  [ e.g. E::value_type is float and M::value_type is complex<float> ]

The (possible) solution.

Expressions of the first type might be changed by casting an
E::value_type to a M::value_type, like in this way:
  if (static_cast<value_type>(v) != value_type/*zero*/()) ... // where
v is of type E::value_type
Obviously, this does not work when E::value_type and M::value_type are
not castable (e.g., std::complex and double, respectively)

Expressions of the second type might be changed in this 2 ways:
Option A
  typedef typename matrix_traits<E>::value_type expr_value_type;
  functor_type::apply(*it, expr_value_type/*zero*/()); // NOTE: use of
E::value_type in place of M::value_type
or alternatively:
Option B
  typedef F<typename M::..., value_type> functor_type; // NOTE: use
M::value_type instead of E::value_type
  functor_type::apply(*it, value_type/*zero*/()); // unchanged

I have opened a new ticket
  https://svn.boost.org/trac/boost/ticket/4410

Cheers,

-- Marco