|
|
Ublas : |
From: Gunter Winkler (guwi17_at_[hidden])
Date: 2007-03-23 17:38:37
Hello,
below is a (possibly broken) patch to matrix.hpp which shows how to intercept
the expression template machine and map the expression
A.assign( prod(B,C) );
to a call to gemm. The crucial part is to identify the correct type of the
argument of assign. (However, neither noalias(A)=prod(B,C) nor A=prod(B,C)
work with this specialization).
mfg
Gunter
Index: matrix.hpp
===================================================================
RCS file: /cvsroot/boost/boost/boost/numeric/ublas/matrix.hpp,v
retrieving revision 1.72
diff -u -p -r1.72 matrix.hpp
--- matrix.hpp 6 Dec 2006 09:34:01 -0000 1.72
+++ matrix.hpp 23 Mar 2007 21:31:17 -0000
@@ -20,6 +20,12 @@
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix_expression.hpp>
#include <boost/numeric/ublas/detail/matrix_assign.hpp>
+
+#ifdef BOOST_UBLAS_AUTOBLAS
+#include <boost/numeric/ublas/functional.hpp>
+#endif
// Iterators based on ideas of Jeremy Siek
@@ -208,6 +214,14 @@ namespace boost { namespace numeric { na
matrix_assign<scalar_assign> (*this, ae);
return *this;
}
+#ifdef BOOST_UBLAS_AUTOBLAS
+ BOOST_UBLAS_INLINE
+ matrix &assign (const matrix_matrix_binary<matrix, matrix, matrix_matrix_prod<matrix, matrix, value_type> > &ae) {
+ std::cout << "gemm(1.0,B,C,0.0,*this)\n";
+ matrix_assign<scalar_assign> (*this, ae);
+ return *this;
+ }
+#endif
template<class AE>
BOOST_UBLAS_INLINE
matrix& operator += (const matrix_expression<AE> &ae) {