--- In Boost-Users@y..., John Lloyd wrote:

Hickman, Greg wrote:

...

Why isn't there an overloaded

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

for peforming matrix multiplication?

Greg Hickman

Info: http://www.boost.org Wiki: <http://www.crystalclearsoftware.com/cgi- bin/boost_wiki/wiki.pl> Unsubscribe: mailto:boost-users-unsubscribe@y...

Your use of Yahoo! Groups is subject to http://docs.yahoo.com/info/terms/

I think it has something to do with the fact that expression templates can be very inefficient for certain chained operations, e.g., A * B * v, where A and B are matrices and v is a vector.

Suppose you have a matrix expression like A = B * C * D. Using expression templates, this could result in an order n^5 time complexity as opposed to order n^3 complexity.

Exactly. I believe the theoretically more important signature is prod (A, B), although prod (A, B) will be more frequently used practically (BTW: the complexity of prod (A, B) (i, j) is O(n) as opposed to O(n^3) for prod (A, B) (i, j)). Best regards Joerg

----- Original Message ----- From: alexei_novakov To: Boost-Users@yahoogroups.com Sent: Sunday, November 10, 2002 11:29 PM Subject: [Boost-Users] Re: uBLAS Matrix Multiplication [snip]

Here are couple of thoughts: 1) matrix-matrix multiplication and matrix-vector multiplication are ones of most used operators in matrix computations (especially second), if these are not presented

Do you mean in form of operator*()?

the value of expression framework itself becomes questionable.

Which signature would you propose for blocked matrix multiplications (with additional parameters like block size) then?

2) To avoid extra computational compexity temporaries can be used. In some cases temporaries are presented in expressions already, like: A = B * C; Matrix A can be used as temporary.

If you write A.assign (prod (B, C)) in uBLAS, exactly this happens.

More complex expressions like: A = B * C * D; will require additional temp matrix.

And a complex ET analysis (maybe even transformation) of the whole statement.

3) To minimize performance penalty template specialization can be used. For example: matrix_matrix_prod { ... }; matrix_matrix_prod { ... }; etc.

Either etc. is infinite or we'd have to impose assumptions which restrict programmer's control over uBLAS. Best regards Joerg

(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 #include namespace boost { namespace numeric { namespace ublas { namespace normal { template 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

::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 #include namespace boost { namespace numeric { namespace ublas { namespace elemwise { // there are no elementwise multiplications in ublas, // at least not yet } namespace normal { template 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 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 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 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