Re: [Boost-Users] uBLAS Matrix Multiplication

From: Kresimir Fresl (fresl_at_[hidden])
Date: 2002-11-07 07:03:00

(I added `Cc: ublas-dev' because I think that this can be interesting
there, too.)

Toon Knapen wrote:

> On Wednesday 06 November 2002 16:41, Hickman, Greg wrote:

>>Why isn't there an overloaded

>>matrix_expression
>>operator* (const matrix_expression&, const matrix_expression&);

>>for peforming matrix multiplication?

> This issue has been discussed before and the trouble IIRC is that operator* 
> could also mean element-wise matrix multiplication. IMHO the optimal solution 
> would be to take over all conventions from Matlab, a convention most are 
> familiar with, but the element-wise multiplication operator '.*' can't be 
> used in C++ ;( Another suggestion was to defined the element-wise and the 
> normal multiplication in a different namespace (I still like this idea very 
> much)

Let's try:
======================================================
#include <iostream>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>

namespace boost { namespace numeric { namespace ublas {

     namespace normal {
           template<class E1, class E2>
           BOOST_UBLAS_INLINE
           typename matrix_matrix_binary_traits<
                 typename E1::value_type, E1, typename E2::value_type, E2
           >::result_type
           operator * (const matrix_expression<E1> &e1,
                            const matrix_expression<E2> &e2) {
                  return prod (e1, e2);
            }
      }
}}}

int main() {
      boost::numeric::ublas::matrix<double> m1 (2, 2), m2 (2, 2);
      // init m1 & m2
      using namespace boost::numeric::ublas::normal;
      std::cout << m1 * m2 << std::endl;
}
======================================================

g++ 3.2 complains (`slightly' edited ;o) :
======================================================
ambiguous overload for `matrix<double>& * matrix<double>&' operator

candidates are:

(defined here):
matrix_matrix_binary_traits<>::result_type
operator*(const matrix_expression<E>&, const matrix_expression<E2>&)
[with E1 = matrix<double>, E2 = matrix<double>]

(defined in `matrix_expression.hpp'):
matrix_binary_scalar2_traits<>::result_type
operator*(const matrix_expression<E>&, const T2&)
[with E1 = matrix<double>, T2 = matrix<double>]

(defined in `matrix_expression.hpp'):
matrix_binary_scalar1_traits<>::result_type
operator*(const T1&, matrix_expression<E2>&)
[with T1 = matrix<double>, E2 = matrix<double>]
======================================================
Comeau's compiler 4.3.0.1 gives similar diagnostics.

In `operator*(matrix_expression, T)'  and `operator*(T, matrix_expression)'
T is meant to be scalar, but compilers do not know this.

Let's try to tell them (probably can be simplified):
====================================================
template<class E2>
BOOST_UBLAS_INLINE
typename matrix_binary_scalar1_traits<
     typename E2::value_type, E2,
     scalar_multiplies<typename E2::value_type, typename E2::value_type>
 >::result_type
operator * (const typename E2::value_type &e1,
                   const matrix_expression<E2> &e2) {   //  etc.
====================================================

It works; not only `m1 * m2', but also `m1 * trans (m2)',
`m1 * (m1 + m2)' etc.

Interesting thing is that if first parameter is declared as
`const typename matrix_expression<E2>::value_type &e1'
neither como nor g++ can find appropriate operator*
for e.g. `2. * m1'.

So, header `operator_mult.hpp' can be:
=================================================
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
namespace boost { namespace numeric { namespace ublas {

     namespace elemwise {
         // there are no elementwise multiplications in ublas,
         // at least not yet
     }

     namespace normal {

         template<class E1, class E2>
         BOOST_UBLAS_INLINE
         typename vector_scalar_binary_traits<
              E1,
              E2,
              vector_inner_prod<
                    typename E1::value_type,
                    typename E2::value_type,
                    typename promote_traits<
                          typename E1::value_type,
                          typename E2::value_type
                    >::promote_type
              >
         >::result_type
         operator * (const vector_expression<E1> &e1,
                            const vector_expression<E2> &e2)
         {
                return inner_prod (e1, e2);
         }

         template<class E1, class E2>
         BOOST_UBLAS_INLINE
         typename matrix_vector_binary1_traits<
                typename E1::value_type,
                E1,
                typename E2::value_type,
                E2
         >::result_type
         operator * (const matrix_expression<E1> &e1,
                           const vector_expression<E2> &e2)
         {
               return prod (e1, e2);
         }

         template<class E1, class E2>
         BOOST_UBLAS_INLINE
         typename matrix_vector_binary2_traits<
               typename E1::value_type,
               E1,
               typename E2::value_type,
               E2
         >::result_type
         operator * (const vector_expression<E1> &e1,
                            const matrix_expression<E2> &e2)
         {
               return prod (e1, e2);
         }

         template<class E1, class E2>
         BOOST_UBLAS_INLINE
         typename matrix_matrix_binary_traits<
                typename E1::value_type,
                E1,
                typename E2::value_type,
                E2
         >::result_type
         operator * (const matrix_expression<E1> &e1,
                            const matrix_expression<E2> &e2)
         {
               return prod (e1, e2);
         }

     }
}}}
================================================
but, before it can be used, interfaces of existings operators * in
matrix_expression.hpp and vector_expression.hpp must be
adapted as described.

There's one more, independant reason to change the interface
of the operator* (scalar, Matrix_or_Vector_expression) and
related/similar operators to something like

    operator * (const typename E2::value_type &e1,
                       const M_or_V_expression<E2> &e2)

Currently, one sometimes gets rather surprising and confusing
results:
==============================================
   boost::numeric::ublas::matrix<double> m1 (2, 2);
   m1 (0, 0) = 0.1; m1 (0, 1) = 0.2;
   m1 (1, 0) = 0.3; m1 (1, 1) = 0.4;
   std::cout << 2 * m1 << std::endl;  //  note: `2', not `2.',
   std::cout << m1 * 2 << std::endl;  //  i.e. `int', not `double'
==============================================
Output is:

   [2,2]((0,0),(0,0))
   [2,2]((0.2,0.4),(0.6,0.8))

Sincerely,

fres