AMDiS/lib/mtl4/boost/numeric/mtl/operation/mult.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef MTL_MULT_INCLUDE
+#define MTL_MULT_INCLUDE
+
+#include <boost/numeric/mtl/config.hpp>
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/utility/category.hpp>
+#include <boost/numeric/mtl/utility/flatcat.hpp>
+#include <boost/numeric/mtl/utility/ashape.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/mtl/utility/static_assert.hpp>
+#include <boost/numeric/mtl/operation/dmat_dmat_mult.hpp>
+#include <boost/numeric/mtl/operation/smat_smat_mult.hpp>
+#include <boost/numeric/mtl/operation/smat_dmat_mult.hpp>
+#include <boost/numeric/mtl/operation/mat_vec_mult.hpp>
+#include <boost/numeric/mtl/operation/rvec_mat_mult.hpp> // Row vector times matrix
+#include <boost/numeric/mtl/operation/mult_specialize.hpp>
+#include <boost/numeric/mtl/operation/assign_mode.hpp>
+#include <boost/numeric/mtl/operation/mult_assign_mode.hpp>
+#include <boost/numeric/mtl/utility/enable_if.hpp>
+
+#include <boost/mpl/if.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+
+namespace mtl { namespace matrix {
+
+
+/// Multiplication: mult(a, b, c) computes c= a * b;
+/** The 3 types must be compatible, e.g. all three matrices or b and c are column vectors and a is a matrix.
+    The dimensions are checked at compile time. **/
+template <typename A, typename B, typename C>
+typename mtl::traits::enable_if_matrix<A>::type
+inline mult(const A& a, const B& b, C& c)
+{
+    vampir_trace<4010> tracer;
+    MTL_DEBUG_THROW_IF(static_cast<const void*>(&a) == static_cast<const void*>(&c)
+		       || static_cast<const void*>(&b) == static_cast<const void*>(&c),
+		       argument_result_conflict());
+
+    // dispatch between matrices, vectors, and scalars
+    using mtl::traits::shape_flatcat;
+    gen_mult(a, b, c, assign::assign_sum(), shape_flatcat<A>(), shape_flatcat<B>(), shape_flatcat<C>());
+                                            // typename category<A>::type(), typename category<B>::type(), typename category<C>::type());
+}
+
+
+/// Multiplication: mult_add(a, b, c) computes c+= a * b;
+/** The 3 types must be compatible, e.g. all three matrices or b and c are column vectors and a is a matrix.
+    The dimensions are checked at compile time. **/
+template <typename A, typename B, typename C>
+typename mtl::traits::enable_if_matrix<A>::type
+inline mult_add(const A& a, const B& b, C& c)
+{
+    vampir_trace<4010> tracer;
+    // dispatch between matrices, vectors, and scalars
+    using mtl::traits::shape_flatcat;
+    gen_mult(a, b, c, assign::plus_sum(), shape_flatcat<A>(), shape_flatcat<B>(), shape_flatcat<C>());
+                                          // typename category<A>::type(), typename category<B>::type(), typename category<C>::type());
+}
+
+
+/// Four term multiplication: mult(a, x, y, z) computes z= a * x + y;
+/** The 4 types must be compatible, i.e. a*x must be assignable to z and z must be incrementable by y.
+    Right now, it is not more efficient than z= a * x; z+= y. For compatibility with MTL2. **/
+template <typename A, typename X, typename Y, typename Z>
+inline void mult(const A& a, const X& x, const Y& y, Z& z)
+{
+    vampir_trace<4010> tracer;
+    mult(a, x, z);
+    z+= y;
+}
+
+
+// Matrix multiplication
+template <typename MatrixA, typename MatrixB, typename MatrixC, typename Assign>
+inline void gen_mult(const MatrixA& a, const MatrixB& b, MatrixC& c, Assign, tag::flat<tag::matrix>, tag::flat<tag::matrix>, tag::flat<tag::matrix>)
+{
+    vampir_trace<4011> tracer;
+#if 1
+    MTL_DEBUG_THROW_IF((const void*)&a == (const void*)&c || (const void*)&b == (const void*)&c, argument_result_conflict());
+#else
+    if ((const void*)&a == (const void*)&c || (const void*)&b == (const void*)&c) {
+	C tmp(num_rows(c), num_cols(c));
+	mult(a, b, tmp);
+	swap(C, tmp);
+	return;
+    }
+#endif
+
+    MTL_DEBUG_THROW_IF(num_rows(a) != num_rows(c) || num_cols(a) != num_rows(b) || num_cols(b) != num_cols(c), incompatible_size());
+    // dispatch between dense and sparse
+    using mtl::traits::sparsity_flatcat;
+    mat_mat_mult(a, b, c, Assign(), sparsity_flatcat<MatrixA>(), sparsity_flatcat<MatrixB>(), sparsity_flatcat<MatrixC>());
+		 // typename category<MatrixA>::type(), typename category<MatrixB>::type(), typename category<MatrixC>::type());
+}
+
+
+/// Dense matrix multiplication
+/**  The function for dense matrix multiplication defines a default multiplication functor.
+     Alternatively the user can define his own functors for specific triplets of matrix types,
+     see detail::dmat_dmat_mult_specialize.
+     The default functor for dense matrix multiplication is:
+     -# Use BLAS   if available, otherwise
+     -# Recursive multiplication with:
+        -# Platform optimized mult on blocks   if available, otherwise
+        -# Tiled multiplication on blocks      if available, otherwise
+        -# Naive multiplication on blocks
+     -# Naive multiplication on entire matrices if recursion is not available
+**/
+template <typename MatrixA, typename MatrixB, typename MatrixC, typename Assign>
+inline void mat_mat_mult(const MatrixA& A, const MatrixB& b, MatrixC& c, Assign, tag::flat<tag::dense>, tag::flat<tag::dense>, tag::flat<tag::dense>)
+{
+    vampir_trace<4012> tracer;
+    using assign::plus_sum; using assign::assign_sum;
+
+    static const unsigned long tiling1= detail::dmat_dmat_mult_tiling1<MatrixA, MatrixB, MatrixC>::value;
+    static const unsigned long tiling2= detail::dmat_dmat_mult_tiling2<MatrixA, MatrixB, MatrixC>::value;
+    typedef gen_tiling_dmat_dmat_mult_t<tiling1, tiling2, plus_sum>    tiling_mult_t;
+
+    typedef gen_platform_dmat_dmat_mult_t<plus_sum, tiling_mult_t>     platform_mult_t;
+    typedef gen_recursive_dmat_dmat_mult_t<platform_mult_t>            recursive_mult_t;
+    typedef gen_blas_dmat_dmat_mult_t<assign_sum, recursive_mult_t>    blas_mult_t;
+    typedef size_switch_dmat_dmat_mult_t<straight_dmat_dmat_mult_limit, tiling_mult_t, blas_mult_t>   variable_size_t;
+
+    typedef fully_unroll_fixed_size_dmat_dmat_mult_t<Assign>           fully_unroll_t;
+    typedef size_switch_dmat_dmat_mult_t<fully_unroll_dmat_dmat_mult_limit, fully_unroll_t, tiling_mult_t> fixed_size_t;
+
+    static const bool all_static= mtl::traits::is_static<MatrixA>::value && mtl::traits::is_static<MatrixB>::value
+	                          && mtl::traits::is_static<MatrixC>::value;
+    typedef static_switch_dmat_dmat_mult_t<all_static, fixed_size_t, variable_size_t>  default_functor_t;
+
+    /// Use user-defined functor if provided (assign mode can be arbitrary)
+    typedef typename boost::mpl::if_<
+	detail::dmat_dmat_mult_specialize<MatrixA, MatrixB, MatrixC>
+      , typename detail::dmat_dmat_mult_specialize<MatrixA, MatrixB, MatrixC>::type
+      , default_functor_t
+    >::type raw_functor_type;
+
+    /// Finally substitute assign mode (consistently)
+    typename assign::mult_assign_mode<raw_functor_type, Assign>::type functor;
+
+    functor(A, b, c);
+}
+
+template <typename MatrixA, typename MatrixB, typename MatrixC, typename Assign>
+inline void mat_mat_mult(const MatrixA& A, const MatrixB& b, MatrixC& c, Assign, tag::flat<tag::dense>, tag::flat<tag::dense>, tag::flat<tag::sparse>)
+{
+    vampir_trace<4012> tracer;
+    // This is a useless and extremely inefficient operation!!!!
+    // We compute this with a dense matrix and copy the result back
+    dense2D<typename Collection<MatrixC>::value_type, matrix::parameters<> > c_copy(num_rows(c), num_cols(c));
+    c_copy= c;
+    mat_mat_mult(A, b, c_copy, Assign(), tag::flat<tag::dense>(), tag::flat<tag::dense>(), tag::flat<tag::dense>());
+    c= c_copy;
+}
+
+/// Sparse matrix multiplication
+template <typename MatrixA, typename MatrixB, typename MatrixC, typename Assign>
+inline void mat_mat_mult(const MatrixA& A, const MatrixB& b, MatrixC& c, Assign, tag::flat<tag::sparse>, tag::flat<tag::sparse>, tag::flat<tag::sparse>)
+{
+    vampir_trace<4012> tracer;
+    smat_smat_mult(A, b, c, Assign(), typename OrientedCollection<MatrixA>::orientation(),
+		   typename OrientedCollection<MatrixB>::orientation());
+}
+
+template <typename MatrixA, typename MatrixB, typename MatrixC, typename Assign>
+inline void mat_mat_mult(const MatrixA& A, const MatrixB& b, MatrixC& c, Assign, tag::flat<tag::sparse>, tag::flat<tag::sparse>, tag::flat<tag::dense>)
+{
+    vampir_trace<4012> tracer;
+    // This is a useless and extremely inefficient operation!!!!
+    // We compute this with a sparse matrix and copy the result back
+    compressed2D<typename Collection<MatrixC>::value_type, matrix::parameters<> > c_copy(num_rows(c), num_cols(c));
+    c_copy= c;
+    smat_smat_mult(A, b, c_copy, Assign(), typename OrientedCollection<MatrixA>::orientation(),
+		   typename OrientedCollection<MatrixB>::orientation());
+    c= c_copy;
+}
+
+/// Product of sparse times dense matrix
+/**  This function (specialization of mult) is intended to multiply sparse matrices with multiple matrices
+     gathered into a dense matrix.  Likewise, the resulting dense matrix corresponds to multiple vectors.
+     The default functor for this operation is:
+     -# Use tiled multiplication      if available, otherwise
+     -# Naive multiplication
+**/
+template <typename MatrixA, typename MatrixB, typename MatrixC, typename Assign>
+inline void mat_mat_mult(const MatrixA& A, const MatrixB& b, MatrixC& c, Assign, tag::flat<tag::sparse>, tag::flat<tag::dense>, tag::flat<tag::dense>)
+{
+    vampir_trace<4012> tracer;
+    using assign::plus_sum; using assign::assign_sum;
+    using namespace functor;
+
+    // static const unsigned long tiling1= detail::dmat_dmat_mult_tiling1<MatrixA, MatrixB, MatrixC>::value;
+
+    //typedef gen_smat_dmat_mult<Assign>                         default_functor_t;
+    typedef gen_tiling_smat_dmat_mult<8, Assign>                         default_functor_t;
+
+    // Finally substitute assign mode (consistently)
+    // typename assign::mult_assign_mode<raw_functor_type, Assign>::type functor;
+
+    default_functor_t functor;
+    functor(A, b, c);
+}
+
+template <typename MatrixA, typename MatrixB, typename MatrixC, typename Assign>
+inline void mat_mat_mult(const MatrixA& A, const MatrixB& b, MatrixC& c, Assign, tag::flat<tag::sparse>, tag::flat<tag::dense>, tag::flat<tag::sparse>)
+{
+    vampir_trace<4012> tracer;
+    // This is a useless and extremely inefficient operation!!!!
+    // We compute this with a sparse matrix and copy the result back
+    dense2D<typename Collection<MatrixC>::value_type, matrix::parameters<> > c_copy(num_rows(c), num_cols(c));
+    c_copy= c;
+    mat_mat_mult(A, b, c_copy, Assign(), tag::flat<tag::sparse>(), tag::flat<tag::dense>(), tag::flat<tag::dense>());
+    c= c_copy;
+}
+
+
+template <typename MatrixA, typename MatrixB, typename MatrixC, typename Assign>
+inline void mat_mat_mult(const MatrixA& A, const MatrixB& b, MatrixC& c, Assign, tag::flat<tag::dense>, tag::flat<tag::sparse>, tag::flat<tag::dense>)
+{
+    vampir_trace<4012> tracer;
+    // This is could be a usefull operation, i.e. multiplying multiple row vectors with a sparse matrix
+    // Might be supported in future
+    // Now we compute this with a sparse matrix as first argument
+    compressed2D<typename Collection<MatrixA>::value_type, matrix::parameters<> > A_copy(num_rows(A), num_cols(A));
+    A_copy= A;
+    compressed2D<typename Collection<MatrixC>::value_type, matrix::parameters<> > c_copy(num_rows(c), num_cols(c));
+    c_copy= c;
+    mat_mat_mult(A_copy, b, c_copy, Assign(), tag::flat<tag::sparse>(), tag::flat<tag::sparse>(), tag::flat<tag::sparse>());
+    c= c_copy;
+}
+
+
+
+template <typename MatrixA, typename MatrixB, typename MatrixC, typename Assign>
+inline void mat_mat_mult(const MatrixA& A, const MatrixB& b, MatrixC& c, Assign, tag::flat<tag::dense>, tag::flat<tag::sparse>, tag::flat<tag::sparse>)
+{
+    vampir_trace<4012> tracer;
+    // This is not a usefull operation, because the result is dense
+    // Now we compute this with a sparse matrix as first argument
+    compressed2D<typename Collection<MatrixA>::value_type, matrix::parameters<> > A_copy(num_rows(A), num_cols(A));
+    A_copy= A;
+    mat_mat_mult(A_copy, b, c, Assign(), tag::flat<tag::sparse>(), tag::flat<tag::sparse>(), tag::flat<tag::sparse>());
+}
+
+
+
+// Matrix vector multiplication
+template <typename Matrix, typename VectorIn, typename VectorOut, typename Assign>
+inline void gen_mult(const Matrix& A, const VectorIn& v, VectorOut& w, Assign, tag::flat<tag::matrix>, tag::flat<tag::col_vector>, tag::flat<tag::col_vector>)
+{
+    vampir_trace<4011> tracer;
+    // Vector must be column vector
+    // If vector is row vector then matrix must have one column and the operation is a outer product
+    //   -> result should be a matrix too
+
+    // Check if element types are compatible (in contrast to tag dispatching, nesting is considered here)
+    MTL_STATIC_ASSERT((boost::is_same< typename ashape::mult_op<typename ashape::ashape<Matrix>::type,
+			                                          typename ashape::ashape<VectorIn>::type >::type,
+			                 ::mtl::ashape::mat_cvec_mult
+				     >::value),
+		      "The type nesting of the arguments does not allow for a consistent matrix vector product.");
+
+
+#if 1
+    MTL_DEBUG_THROW_IF((void*)&v == (void*)&w, argument_result_conflict());
+#else
+    if ((void*)&v == (void*)&w) {
+	VectorOut tmp(size(w));
+	mult(A, b, tmp);
+	swap(w, tmp);
+	return;
+    }
+#endif
+    // w.checked_change_dim(num_rows(A)); // destroys distribution in parallel -> dimension changed in assignment
+    MTL_DEBUG_THROW_IF(num_rows(A) != mtl::vector::size(w), incompatible_size());
+    MTL_DEBUG_THROW_IF(num_cols(A) != mtl::vector::size(v), incompatible_size());
+
+    mat_cvec_mult(A, v, w, Assign(), mtl::traits::mat_cvec_flatcat<Matrix>());
+}
+
+
+// Vector matrix multiplication
+template <typename VectorIn, typename Matrix, typename VectorOut, typename Assign>
+inline void gen_mult(const VectorIn& v, const Matrix& A, VectorOut& w, Assign, tag::flat<tag::row_vector>, tag::flat<tag::matrix>, tag::flat<tag::row_vector>)
+{
+    vampir_trace<4011> tracer;
+    // Vector must be column vector
+    // If vector is row vector then matrix must have one column and the operation is a outer product
+    //   -> result should be a matrix too
+
+    // Check if element types are compatible (in contrast to tag dispatching, nesting is considered here)
+    MTL_STATIC_ASSERT((boost::is_same< typename ashape::mult_op<typename ashape::ashape<VectorIn>::type,
+			                                          typename ashape::ashape<Matrix>::type >::type,
+			                 ::mtl::ashape::rvec_mat_mult
+			             >::value),
+		      "The type nesting of the arguments does not allow for a consistent matrix vector product.");
+
+#if 1
+    MTL_DEBUG_THROW_IF((void*)&v == (void*)&w, argument_result_conflict());
+#else
+    if ((void*)&v == (void*)&w) {
+	VectorOut tmp(size(w));
+	mult(A, b, tmp);
+	swap(w, tmp);
+	return;
+    }
+#endif
+    // w.checked_change_dim(num_cols(A));
+    w.checked_change_resource(v);
+    MTL_DEBUG_THROW_IF(num_cols(v) != num_rows(A), incompatible_size());
+
+    // same dispatching criterion as mat_cvec_mult (until now)
+    rvec_mat_mult(v, A, w, Assign(), typename mtl::traits::mat_cvec_flatcat<Matrix>::type());
+}
+
+
+
+
+
+}} // namespace mtl::matrix
+
+#endif // MTL_MULT_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/mult_assign_mode.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef MTL_MULT_ASSIGN_MODE_INCLUDE
+#define MTL_MULT_ASSIGN_MODE_INCLUDE
+
+#include <boost/numeric/mtl/operation/assign_mode.hpp>
+#include <boost/numeric/mtl/operation/dmat_dmat_mult.hpp>
+#include <boost/numeric/mtl/operation/no_op.hpp>
+
+namespace mtl { namespace assign {
+
+namespace detail {
+
+    template <typename Assign>
+    struct subm_assign
+    {
+	typedef Assign type;
+    };
+
+    template<>
+    struct subm_assign<assign_sum>
+    {
+	typedef plus_sum type;
+    };
+
+}
+
+// Set assign_mode of functor 'Mult' to 'Assign'
+// including assign_mode of backup functors and functors for sub-matrices
+template <typename Mult, typename Assign>
+struct mult_assign_mode
+{};
+
+
+#if 0
+// Omit the fully typed functors; they shouldn't be used directly
+template <typename MatrixA, typename MatrixB, typename MatrixC, typename OldAssign, typename Backup,
+	  typename Assign>
+struct mult_assign_mode<gen_dmat_dmat_mult_ft<MatrixA, MatrixB, MatrixC, OldAssign, Backup>,
+			Assign>
+{
+    typedef gen_dmat_dmat_mult_ft<MatrixA, MatrixB, MatrixC, Assign, Backup> type;
+};
+#endif
+
+
+template <typename Assign>
+struct mult_assign_mode<no_op, Assign>
+{
+    typedef no_op type;
+};
+
+template <typename OldAssign, typename Backup, typename Assign>
+struct mult_assign_mode<gen_dmat_dmat_mult_t<OldAssign, Backup>, Assign>
+{
+    typedef typename mult_assign_mode<Backup, Assign>::type                      backup_type;
+    typedef gen_dmat_dmat_mult_t<Assign, backup_type> type;
+};
+
+template <typename OldAssign, typename Backup, typename Assign>
+struct mult_assign_mode<gen_cursor_dmat_dmat_mult_t<OldAssign, Backup>, Assign>
+{
+    typedef typename mult_assign_mode<Backup, Assign>::type                      backup_type;
+    typedef gen_cursor_dmat_dmat_mult_t<Assign, backup_type> type;
+};
+
+template <unsigned long Tiling1, unsigned long Tiling2, typename OldAssign, typename Backup, typename Assign>
+struct mult_assign_mode<gen_tiling_dmat_dmat_mult_t<Tiling1, Tiling2, OldAssign, Backup>, Assign>
+{
+    typedef typename mult_assign_mode<Backup, Assign>::type                      backup_type;
+    typedef gen_tiling_dmat_dmat_mult_t<Tiling1, Tiling2, Assign, backup_type>   type;
+};
+
+template <typename OldAssign, typename Backup, typename Assign>
+struct mult_assign_mode<gen_tiling_44_dmat_dmat_mult_t<OldAssign, Backup>, Assign>
+{
+    typedef typename mult_assign_mode<Backup, Assign>::type                      backup_type;
+    typedef gen_tiling_44_dmat_dmat_mult_t<Assign, backup_type>                  type;
+};
+
+template <typename OldAssign, typename Backup, typename Assign>
+struct mult_assign_mode<gen_tiling_22_dmat_dmat_mult_t<OldAssign, Backup>, Assign>
+{
+    typedef typename mult_assign_mode<Backup, Assign>::type                      backup_type;
+    typedef gen_tiling_22_dmat_dmat_mult_t<Assign, backup_type>                  type;
+};
+
+template <typename BaseMult, typename BaseTest, typename OldAssign, typename Backup, typename Assign>
+struct mult_assign_mode<gen_recursive_dmat_dmat_mult_t<BaseMult, BaseTest, OldAssign, Backup>, Assign>
+{
+    typedef typename mult_assign_mode<Backup, Assign>::type                      backup_type;
+
+    // Corresponding assignment type for sub-matrices
+    typedef typename detail::subm_assign<Assign>::type                           base_assign_type;
+    typedef typename mult_assign_mode<BaseMult, base_assign_type>::type          base_mult_type;
+
+    typedef gen_recursive_dmat_dmat_mult_t<base_mult_type, BaseTest, Assign, backup_type>  type;
+};
+
+template <typename OldAssign, typename Backup, typename Assign>
+struct mult_assign_mode<gen_platform_dmat_dmat_mult_t<OldAssign, Backup>, Assign>
+{
+    typedef typename mult_assign_mode<Backup, Assign>::type                      backup_type;
+    typedef gen_platform_dmat_dmat_mult_t<Assign, backup_type>                   type;
+};
+
+template <typename OldAssign, typename Backup, typename Assign>
+struct mult_assign_mode<gen_blas_dmat_dmat_mult_t<OldAssign, Backup>, Assign>
+{
+    typedef typename mult_assign_mode<Backup, Assign>::type                      backup_type;
+    typedef gen_blas_dmat_dmat_mult_t<Assign, backup_type>                       type;
+};
+
+template <std::size_t SizeLimit, typename FunctorSmall, typename FunctorLarge, typename Assign>
+struct mult_assign_mode<size_switch_dmat_dmat_mult_t<SizeLimit, FunctorSmall, FunctorLarge>, Assign>
+{
+    typedef typename mult_assign_mode<FunctorSmall, Assign>::type                small_type;
+    typedef typename mult_assign_mode<FunctorLarge, Assign>::type                large_type;
+    typedef size_switch_dmat_dmat_mult_t<SizeLimit, small_type, large_type>      type;
+};
+
+template <bool IsStatic, typename FunctorStatic, typename FunctorDynamic, typename Assign>
+struct mult_assign_mode<static_switch_dmat_dmat_mult_t<IsStatic, FunctorStatic, FunctorDynamic>, Assign>
+{
+    typedef typename mult_assign_mode<FunctorStatic, Assign>::type               static_type;
+    typedef typename mult_assign_mode<FunctorDynamic, Assign>::type              dynamic_type;
+    typedef static_switch_dmat_dmat_mult_t<IsStatic, static_type, dynamic_type>  type;
+};
+
+template <typename OldAssign, typename Backup, typename Assign>
+struct mult_assign_mode<fully_unroll_fixed_size_dmat_dmat_mult_t<OldAssign, Backup>, Assign>
+{
+    typedef fully_unroll_fixed_size_dmat_dmat_mult_t<Assign, Backup>             type;
+};
+
+}} // namespace mtl::assign
+
+#endif // MTL_MULT_ASSIGN_MODE_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/mult_result.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef MTL_MULT_RESULT_INCLUDE
+#define MTL_MULT_RESULT_INCLUDE
+
+#include <boost/numeric/mtl/utility/ashape.hpp>
+#include <boost/numeric/mtl/mtl_fwd.hpp>
+
+#if 0
+#include <boost/numeric/mtl/matrix/map_view.hpp>
+#include <boost/numeric/mtl/matrix/mat_mat_times_expr.hpp>
+#include <boost/numeric/mtl/matrix/implicit_dense.hpp>
+#include <boost/numeric/mtl/vector/map_view.hpp>
+#include <boost/numeric/mtl/operation/mat_cvec_times_expr.hpp>
+#include <boost/numeric/mtl/vector/rvec_mat_times_expr.hpp>
+#endif
+
+namespace mtl { namespace traits {
+
+template <typename Op1, typename Op2, typename MultOp> struct mult_result_aux;
+template <typename Op1, typename Op2, typename MultOp> struct vec_mult_result_aux;
+    //template <typename Op1, typename Op2, typename MultOp1, typename MultOp2> struct mult_result_if_equal_aux;
+
+/// Result type for multiplying arguments of types Op1 and Op2
+/** Can be used in enable-if-style as type is only defined when appropriate.
+    This one is used if at least one argument is a matrix.
+**/
+template <typename Op1, typename Op2>
+struct mult_result
+  : public mult_result_aux<Op1, Op2, typename ashape::mult_op<typename ashape::ashape<Op1>::type,
+							      typename ashape::ashape<Op2>::type >::type>
+{};
+
+
+/// Result type for multiplying arguments of types Op1 and Op2
+/** Can be used in enable-if-style as type is only defined when appropriate.
+    This one is used if at least one argument is a vector and none is a matrix.
+**/
+template <typename Op1, typename Op2>
+struct vec_mult_result
+  : public vec_mult_result_aux<Op1, Op2, typename ashape::mult_op<typename ashape::ashape<Op1>::type,
+								  typename ashape::ashape<Op2>::type >::type>
+{};
+
+
+/// Result type for multiplying arguments of types Op1 and Op2
+/** MultOp according to the algebraic shapes **/
+template <typename Op1, typename Op2, typename MultOp>
+struct mult_result_aux {};
+
+/// Scale matrix from left
+template <typename Op1, typename Op2>
+struct mult_result_aux<Op1, Op2, ::mtl::ashape::scal_mat_mult>
+{
+    typedef matrix::scaled_view<Op1, Op2> type;
+};
+
+
+/// Scale matrix from right needs functor for scaling from right
+template <typename Op1, typename Op2>
+struct mult_result_aux<Op1, Op2, ::mtl::ashape::mat_scal_mult>
+{
+    typedef matrix::rscaled_view<Op1, Op2> type;
+};
+
+/// Multiply matrices
+template <typename Op1, typename Op2>
+struct mult_result_aux<Op1, Op2, ::mtl::ashape::mat_mat_mult>
+{
+    typedef matrix::mat_mat_times_expr<Op1, Op2> type;
+};
+
+/// Multiply matrix with column vector
+template <typename Op1, typename Op2>
+struct mult_result_aux<Op1, Op2, ::mtl::ashape::mat_cvec_mult>
+{
+    typedef mat_cvec_times_expr<Op1, Op2> type;
+};
+
+/// Multiply column with row vector and return implicit matrix
+template <typename Op1, typename Op2>
+struct vec_mult_result_aux<Op1, Op2, ::mtl::ashape::cvec_rvec_mult>
+{
+    typedef mtl::matrix::outer_product_matrix<Op1, Op2> type;
+};
+
+/// Result type for multiplying arguments of types Op1 and Op2
+/** MultOp according to the algebraic shapes **/
+template <typename Op1, typename Op2, typename MultOp>
+struct vec_mult_result_aux {};
+
+/// Scale row vector from left
+template <typename Op1, typename Op2>
+struct vec_mult_result_aux<Op1, Op2, ::mtl::ashape::scal_rvec_mult>
+{
+    typedef vector::scaled_view<Op1, Op2> type;
+};
+
+
+
+
+/// Scale column vector from left
+template <typename Op1, typename Op2>
+struct vec_mult_result_aux<Op1, Op2, ::mtl::ashape::scal_cvec_mult>
+{
+    typedef vector::scaled_view<Op1, Op2> type;
+};
+
+/// Multiply row vector with matrix
+// added by Cornelius and Peter
+template <typename Op1, typename Op2>
+struct vec_mult_result_aux<Op1, Op2, ::mtl::ashape::rvec_mat_mult>
+{
+    typedef vector::rvec_mat_times_expr<Op1, Op2> type;
+};
+
+/// Scale row vector from right
+// added by Hui Li
+template <typename Op1, typename Op2>
+struct vec_mult_result_aux<Op1, Op2, ::mtl::ashape::rvec_scal_mult>
+{
+    typedef vector::rscaled_view<Op1, Op2> type;
+};
+
+/// Scale column vector from right
+// added by Hui Li
+template <typename Op1, typename Op2>
+struct vec_mult_result_aux<Op1, Op2, ::mtl::ashape::cvec_scal_mult>
+{
+    typedef vector::rscaled_view<Op1, Op2> type;
+};
+
+
+/// Enabler if operation is rvec_cvec_mult
+template <typename Op1, typename Op2, typename Result>
+struct lazy_enable_if_rvec_cvec_mult
+  : boost::lazy_enable_if<boost::is_same<typename ashape::mult_op<typename ashape::ashape<Op1>::type,
+								  typename ashape::ashape<Op2>::type >::type,
+					 ashape::rvec_cvec_mult>,
+			  Result>
+{};
+
+
+}} // namespace mtl::traits
+
+#endif // MTL_MULT_RESULT_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/mult_specialize.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef MTL_MULT_SPECIALIZE_INCLUDE
+#define MTL_MULT_SPECIALIZE_INCLUDE
+
+#include <boost/numeric/mtl/operation/dmat_dmat_mult.hpp>
+#include <boost/mpl/bool.hpp>
+
+namespace mtl { namespace matrix {namespace detail {
+
+template <typename MatrixA, typename MatrixB, typename MatrixC>
+struct dmat_dmat_mult_tiling1
+{
+    static const unsigned long value= 2;
+};
+
+template <typename MatrixA, typename MatrixB, typename MatrixC>
+struct dmat_dmat_mult_tiling2
+{
+    static const unsigned long value= 4;
+};
+
+template <typename MatrixA, typename MatrixB, typename MatrixC>
+struct dmat_dmat_mult_specialize
+    : public boost::mpl::false_
+{};
+
+/*
+   In order to specialize the functor, write for instance:
+
+template <typename MatrixA, typename MatrixB, typename MatrixC>
+struct dmat_dmat_mult_specialize
+    : public boost::mpl::true_
+{
+    typedef gen_dmat_dmat_mult_t<> type;
+};
+*/
+
+
+
+
+}}} // namespace mtl::matrix::detail
+
+#endif // MTL_MULT_SPECIALIZE_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/multi_action_block.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef MTL_MULTI_ACTION_BLOCK_INCLUDE
+#define MTL_MULTI_ACTION_BLOCK_INCLUDE
+
+namespace mtl {
+
+/*
+     Functor for unrolling arbitrary loops with independent operations
+     ...
+
+ */
+template <typename MultiAction, unsigned Steps> struct multi_action_block;
+
+
+template <typename MultiAction, unsigned MaxSteps, unsigned RemainingSteps>
+struct multi_action_helper
+{
+    static unsigned const step= MaxSteps - RemainingSteps;
+
+    void operator() (MultiAction const& action) const
+    {
+	action(step);
+	multi_action_helper<MultiAction, MaxSteps, RemainingSteps-1>()(action);
+    }
+
+    void operator() (MultiAction& action) const
+    {
+	action(step);
+	multi_action_helper<MultiAction, MaxSteps, RemainingSteps-1>()(action);
+    }
+};
+
+
+template <typename MultiAction, unsigned MaxSteps>
+struct multi_action_helper<MultiAction, MaxSteps, 1>
+{
+    static unsigned const step= MaxSteps - 1;
+
+    void operator() (MultiAction const& action) const
+    {
+	action(step);
+    }
+
+    void operator() (MultiAction& action) const
+    {
+	action(step);
+    }
+};
+
+
+template <typename MultiAction, unsigned Steps>
+struct multi_action_block
+{
+    void operator() (MultiAction const& action) const
+    {
+	multi_action_helper<MultiAction, Steps, Steps>()(action);
+    }
+
+    void operator() (MultiAction& action) const
+    {
+	multi_action_helper<MultiAction, Steps, Steps>()(action);
+    }
+};
+
+
+
+} // namespace mtl
+
+#endif // MTL_MULTI_ACTION_BLOCK_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/no_op.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef MTL_NO_OP_INCLUDE
+#define MTL_NO_OP_INCLUDE
+
+namespace mtl {
+
+struct no_op {};
+
+} // namespace mtl
+
+#endif // MTL_NO_OP_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/norms.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef MTL_NORMS_INCLUDE
+#define MTL_NORMS_INCLUDE
+
+#include <boost/numeric/mtl/operation/one_norm.hpp>
+#include <boost/numeric/mtl/operation/two_norm.hpp>
+#include <boost/numeric/mtl/operation/infinity_norm.hpp>
+#include <boost/numeric/mtl/operation/frobenius_norm.hpp>
+
+#endif // MTL_NORMS_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/num_cols.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef MTL_NUM_COLS_INCLUDE
+#define MTL_NUM_COLS_INCLUDE
+
+#include <vector>
+
+namespace mtl {
+
+namespace traits {
+
+    /// General declaration, used to disable unsupported types
+    template <typename Collection>
+    struct num_cols {};
+
+    /// num_cols implementation for STL vectors
+    template <typename Value>
+    struct num_cols< std::vector<Value> >
+    {
+	typedef std::size_t   type;
+	type operator()(const std::vector<Value>& ) { return 1; }
+    };
+
+    /// num_cols implementation for (1D) arrays interpreted as vectors
+    template <typename Value, unsigned Size>
+    struct num_cols<Value[Size]>
+    {
+	typedef std::size_t   type;
+	type operator()(const Value[Size]) { return 1; }
+    };
+
+    /// num_cols implementation for (2D and higher) arrays interpreted as matrices
+    template <typename Value, unsigned Rows, unsigned Cols>
+    struct num_cols<Value[Rows][Cols]>
+    {
+	typedef std::size_t   type;
+	type operator()(const Value[Rows][Cols]) { return Cols; }
+    };
+}
+
+
+/// num_cols function for non-MTL types (uses implicit enable_if), 1D interpreted as Column vector
+template <typename Collection>
+typename traits::num_cols<Collection>::type
+inline num_cols(const Collection& c)
+{
+    return traits::num_cols<Collection>()(c);
+}
+
+
+} // namespace mtl
+
+#endif // MTL_NUM_COLS_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/num_rows.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef MTL_NUM_ROWS_INCLUDE
+#define MTL_NUM_ROWS_INCLUDE
+
+#include <vector>
+
+namespace mtl {
+
+namespace traits {
+
+    /// General declaration, used to disable unsupported types
+    template <typename Collection>
+    struct num_rows {};
+
+    /// num_rows implementation for STL vectors
+    template <typename Value>
+    struct num_rows< std::vector<Value> >
+    {
+	typedef std::size_t   type;
+	type operator()(const std::vector<Value>& v) { return v.size(); }
+    };
+
+    /// num_rows implementation for (1D) arrays interpreted as vectors
+    template <typename Value, unsigned Size>
+    struct num_rows<Value[Size]>
+    {
+	typedef std::size_t   type;
+	type operator()(const Value[Size]) { return Size; }
+    };
+
+    /// num_rows implementation for (2D and higher) arrays interpreted as matrices
+    template <typename Value, unsigned Rows, unsigned Cols>
+    struct num_rows<Value[Rows][Cols]>
+    {
+	typedef std::size_t   type;
+	type operator()(const Value[Rows][Cols]) { return Rows; }
+    };
+}
+
+
+/// num_rows function for non-MTL types (uses implicit enable_if), 1D interpreted as Column vector
+template <typename Collection>
+typename traits::num_rows<Collection>::type
+inline num_rows(const Collection& c)
+{
+    return traits::num_rows<Collection>()(c);
+}
+
+
+} // namespace mtl
+
+#endif // MTL_NUM_ROWS_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/one_norm.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef MTL_ONE_NORM_INCLUDE
+#define MTL_ONE_NORM_INCLUDE
+
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/utility/enable_if.hpp>
+#include <boost/numeric/mtl/utility/is_row_major.hpp>
+#include <boost/numeric/mtl/utility/tag.hpp>
+#include <boost/numeric/mtl/utility/category.hpp>
+#include <boost/numeric/mtl/utility/property_map.hpp>
+#include <boost/numeric/mtl/operation/max_of_sums.hpp>
+#include <boost/numeric/mtl/vector/lazy_reduction.hpp>
+#include <boost/numeric/mtl/vector/reduction.hpp>
+#include <boost/numeric/mtl/vector/reduction_functors.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+
+namespace mtl {
+
+    namespace vector {
+
+	template <unsigned long Unroll, typename Vector>
+	typename traits::enable_if_vector<Vector, typename RealMagnitude<typename Collection<Vector>::value_type>::type>::type
+	inline one_norm(const Vector& vector)
+	{
+	    vampir_trace<2014> tracer;
+	    typedef typename RealMagnitude<typename Collection<Vector>::value_type>::type result_type;
+	    return reduction<Unroll, one_norm_functor, result_type>::apply(vector);
+	}
+
+	/*! One-norm for vectors: one_norm(x) \f$\rightarrow |x|_1\f$.
+	    \retval The magnitude type of the respective value type, see Magnitude.
+
+	    The norms are defined as \f$|v|_1=\sum_i |v_i|\f$.
+	    Vector norms are unrolled 8-fold by default.
+	    An n-fold unrolling can be generated with one_norm<n>(x).
+	    The maximum for n is 8 (it might be increased later).
+	**/
+	template <typename Vector>
+	typename traits::enable_if_vector<Vector, typename RealMagnitude<typename Collection<Vector>::value_type>::type>::type
+	inline one_norm(const Vector& vector)
+	{
+	    return one_norm<8>(vector);
+	}
+
+	template <typename Vector>
+	lazy_reduction<Vector, one_norm_functor> inline lazy_one_norm(const Vector& v)
+	{  return lazy_reduction<Vector, one_norm_functor>(v); 	}
+    }
+
+    namespace matrix {
+
+	// Ignore unrolling for matrices
+	template <unsigned long Unroll, typename Matrix>
+	typename mtl::traits::enable_if_matrix<Matrix, typename RealMagnitude<typename Collection<Matrix>::value_type>::type>::type
+	inline one_norm(const Matrix& matrix)
+	{
+	    vampir_trace<3025> tracer;
+	    using mtl::impl::max_of_sums;
+	    typename mtl::traits::col<Matrix>::type                             col(matrix);
+	    return max_of_sums(matrix, !mtl::traits::is_row_major<typename OrientedCollection<Matrix>::orientation>(),
+			       col, num_cols(matrix));
+	}
+
+	/*! One-norm for matrices: one_norm(x) \f$\rightarrow |x|_1\f$.
+	    \retval The magnitude type of the respective value type, see Magnitude.
+
+	    The norms are defined as \f$|A|_1=\max_j\{\sum_i(|A_{ij}|)\}\f$.
+	    Matrix norms are not optimized by unrolling (yet).
+	**/
+	template <typename Matrix>
+	typename mtl::traits::enable_if_matrix<Matrix, typename RealMagnitude<typename Collection<Matrix>::value_type>::type>::type
+	inline one_norm(const Matrix& matrix)
+	{
+	    return one_norm<8>(matrix);
+	}
+    }
+
+    using vector::one_norm;
+    using vector::lazy_one_norm;
+    using matrix::one_norm;
+
+} // namespace mtl
+
+#endif // MTL_ONE_NORM_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/ones.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef MTL_ONES_INCLUDE
+#define MTL_ONES_INCLUDE
+
+#include <boost/numeric/mtl/matrix/implicit_dense.hpp>
+
+namespace mtl {
+
+    namespace matrix {
+
+	/// Return r by c matrix with ones of type Value in all entries
+	template <typename Value>
+	ones_matrix<Value> inline ones(std::size_t r, std::size_t c)
+	{
+	    return ones_matrix<Value>(r, c);
+	}
+
+	/// Return r by c matrix with ones of type Value in all entries
+	ones_matrix<> inline ones(std::size_t r, std::size_t c)
+	{
+	    return ones_matrix<>(r, c);
+	}
+
+    }
+
+    using matrix::ones;
+
+
+} // namespace mtl
+
+#endif // MTL_ONES_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/operators.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef MTL_OPERATORS_INCLUDE
+#define MTL_OPERATORS_INCLUDE
+
+#include <boost/utility/enable_if.hpp>
+#include <boost/mpl/and.hpp>
+
+#include <boost/numeric/mtl/utility/ashape.hpp>
+#include <boost/numeric/mtl/matrix/operators.hpp>
+//#include <boost/numeric/mtl/vector/operators.hpp>
+#include <boost/numeric/mtl/operation/mult_result.hpp>
+#include <boost/numeric/mtl/operation/div_result.hpp>
+#include <boost/numeric/mtl/operation/dot.hpp>
+#include <boost/numeric/mtl/operation/mat_cvec_times_expr.hpp>
+#include <boost/numeric/mtl/matrix/all_mat_expr.hpp>
+#include <boost/numeric/mtl/utility/enable_if.hpp>
+#include <boost/numeric/mtl/utility/category.hpp>
+#include <boost/numeric/mtl/utility/true_copy.hpp>
+#include <boost/numeric/mtl/vector/rvec_mat_times_expr.hpp>
+
+
+namespace mtl {
+
+namespace matrix {
+
+    /// Multiplication for all supported types of operations
+    /** Enable-if-like technique make sure that only called when properly defined **/
+    template <typename Op1, typename Op2>
+    typename mtl::traits::mult_result<typename mtl::traits::true_copy<Op1>::type, Op2>::type
+    inline operator*(const Op1& op1, const Op2& op2)
+    {
+
+        return typename mtl::traits::mult_result<typename mtl::traits::true_copy<Op1>::type, Op2>::type(op1, op2);
+    }
+
+
+
+    /// Division of matrices and vectors by salars
+    /** Enable-if-like technique make sure that only called when properly defined **/
+    // added by Hui Li
+    // enable_if_matrix shouldn't be needed but was nessecary in cuppen.hpp
+    template < typename Op1, typename Op2 >
+    typename mtl::traits::enable_if_matrix<Op1, typename mtl::traits::div_result<Op1,Op2>::type>::type
+    inline operator/(const Op1& op1, const Op2& op2)
+    {
+        return typename mtl::traits::div_result<Op1,Op2>::type(op1,op2);
+    }
+
+} // namespace matrix
+
+
+namespace vector {
+
+    /// Multiplication for all supported types of operations
+    /** Enable-if-like technique make sure that only called when properly defined **/
+    template <typename Op1, typename Op2>
+    typename mtl::traits::vec_mult_result<typename mtl::traits::true_copy<Op1>::type, Op2>::type
+    inline operator*(const Op1& op1, const Op2& op2)
+    {
+        return typename mtl::traits::vec_mult_result<typename mtl::traits::true_copy<Op1>::type, Op2>::type(op1, op2);
+    }
+
+    /// Multiply row vector with column vector; result is scalar
+    template <typename Op1, typename Op2>
+    typename traits::lazy_enable_if_rvec_cvec_mult<Op1, Op2, detail::dot_result<Op1, Op2> >::type
+    inline operator*(const Op1& op1, const Op2& op2)
+    {
+	return dot_real(op1, op2);
+    }
+
+    /// Division of matrices and vectors by salars
+    /** Enable-if-like technique make sure that only called when properly defined **/
+    // added by Hui Li
+    template < typename Op1, typename Op2 >
+    typename traits::div_result<Op1,Op2>::type
+    inline operator/(const Op1& op1, const Op2& op2)
+    {
+        return typename traits::div_result<Op1,Op2>::type(op1,op2);
+    }
+
+    /// Compare two vectors for equality
+    /** Enable-if makes sure that only called when properly defined **/
+    template < typename Op1, typename Op2 >
+    typename boost::enable_if<boost::mpl::and_<mtl::traits::is_vector<Op1>,
+					       mtl::traits::is_vector<Op2> >,
+			      bool>::type
+    inline operator==(const Op1& op1, const Op2& op2)
+    {
+	unsigned s1= num_rows(op1) * num_cols(op1), s2= num_rows(op2) * num_cols(op2); // ugly hack to fight with ADL
+	if (s1 != s2)
+	    return false;
+	for (unsigned i= 0; i < s1; i++)
+	    if (op1[i] != op2[i])
+		return false;
+	return true;
+    }
+
+    /// Compare two vectors for unequality
+    /** Enable-if makes sure that only called when properly defined **/
+    template < typename Op1, typename Op2 >
+    typename boost::enable_if<boost::mpl::and_<mtl::traits::is_vector<Op1>,
+					       mtl::traits::is_vector<Op2> >,
+			      bool>::type
+    inline operator!=(const Op1& op1, const Op2& op2)
+    {
+	return !(op1 == op2);
+    }
+
+} // namespace vector
+
+
+} // namespace mtl
+
+#endif // MTL_OPERATORS_INCLUDE