AMDiS/lib/mtl4/boost/numeric/mtl/operation/compute_factors.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_COMPUTE_FACTORS_INCLUDE
+#define MTL_COMPUTE_FACTORS_INCLUDE
+
+#include <boost/numeric/mtl/mtl_fwd.hpp>
+
+namespace mtl { namespace operation {
+
+// Default is to just refer to the expression
+template <typename Result, typename Expr>
+struct compute_one_factor
+{
+    typedef const Expr&           type;
+    typedef const Expr&           const_reference;
+
+    compute_one_factor(type src) : value(src) {}
+
+    type value;
+};
+
+template <typename Result, typename E1, typename E2>
+struct compute_one_factor<Result, matrix::mat_mat_times_expr<E1, E2> >
+{
+    typedef Result                type;
+    typedef const Result&         const_reference;
+
+    compute_one_factor(const matrix::mat_mat_times_expr<E1, E2>& src)
+	: value(src.first * src.second) {}
+
+    type value;
+};
+
+template <typename Result, typename E1, typename E2>
+struct compute_one_factor<Result, matrix::mat_mat_ele_times_expr<E1, E2> >
+{
+    typedef Result                type;
+    typedef const Result&         const_reference;
+
+    compute_one_factor(const matrix::mat_mat_ele_times_expr<E1, E2>& src)
+	: value(ele_prod(src.first, src.second)) {}
+
+    type value;
+};
+
+
+// Only defined for matrix::mat_mat_times_expr and matrix::mat_mat_ele_times_expr
+template <typename Result, typename Expr>
+struct compute_factors {};
+
+template <typename Result, typename E1, typename E2>
+struct compute_factors<Result, matrix::mat_mat_times_expr<E1, E2> >
+{
+    compute_factors(const matrix::mat_mat_times_expr<E1, E2>& src)
+	: first_factor(src.first), second_factor(src.second),
+	  first(first_factor.value), second(second_factor.value)
+    {}
+
+    compute_one_factor<Result, E1> first_factor;
+    compute_one_factor<Result, E2> second_factor;
+
+    typename compute_one_factor<Result, E1>::const_reference first;
+    typename compute_one_factor<Result, E2>::const_reference second;
+};
+
+// First factor is handled implicitly in the evaluation
+template <typename Result, typename E1, typename E2>
+struct compute_factors<Result, matrix::mat_mat_ele_times_expr<E1, E2> >
+{
+    compute_factors(const matrix::mat_mat_ele_times_expr<E1, E2>& src)
+	: first(src.first), second_factor(src.second),
+	  second(second_factor.value)
+    {}
+
+    const E1&                                                first;
+    compute_one_factor<Result, E2>                           second_factor;
+    typename compute_one_factor<Result, E2>::const_reference second;
+};
+
+}} // namespace mtl::operation
+
+#endif // MTL_COMPUTE_FACTORS_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/compute_summand.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__OPERATION_COMPUTE_SUMMAND_INCLUDE
+#define MTL__OPERATION_COMPUTE_SUMMAND_INCLUDE
+
+#include <boost/numeric/mtl/mtl_fwd.hpp>
+#include <boost/numeric/mtl/utility/copy_expression_const_ref_container.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+namespace mtl { namespace operation {
+
+/// Compute a summand in an expression
+/** For instance matrix vector products are transformed into mult function calls
+    when assigned to a column vector. Adding such a vector to other vector expressions
+    requires to compute (evaluate) the summand first and then add the resulting
+    vector. The current implementation assumes that the result of the operation
+    can be represented by the multiplied column vector. This is for instance wrong
+    when the matrix is complex and the vector real. To handle this is signifantly
+    more complicated and is planned for the future.
+**/
+template <typename Expr>
+struct compute_summand
+{
+    typedef Expr type;
+    compute_summand(const Expr& expr) : value(expr) {}
+    // value is a const& if Expr is a true vector and a copy if it is an expression
+    typename mtl::traits::copy_expression_const_ref_container<Expr>::type value;
+};
+
+
+/// Specialization for matrix vector products
+template <typename Matrix, typename CVector>
+struct compute_summand< mat_cvec_times_expr<Matrix, CVector> >
+{
+    typedef CVector    type;
+
+    compute_summand(const mat_cvec_times_expr<Matrix, CVector>& expr)
+      : value(num_rows(expr.first))
+#ifndef NDEBUG // might be helpful, e.g. for generating AST
+      , first(expr.first), second(expr.second)
+#endif
+    {
+	vampir_trace<3005> tracer;
+	value= expr.first * expr.second;
+    }
+
+    CVector value;
+#ifndef NDEBUG
+    const Matrix&  first;
+    const CVector& second;
+#endif
+};
+
+/// Specialization for matrix vector products
+template <typename Matrix, typename CVector>
+struct compute_summand< vector::mat_cvec_multiplier<Matrix, CVector> >
+{
+    typedef CVector    type;
+
+    compute_summand(const vector::mat_cvec_multiplier<Matrix, CVector>& expr)
+      : value(num_rows(expr.A))
+    {
+	vampir_trace<3005> tracer;
+	expr.assign_to(value);
+    }
+
+    CVector value;
+};
+
+}} // namespace mtl::operation
+
+#endif // MTL__OPERATION_COMPUTE_SUMMAND_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/conj.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_CONJ_INCLUDE
+#define MTL_CONJ_INCLUDE
+
+#include <complex>
+
+#include <boost/numeric/mtl/mtl_fwd.hpp>
+#include <boost/numeric/mtl/utility/enable_if.hpp>
+#include <boost/numeric/mtl/utility/tag.hpp>
+#include <boost/numeric/mtl/utility/category.hpp>
+#include <boost/numeric/mtl/utility/algebraic_category.hpp>
+#include <boost/numeric/mtl/utility/is_what.hpp>
+#include <boost/numeric/mtl/utility/view_code.hpp>
+#include <boost/numeric/mtl/utility/viewed_collection.hpp>
+#include <boost/numeric/mtl/utility/compose_view.hpp>
+#include <boost/numeric/linear_algebra/identity.hpp>
+#include <boost/numeric/mtl/matrix/view_ref.hpp>
+
+namespace mtl {
+
+namespace sfunctor {
+
+    template <typename Value, typename AlgebraicCategory>
+    struct conj_aux
+    {
+	typedef Value result_type;
+
+	static inline result_type apply(const Value& v)
+	{
+	    return v;
+	}
+
+	result_type operator() (const Value& v) const
+	{
+	    return v;
+	}
+    };
+
+
+    template <typename Value, typename AlgebraicCategory>
+    struct conj_aux<std::complex<Value>, AlgebraicCategory>
+    {
+	typedef std::complex<Value> result_type;
+
+	static inline result_type apply(const std::complex<Value>& v)
+	{
+	    return std::conj(v);
+	}
+
+	result_type operator() (const std::complex<Value>& v) const
+	{
+	    return std::conj(v);
+	}
+    };
+
+    // Only declarations here, definitions in matrix::map_view (vector::map_view)
+    template <typename Matrix>
+    struct conj_aux<Matrix, tag::matrix>;
+
+    template <typename Vector>
+    struct conj_aux<Vector, tag::vector>;
+
+    // Short cut for result type
+    template <typename Value>
+    struct conj
+	: public conj_aux<Value, typename mtl::traits::algebraic_category<Value>::type>
+    {};
+
+} // namespace sfunctor
+
+    namespace vector {
+
+	/// Conjugate of an vector
+	template <typename Vector>
+	typename mtl::traits::enable_if_vector<Vector, conj_view<Vector> >::type
+	inline conj(const Vector& v)
+	{
+	    return conj_view<Vector>(v);
+	}
+    }
+
+    namespace matrix {
+
+	namespace detail {
+
+	    template <typename Matrix>
+	    struct conj_trait
+	    {
+		static const unsigned code= mtl::traits::view_toggle_conj<mtl::traits::view_code<Matrix> >::value;
+		typedef typename mtl::traits::compose_view<code, typename mtl::traits::viewed_collection<Matrix>::type>::type type;
+
+		static inline type apply(const Matrix& A)
+		{
+		    return type(view_ref(A));
+		}
+	    };
+
+	}
+
+	/// Conjugate of a matrix
+	template <typename Matrix>
+	typename mtl::traits::lazy_enable_if_matrix<Matrix, detail::conj_trait<Matrix> >::type
+	inline conj(const Matrix& A)
+	{
+	    return detail::conj_trait<Matrix>::apply(A);
+	}
+    }
+
+    namespace scalar {
+
+	// Only scalar values remain here
+	template <typename Value>
+	typename mtl::traits::enable_if_scalar<
+	    Value
+	  , typename sfunctor::conj<Value>::result_type
+	>::type
+	inline conj(const Value& v)
+	{
+	    return mtl::sfunctor::conj<Value>::apply(v);
+	}
+
+	float inline conj(float v) { return v; }
+	double inline conj(double v) { return v; }
+	long double inline conj(long double v) { return v; }
+    }
+
+    /// Conjugate of vector, matrix, or scalar
+    using vector::conj;
+    using matrix::conj;
+    using scalar::conj;
+    // using std::conj;
+
+} // namespace mtl
+
+#endif // MTL_CONJ_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/copy.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_COPY_INCLUDE
+#define MTL_COPY_INCLUDE
+
+#include <boost/numeric/mtl/mtl_fwd.hpp>
+#include <boost/numeric/mtl/detail/index.hpp>
+#include <boost/numeric/mtl/utility/category.hpp>
+#include <boost/numeric/mtl/utility/flatcat.hpp>
+#include <boost/numeric/mtl/utility/tag.hpp>
+#include <boost/numeric/mtl/utility/is_row_major.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/mtl/utility/range_generator.hpp>
+#include <boost/numeric/mtl/utility/ashape.hpp>
+#include <boost/numeric/mtl/utility/property_map.hpp>
+#include <boost/numeric/mtl/utility/static_assert.hpp>
+#include <boost/numeric/mtl/utility/updater_to_assigner.hpp>
+#include <boost/numeric/mtl/matrix/inserter.hpp>
+#include <boost/numeric/mtl/operation/set_to_zero.hpp>
+#include <boost/numeric/mtl/operation/update.hpp>
+#include <boost/numeric/mtl/operation/print.hpp>
+#include <boost/numeric/mtl/operation/crop.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+#include <boost/type_traits/is_same.hpp>
+#include <boost/utility/enable_if.hpp>
+#include <iostream>
+#include <limits>
+
+namespace mtl {
+
+    namespace detail {
+
+	// Set Destination matrix to zero when source is sparse
+	// (otherwise everything is overwritten anyway)
+	template <typename MatrixDest>
+	inline void zero_with_sparse_src(MatrixDest& dest, tag::flat<tag::sparse>)
+	{
+	    set_to_zero(dest);
+	}
+
+	template <typename MatrixDest>
+	inline void zero_with_sparse_src(MatrixDest&, tag::universe) {}
+
+	// Adapt inserter size to operation
+	template <typename Updater> struct copy_inserter_size {};
+
+	// Specialization for store
+	template <typename Value>
+	struct copy_inserter_size< operations::update_store<Value> >
+	{
+	    template <typename MatrixSrc, typename MatrixDest>
+	    static inline int apply(const MatrixSrc& src, const MatrixDest& dest)
+	    {
+		// std::cout << "nnz = " << src.nnz() << ", dim1 = " << dest.dim1() << "\n";
+		return int(src.nnz() * 1.2 / dest.dim1());
+	    }
+	};
+
+	struct sum_of_sizes
+	{
+	    template <typename MatrixSrc, typename MatrixDest>
+	    static inline int apply(const MatrixSrc& src, const MatrixDest& dest)
+	    {	return int((src.nnz() + dest.nnz()) * 1.2 / dest.dim1()); }
+	};
+
+	// Specialization for plus and minus
+	template <typename Value> struct copy_inserter_size< operations::update_plus<Value> > : sum_of_sizes {};
+	template <typename Value> struct copy_inserter_size< operations::update_minus<Value> > : sum_of_sizes {};
+
+    } // namespace detail
+
+
+    template <typename Updater, typename MatrixSrc, typename MatrixDest>
+    inline void gen_matrix_copy(const MatrixSrc& src, MatrixDest& dest, bool with_reset)
+    {
+	vampir_trace<3002> tracer;
+	MTL_THROW_IF(num_rows(src) != num_rows(dest) || num_cols(src) != num_cols(dest), incompatible_size());
+
+	if (with_reset)
+	    detail::zero_with_sparse_src(dest, traits::sparsity_flatcat<MatrixSrc>());
+
+	typename traits::row<MatrixSrc>::type             row(src);
+	typename traits::col<MatrixSrc>::type             col(src);
+	typename traits::const_value<MatrixSrc>::type     value(src);
+	typedef typename traits::range_generator<tag::major, MatrixSrc>::type  cursor_type;
+
+	// std::cout << "Slot size is " << detail::copy_inserter_size<Updater>::apply(src, dest) << "\n";
+	matrix::inserter<MatrixDest, Updater>   ins(dest, detail::copy_inserter_size<Updater>::apply(src, dest));
+	for (cursor_type cursor = mtl::begin<tag::major>(src), cend = mtl::end<tag::major>(src);
+	     cursor != cend; ++cursor) {
+	    // std::cout << dest << '\n';
+
+	    typedef typename traits::range_generator<tag::nz, cursor_type>::type icursor_type;
+	    for (icursor_type icursor = mtl::begin<tag::nz>(cursor), icend = mtl::end<tag::nz>(cursor);
+		 icursor != icend; ++icursor) {
+		//std::cout << "in " << row(*icursor) << ", " << col(*icursor) << " insert " << value(*icursor) << '\n';
+		ins(row(*icursor), col(*icursor)) << value(*icursor); }
+	}
+    }
+
+    // Specialization for multi_vector
+    template <typename Updater, typename MatrixSrc, typename Vector>
+    inline void gen_matrix_copy(const MatrixSrc& src, mtl::matrix::multi_vector<Vector>& dest, bool)
+    {
+	MTL_THROW_IF(num_rows(src) != num_rows(dest) || num_cols(src) != num_cols(dest), incompatible_size());
+	typedef typename mtl::traits::updater_to_assigner<Updater>::type Assigner;
+
+	for (std::size_t i= 0, n= num_cols(src); i < n; ++i)
+	    Assigner::first_update(dest.vector(i), src.vector(i));
+    }
+
+    namespace {
+#    ifdef __clang__
+#     pragma clang diagnostic ignored "-Wunneeded-internal-declaration"
+#     pragma clang diagnostic ignored "-Wunused-function"
+#    endif
+	inline long inc_wo_over(long i) { return i == std::numeric_limits<long>::max() ? i : i+1; }
+	inline long negate_wo_over(long i) { return i == std::numeric_limits<long>::min() ? std::numeric_limits<long>::max() : -i; }
+    }
+
+    template <typename Updater, typename ValueSrc, typename Para, typename ValueDest>
+    typename boost::enable_if<boost::is_same<Updater, operations::update_store<ValueDest> > >::type
+    inline gen_matrix_copy(const matrix::banded_view<mtl::matrix::compressed2D<ValueSrc, Para> >& src, mtl::matrix::compressed2D<ValueDest, Para>& dest, bool)
+    {
+	vampir_trace<3061> tracer;
+	typedef typename Para::size_type size_type;
+	dest.change_dim(num_rows(src), num_cols(src)); // contains make_empty
+	set_to_zero(dest);
+	const mtl::matrix::compressed2D<ValueSrc, Para>  &sref= src.ref;
+	const std::vector<size_type>        &sstarts= sref.ref_major(), &sindices= sref.ref_minor();
+	long first, last;
+	if (traits::is_row_major<Para>::value) {
+	    first= src.get_begin();
+	    last= src.get_end();
+	} else {
+	    first= inc_wo_over(negate_wo_over(src.get_end()));
+	    last=  inc_wo_over(negate_wo_over(src.get_begin()));
+	}
+
+	long jd= 0, j_end= sstarts[0];
+	for (long i= 0, i_end= src.dim1(), f= first, l= last; i < i_end; ++i) {
+	    dest.ref_major()[i]= jd;
+	    long j= j_end;
+	    j_end= sstarts[i+1];
+	    while (j < j_end && long(sindices[j]) < f) j++;
+	    while (j < j_end && long(sindices[j]) < l) jd++, j++;
+	    f= inc_wo_over(f);
+	    l= inc_wo_over(l);
+	}
+	dest.ref_major()[src.dim1()]= jd;
+	dest.set_nnz(jd); // resizes indices and data
+
+	for (long i= 0, i_end= src.dim1(), jd= 0, j_end= sstarts[0]; i < i_end; ++i) {
+	    dest.ref_major()[i]= jd;
+	    long j= j_end;
+	    j_end= sstarts[i+1];
+	    while (j < j_end && long(sindices[j]) < first) j++;
+	    while (j < j_end && long(sindices[j]) < last) {
+		dest.ref_minor()[jd]= sindices[j];
+		dest.data[jd++]= sref.data[j++];
+	    }
+	    first= inc_wo_over(first);
+	    last= inc_wo_over(last);
+	}
+    }
+
+
+    /// Copy matrix \p src into matrix \p dest
+    template <typename MatrixSrc, typename MatrixDest>
+    inline void matrix_copy(const MatrixSrc& src, MatrixDest& dest)
+    {
+	gen_matrix_copy< operations::update_store<typename MatrixDest::value_type> >(src, dest, true);
+    }
+
+
+    /// Add matrix \p src to matrix \p dest in copy-like style
+    template <typename MatrixSrc, typename MatrixDest>
+    inline void matrix_copy_plus(const MatrixSrc& src, MatrixDest& dest)
+    {
+	gen_matrix_copy< operations::update_plus<typename MatrixDest::value_type> >(src, dest, false);
+    }
+
+    /// Subtract matrix \p src from matrix \p dest in copy-like style
+    template <typename MatrixSrc, typename MatrixDest>
+    inline void matrix_copy_minus(const MatrixSrc& src, MatrixDest& dest)
+    {
+	gen_matrix_copy< operations::update_minus<typename MatrixDest::value_type> >(src, dest, false);
+    }
+
+    /// Multiply matrix \p src element-wise with matrix \p dest in copy-like style
+    template <typename MatrixSrc, typename MatrixDest>
+    inline void matrix_copy_ele_times(const MatrixSrc& src, MatrixDest& dest)
+    {
+	vampir_trace<3001> tracer;
+	MTL_THROW_IF(num_rows(src) != num_rows(dest) || num_cols(src) != num_cols(dest), incompatible_size());
+
+	typename traits::row<MatrixDest>::type             row(dest);
+	typename traits::col<MatrixDest>::type             col(dest);
+	typename traits::value<MatrixDest>::type           value(dest);
+	typedef typename traits::range_generator<tag::major, MatrixDest>::type  cursor_type;
+	typedef typename traits::range_generator<tag::nz, cursor_type>::type icursor_type;
+
+	for (cursor_type cursor = begin<tag::major>(dest), cend = end<tag::major>(dest); cursor != cend; ++cursor)
+	    for (icursor_type icursor = begin<tag::nz>(cursor), icend = end<tag::nz>(cursor); icursor != icend; ++icursor)
+		value(*icursor, value(*icursor) * src[row(*icursor)][col(*icursor)]);
+#if 0   // copy would result in a*0 = a and 0*b = b!!!!
+	gen_matrix_copy< operations::update_times<typename MatrixDest::value_type> >(src, dest, false);
+#endif
+	crop(dest);
+    }
+
+
+    template <typename MatrixSrc, typename MatrixDest>
+    inline void copy(const MatrixSrc& src, tag::flat<tag::matrix>, MatrixDest& dest, tag::flat<tag::matrix>)
+	// inline void copy(const MatrixSrc& src, tag::matrix_expr, MatrixDest& dest, tag::matrix)
+    {
+	return matrix_copy(src, dest);
+    }
+
+
+
+    template <typename Updater, typename VectorSrc, typename VectorDest>
+    inline void gen_vector_copy(const VectorSrc& src, VectorDest& dest, bool with_reset)
+    {
+	// Works only with dense vectors as dest !!!!! (source could be sparse)
+	// Needs vector inserter
+	vampir_trace<2001> tracer;
+
+	MTL_STATIC_ASSERT((boost::is_same<typename ashape::ashape<VectorSrc>::type,
+					  typename ashape::ashape<VectorDest>::type>::value), "Source and target must have the same algebraic shape.");
+
+	MTL_THROW_IF(size(src) != size(dest), incompatible_size());
+
+	if (with_reset)
+	    detail::zero_with_sparse_src(dest, typename traits::category<VectorSrc>::type());
+
+	typename traits::index<VectorSrc>::type           index(src);
+	typename traits::const_value<VectorSrc>::type     value(src);
+
+	typedef typename traits::range_generator<tag::nz, VectorSrc>::type  cursor_type;
+	for (cursor_type cursor = begin<tag::nz>(src), cend = end<tag::nz>(src);
+	     cursor != cend; ++cursor)
+	    Updater()(dest[index(*cursor)], value(*cursor));
+    }
+
+    /// Copy vector \p src into vector \p dest
+    template <typename VectorSrc, typename VectorDest>
+    inline void vector_copy(const VectorSrc& src, VectorDest& dest)
+    {
+	gen_vector_copy< operations::update_store<typename VectorDest::value_type> >(src, dest, true);
+    }
+
+
+    /// Add vector \p src to vector \p dest in copy-like style
+    template <typename VectorSrc, typename VectorDest>
+    inline void vector_copy_plus(const VectorSrc& src, VectorDest& dest)
+    {
+	gen_vector_copy< operations::update_plus<typename VectorDest::value_type> >(src, dest, false);
+    }
+
+    /// Subtract vector \p src from vector \p dest in copy-like style
+    template <typename VectorSrc, typename VectorDest>
+    inline void vector_copy_minus(const VectorSrc& src, VectorDest& dest)
+    {
+	gen_vector_copy< operations::update_minus<typename VectorDest::value_type> >(src, dest, false);
+    }
+
+
+
+    template <typename VectorSrc, typename VectorDest>
+    inline void copy(const VectorSrc& src, tag::flat<tag::vector>, VectorDest& dest, tag::flat<tag::vector>)
+    {
+	return vector_copy(src, dest);
+    }
+
+
+    template <typename CollSrc, typename CollDest>
+    inline void copy(const CollSrc& src, CollDest& dest)
+    {
+	vampir_trace<3003> tracer;
+	return copy(src, traits::flatcat2<CollSrc, tag::matrix, tag::vector>(),
+		    dest, traits::flatcat2<CollDest, tag::matrix, tag::vector>());
+    }
+
+
+} // namespace mtl
+
+#endif // MTL_COPY_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/copysign.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_COPYSIGN_INCLUDE
+#define MTL_COPYSIGN_INCLUDE
+
+#include <complex>
+#include <cmath>
+#include <boost/numeric/mtl/operation/real.hpp>
+#include <boost/numeric/mtl/operation/imag.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+namespace mtl {
+
+namespace sfunctor {
+
+    template <typename Value1, typename Value2>
+    struct copysign
+    {
+	static inline Value1 apply(const Value1& v, const Value2& s)
+	{
+	    using math::zero; using std::abs;
+
+	    Value1 a(abs(v));
+	    return s < zero(s) ? -a : a;
+	}
+    };
+
+    // This specialization is questionable and thus subject to elimination
+    template <typename Value1, typename Value2>
+    struct copysign<std::complex<Value1>, Value2>
+    {
+	static inline Value1 apply(const std::complex<Value1>& v, const Value2& s)
+	{
+	    using mtl::real; using mtl::imag;
+	    return std::complex<Value1>(copysign<Value1, Value2>::apply(real(v), s),
+					copysign<Value1, Value2>::apply(imag(v), s));
+	}
+    };
+#if 0 // ndef _MSC_VER doesn't work on BigRed with gcc 4.2.2 (????) and isn't used anyway
+    template <>
+    struct copysign<float, float>
+    {
+	static inline float apply(float v, float s)
+	{ return ::copysignf(v, s); }
+    };
+
+    template <>
+    struct copysign<double, double>
+    {
+	static inline double apply(double v, double s)
+	{ return ::copysign(v, s); }
+    };
+
+    template <>
+    struct copysign<long double, long double>
+    {
+	static inline long double apply(long double v, long double s)
+	{ return copysignl(v, s); }
+    };
+#endif // _MSC_VER
+}
+
+/// sign of scalars; for complex numbers sign of real part
+template <typename Value1, typename Value2>
+inline Value1 copysign(const Value1& v, const Value2& s)
+{
+    vampir_trace<1001> tracer;
+    return sfunctor::copysign<Value1, Value2>::apply(v, s);
+}
+
+
+
+
+} // namespace mtl
+
+#endif // MTL_COPYSIGN_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/crop.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_CROP_INCLUDE
+#define MTL_CROP_INCLUDE
+
+#include <boost/numeric/mtl/utility/enable_if.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+namespace mtl {
+
+    namespace vector {
+
+	/// Remove all zero entries from a collection
+	/** Does nothing for dense collections **/
+	template <typename T>
+	typename mtl::traits::enable_if_vector<T, T&>::type inline crop(T& x)
+	{
+	    vampir_trace<3006> tracer;
+	    x.crop(); return x;
+	}
+    }
+
+    namespace matrix {
+
+	/// Remove all zero entries from a collection
+	/** Does nothing for dense collections **/
+	template <typename T>
+	typename mtl::traits::enable_if_matrix<T, T&>::type inline crop(T& x)
+	{
+	    vampir_trace<3006> tracer;
+	    x.crop(); return x;
+	}
+    }
+
+    using vector::crop;
+    using matrix::crop;
+
+} // namespace mtl
+
+#endif // MTL_CROP_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/cross.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_VECTOR_CROSS_INCLUDE
+#define MTL_VECTOR_CROSS_INCLUDE
+
+#include <boost/numeric/mtl/concept/std_concept.hpp>
+#include <boost/numeric/mtl/vector/dense_vector.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+namespace mtl { namespace vector {
+
+    namespace detail {
+
+	// Result type of cross product
+	template <typename Vector1, typename Vector2>
+	struct cross_result
+	{
+	    typedef typename Multiplicable<typename Collection<Vector1>::value_type,
+					   typename Collection<Vector2>::value_type>::result_type value;
+	    typedef dense_vector<value>  type;
+	};
+    }
+
+/// Cross product
+/** Only exists for 3 and 7 dimensions.
+    Consider specialization for fixed-size types
+ **/
+template <typename Vector1, typename Vector2>
+typename detail::cross_result<Vector1, Vector2>::type
+inline cross(const Vector1& v1, const Vector2& v2)
+{
+    vampir_trace<2002> tracer;
+    MTL_THROW_IF((size(v1) != 3 && size(v1) != 7 ) || size(v1) != size(v2), incompatible_size());
+
+    typename detail::cross_result<Vector1, Vector2>::type result(size(v1));
+
+    if (size(v1) == 3)
+	for (unsigned i= 0; i < 3; i++) {
+	    unsigned k= (i+1) % 3, l= (i+2) % 3;
+	    result[i]= v1[k] * v2[l] - v1[l] * v2[k];
+	}
+    else  // must be 7 thus
+	for (unsigned i= 0; i < 7; i++) {
+	    unsigned k= (i+1) % 7, l= (i+3) % 7;
+	    result[i]= v1[k] * v2[l] - v1[l] * v2[k];
+
+	    k= (i+2) % 7, l= (i+6) % 7;
+	    result[i]+= v1[k] * v2[l] - v1[l] * v2[k];
+
+	    k= (i+4) % 7, l= (i+5) % 7;
+	    result[i]+= v1[k] * v2[l] - v1[l] * v2[k];
+	}
+    return result;
+}
+
+
+}} // namespace mtl::vector
+
+#endif // MTL_VECTOR_CROSS_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/cuppen.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.
+
+// With contributions from Cornelius Steinhardt
+
+#ifndef MTL_MATRIX_CUPPEN_INCLUDE
+#define MTL_MATRIX_CUPPEN_INCLUDE
+
+#include <cmath>
+#include <boost/numeric/linear_algebra/identity.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/operation/iota.hpp>
+#include <boost/numeric/mtl/operation/secular.hpp>
+#include <boost/numeric/mtl/operation/sort.hpp>
+#include <boost/numeric/mtl/operation/trans.hpp>
+#include <boost/numeric/mtl/utility/domain.hpp>
+#include <boost/numeric/mtl/matrix/permutation.hpp>
+
+#include <boost/numeric/mtl/vector/dense_vector.hpp>
+#include <boost/numeric/itl/iteration/basic_iteration.hpp>
+#include <boost/numeric/itl/krylov/fsm.hpp>
+
+namespace mtl { namespace matrix {
+
+/// Eigenvalues of triangle matrix A with Cuppen's divide and conquer algorithm
+/** Eigenvalues are returned in vector lambda. A is overwritten. **/
+template <typename Matrix, typename Vector>
+void inline cuppen_inplace(Matrix& A, Matrix& Q, Vector& lambda)
+{
+    using std::abs; using mtl::irange; using mtl::imax; using mtl::iall; using vector::iota;
+
+    typedef typename Collection<Matrix>::value_type     value_type;
+    typedef typename Collection<Matrix>::size_type      size_type;
+    typedef mtl::vector::dense_vector<size_type, vector::parameters<> >   size_vector; // todo: with type trait
+
+    size_type        nrows= num_rows(A);
+    MTL_THROW_IF(nrows != num_cols(A), matrix_not_square());
+    const value_type zero= 0, one= 1;
+
+    if (nrows == 1){
+	lambda[0]= A[0][0];
+	Q= one;
+    } else {
+	size_type     m= size_type(nrows/2);
+	irange        till_m(m), from_m(m, imax);
+
+	size_vector   perm(nrows);
+	Matrix        T1(A[till_m][till_m]), T2(A[from_m][from_m]),                               // sub-matrices of A
+	              Q0(nrows, nrows), Q1(Q0[till_m][till_m]), Q2(Q0[from_m][from_m]);           // Q0 and sub-matrices
+	Vector        v(nrows, zero), diag(nrows), lambda1(diag[till_m]), lambda2(diag[from_m]);  // sub-vectors of diag
+
+	//DIVIDE
+	value_type    b= A[m-1][m];
+	T1[m-1][m-1]-= abs(b);
+	T2[0][0]-= abs(b);
+
+	v[m-1]= b > zero ? one : -one;
+	v[m]= one;
+
+	cuppen_inplace(T1, Q1, lambda1);
+	cuppen_inplace(T2, Q2, lambda2);
+
+	Q0[till_m][from_m]= zero; Q0[from_m][till_m]= zero; // zero out non-diagonal blocks
+
+	T1[m-1][m-1]+= abs(b);
+	T2[0][0]+= abs(b);
+
+	iota(perm);
+	sort(diag, perm);
+
+	// CONQUER, start with eq. (3.0.2) using rows (not columns)
+	v[till_m]= b < zero ? Vector(-trans(Q1[m-1][iall])) : trans(Q1[m-1][iall]);
+	v[from_m]= trans(Q2[0][iall]);
+
+	// permutation on v
+	mtl::matrix::traits::permutation<>::type P= mtl::matrix::permutation(perm);
+	Vector v1(P * v);
+
+	lambda= secular(v1, diag, abs(b));   // solve secular equation
+
+	// std::cout << "lambda is " << lambda << "\ndiag is " << diag << '\n';
+	//Lemma 3.0.2  ... calculate eigenvectors
+	Matrix Q_tilde(nrows, nrows);
+	for (size_type i = 0; i < nrows; i++) {
+	    for (size_type j= 0; j < size(diag); ++j)
+		MTL_DEBUG_THROW_IF (diag[j] == lambda[i], logic_error("Can't compute eigenvector, probably due to double eigenvalue."));
+	    Vector    li(nrows, lambda[i]), lambda_i(ele_quot(v1, diag - li));
+	    Q_tilde[iall][i]= lambda_i / two_norm(lambda_i); // normalized eigenvector in Matrix Q
+	}
+
+	Q= Q0 * P * Q_tilde;
+
+#if 0
+	for (size_type i = 0; i < nrows; i++) {
+	    // Vector    qi(Q[iall][i]); // Todo: find out valgrind complains about the memory of qi for reasons inexplicable
+	    Vector qi(nrows);
+	    for (size_type j= 0; j < nrows; j++)
+		qi[j]= Q[j][i];
+
+	    std::cout << "q[" << i << "] = " << qi << ", lambda[i] = " << lambda[i]
+		      << ", diff = " << two_norm(Vector(A*qi - lambda[i]*qi)) << ", A*qi = " << Vector(A*qi) << ", li*qi = " << Vector(lambda[i]*qi) << '\n';
+	    itl::basic_iteration<double>   iter(1.0, 20, 1e-5, 1e-5);
+	    fsm(A, qi, lambda[i], 0.1, iter);
+	    std::cout << "q[" << i << "] = " << qi << ", diff = " << two_norm(Vector(A*qi - lambda[i]*qi)) << '\n';
+	    Q[iall][i]= qi;
+	}
+#endif
+    }
+}
+
+/// Eigenvalues of triangle matrix A with Cuppen's divide and conquer algorithm
+/** Eigenvalues are returned in vector lambda. A is copied. **/
+// A not as reference to force copy
+template <typename Matrix, typename Vector>
+void inline cuppen(Matrix A, Matrix& Q, Vector& lambda)
+{
+    Q= 0.0;
+    cuppen_inplace(A, Q, lambda);
+}
+
+}} // namespace mtl::matrix
+
+#endif // MTL_MATRIX_CUPPEN_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/cursor_pseudo_dot.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_CURSOR_PSEUDO_DOT_INCLUDE
+#define MTL_CURSOR_PSEUDO_DOT_INCLUDE
+
+namespace mtl {
+
+
+
+namespace functor {
+
+    template <unsigned MaxDepth, typename Value, typename Cursor1, typename Prop1, typename Cursor2, typename Prop2, unsigned Depth>
+    struct cursor_pseudo_dot_block
+    {
+	static unsigned const offset= MaxDepth - Depth;
+
+	void operator() (Cursor1 i1, Prop1& prop1, Cursor2 i2, Prop2& prop2,
+			 Value& s0, Value& s1, Value& s2, Value& s3,
+			 Value& s4, Value& s5, Value& s6, Value& s7)
+	{
+	    Cursor1 tmp1(i1); tmp1+= offset;
+	    Cursor2 tmp2(i2); tmp2+= offset;
+	    s0+= prop1(*tmp1) * prop2(*tmp2);
+	    // s0+= prop1(i1 + offset) * prop2(i2 + offset);
+	    typedef cursor_pseudo_dot_block<MaxDepth, Value, Cursor1, Prop1, Cursor2, Prop2, Depth-1> block_rest;
+	    block_rest() (i1, prop1, i2, prop2, s1, s2, s3, s4, s5, s6, s7, s0);
+	}
+    };
+
+    //template <>
+    template <unsigned MaxDepth, typename Value, typename Cursor1, typename Prop1, typename Cursor2, typename Prop2>
+    struct cursor_pseudo_dot_block<MaxDepth, Value, Cursor1, Prop1, Cursor2, Prop2, 1>
+    {
+	static unsigned const offset= MaxDepth - 1;
+
+	void operator() (Cursor1 i1, Prop1& prop1, Cursor2 i2, Prop2& prop2,
+			 Value& s0, Value&, Value&, Value&,
+			 Value&, Value&, Value&, Value&)
+	{
+	    s0+= prop1(*(i1 + offset)) * prop2(*(i2 + offset));
+	}
+    };
+
+    template <unsigned MaxDepth, typename Value, typename Cursor1, typename Prop1, typename Cursor2, typename Prop2>
+    struct cursor_pseudo_dot_t
+    {
+	Value operator() (Cursor1 i1, Cursor1 end1, Prop1& prop1, Cursor2 i2, Prop2& prop2)
+	{
+	    using math::zero;
+	    Value         ref, my_zero(zero(ref)),
+                          s0= my_zero, s1= my_zero, s2= my_zero, s3= my_zero,
+		          s4= my_zero, s5= my_zero, s6= my_zero, s7= my_zero;
+	    std::size_t size= end1 - i1, blocks= size / MaxDepth, blocked_size= blocks * MaxDepth;
+
+	    typedef cursor_pseudo_dot_block<MaxDepth, Value, Cursor1, Prop1, Cursor2, Prop2, MaxDepth> dot_block_type;
+	    for (unsigned i= 0; i < blocked_size; i+= MaxDepth, i1+= MaxDepth, i2+= MaxDepth) {
+		dot_block_type()(i1, prop1, i2, prop2, s0, s1, s2, s3, s4, s5, s6, s7);
+	    }
+
+	    typedef cursor_pseudo_dot_block<MaxDepth, Value, Cursor1, Prop1, Cursor2, Prop2, MaxDepth> dot_single_type;
+	    s0+= s1 + s2 + s3 + s4 + s5 + s6 + s7;
+	    for (unsigned i= blocked_size; i < size; ++i, ++i1, ++i2)
+		dot_single_type()(i1, prop1, i2, prop2, s0, s1, s2, s3, s4, s5, s6, s7);
+	    return s0;
+	}
+    };
+
+} // namespace functor
+
+template <unsigned MaxDepth, typename Value, typename Cursor1, typename Prop1, typename Cursor2, typename Prop2>
+Value cursor_pseudo_dot(Cursor1 i1, Cursor1 end1, Prop1 prop1, Cursor2 i2, Prop2 prop2, Value)
+{
+    return functor::cursor_pseudo_dot_t<MaxDepth, Value, Cursor1, Prop1, Cursor2, Prop2>()(i1, end1, prop1, i2, prop2);
+}
+
+
+
+} // namespace mtl
+
+#endif // MTL_CURSOR_PSEUDO_DOT_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/diagonal.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_DIAGONAL_INCLUDE
+#define MTL_DIAGONAL_INCLUDE
+
+
+#include <boost/numeric/mtl/utility/enable_if.hpp>
+#include <boost/numeric/mtl/matrix/parameter.hpp>
+#include <boost/numeric/mtl/matrix/compressed2D.hpp>
+#include <boost/numeric/mtl/matrix/inserter.hpp>
+#include <boost/numeric/mtl/vector/dense_vector.hpp>
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+namespace mtl {
+
+    namespace vector {
+
+	/// Transform a vector into a diagonal matrix
+	template <typename Vector>
+	mtl::matrix::compressed2D<typename Collection<Vector>::value_type, matrix::parameters<> >
+	// typename mtl::traits::enable_if_vector<Vector, mtl::matrix::compressed2D<typename Collection<Vector>::value_type> >::type
+	inline diagonal(const Vector& v)
+	{
+	    vampir_trace<2016> tracer;
+	    typedef mtl::matrix::compressed2D<typename Collection<Vector>::value_type, matrix::parameters<> > matrix_type;
+	    matrix_type                           D(size(v), size(v));
+	    D= 0;
+	    mtl::matrix::inserter<matrix_type>    ins(D, 1);
+	    for (typename Collection<Vector>::size_type i= 0; i < size(v); ++i)
+		ins[i][i] << v[i];
+
+	    return D;
+	}
+    }
+
+    namespace matrix {
+
+	/// Return the vector with the diagonal of the matrix
+	template <typename Matrix>
+	// typename mtl::traits::enable_if_matrix<Matrix, conj_view<Matrix> >::type
+	mtl::vector::dense_vector<typename Collection<Matrix>::value_type, vector::parameters<> >
+	inline diagonal(const Matrix& A)
+	{
+	    vampir_trace<3007> tracer;
+	    using std::min;
+	    typedef typename Collection<Matrix>::size_type size_type;
+	    size_type n= min(num_rows(A), num_cols(A));
+		mtl::vector::dense_vector<typename Collection<Matrix>::value_type, vector::parameters<> > v(n);
+
+	    for (size_type i= 0; i < n; ++i)
+		v[i]= A[i][i];
+	    return v;
+	}
+    }
+
+} // namespace mtl
+
+#endif // MTL_DIAGONAL_INCLUDE

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

+/*
+ *  div_result.h
+ *  MTL
+ *
+ *  Created by Hui Li (huil@Princeton.EDU)
+ *
+ */
+
+#ifndef MTL_DIV_RESULT_INCLUDE
+#define MTL_DIV_RESULT_INCLUDE
+
+#include <boost/numeric/mtl/utility/ashape.hpp>
+#include <boost/numeric/mtl/matrix/map_view.hpp>
+#include <boost/numeric/mtl/vector/map_view.hpp>
+
+namespace mtl { namespace traits {
+
+template < typename Op1, typename Op2, typename DivOp > struct div_result_aux {};
+
+/// Result type for dividing Op1 by Op2
+/** Can be used in enable-if-style as type is only defined when appropriate **/
+template < typename Op1, typename Op2 >
+struct div_result
+  : public div_result_aux < Op1, Op2, typename ashape::div_op<typename ashape::ashape<Op1>::type,
+							      typename ashape::ashape<Op2>::type >::type >
+{};
+
+/// Divide column vector by scalar
+template < typename Op1, typename Op2 >
+struct div_result_aux < Op1, Op2, ::mtl::ashape::cvec_scal_div >
+{
+    typedef typename vector::divide_by_view<Op1,Op2> type;
+};
+
+/// Divide row vector by scalar
+template < typename Op1, typename Op2 >
+struct div_result_aux < Op1, Op2, ::mtl::ashape::rvec_scal_div >
+{
+    typedef typename vector::divide_by_view<Op1,Op2> type;
+};
+
+/// Divide matrix by scalar
+template < typename Op1, typename Op2 >
+struct div_result_aux < Op1, Op2, ::mtl::ashape::mat_scal_div >
+{
+    typedef typename matrix::divide_by_view<Op1,Op2> type;
+};
+
+}} // namespace mtl::traits
+
+
+#endif

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

+/*
+ *  divide_by.h
+ *  MTL4
+ *
+ *  Created by Hui Li (huil@Princeton.EDU)
+ *
+ */
+
+#ifndef MTL_DIVIDE_BY_INCLUDE
+#define MTL_DIVIDE_BY_INCLUDE
+
+#include <boost/numeric/mtl/mtl_fwd.hpp>
+#include <boost/numeric/mtl/concept/std_concept.hpp>
+#include <boost/numeric/mtl/matrix/map_view.hpp>
+#include <boost/numeric/mtl/vector/map_view.hpp>
+#include <boost/numeric/mtl/utility/algebraic_category.hpp>
+
+namespace mtl { namespace tfunctor {
+
+    // AlgebraicCategory is by default tag::scalar
+    template <typename Value1, typename Value2, typename AlgebraicCategory>
+    struct divide_by
+    {
+	typedef typename Multiplicable<Value1, Value2>::result_type result_type;
+
+	explicit divide_by(const Value2& v2) : v2(v2) {}
+
+	result_type operator() (const Value1& v1) const
+	{
+	    return v1 / v2;
+	}
+    private:
+	Value2 v2;
+    };
+
+
+    template <typename Matrix, typename Value2>
+    struct divide_by<Matrix, Value2, tag::matrix>
+    {
+	typedef matrix::divide_by_view<Matrix,Value2> result_type;
+
+	explicit divide_by(const Value2& v2) : v2(v2) {}
+
+	result_type operator() (const Matrix& matrix) const
+	{
+	    return result_type(matrix, v2);
+	}
+    private:
+	Value2 v2;
+    };
+
+
+    template <typename Vector, typename Value2>
+    struct divide_by<Vector, Value2, tag::vector>
+    {
+	typedef vector::divide_by_view<Vector, Value2> result_type;
+
+	explicit divide_by(const Value2& v2) : v2(v2) {}
+
+	result_type operator() (const Vector& vector) const
+	{
+	    return result_type(vector, v2);
+	}
+    private:
+	Value2 v2;
+    };
+
+
+} // namespace tfunctor
+
+
+namespace matrix {
+
+    template <typename Value1, typename Value2>
+    typename tfunctor::divide_by<Value1, Value2, typename traits::algebraic_category<Value1>::type>::result_type
+    inline divide_by(const Value1& value1, const Value2& value2)
+    {
+	return tfunctor::divide_by<Value1, Value2, typename traits::algebraic_category<Value1>::type>(value2)(value1);
+    }
+}
+
+namespace vector {
+
+    template <typename Value1, typename Value2>
+    typename tfunctor::divide_by<Value1, Value2, typename traits::algebraic_category<Value1>::type>::result_type
+    inline divide_by(const Value1& value1, const Value2& value2)
+    {
+	return tfunctor::divide_by<Value1, Value2, typename traits::algebraic_category<Value1>::type>(value2)(value1);
+    }
+}
+
+} // namespace mtl
+
+#endif // MTL_DIVIDE_BY_INCLUDE

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

+/*
+ *  divide_by_inplace.hpp
+ *  MTL4
+ *
+ *  Created by Hui Li (huil@Princeton.EDU)
+ *
+ */
+
+#ifndef MTL_DIVIDE_BY_INPLACE_INCLUDE
+#define MTL_DIVIDE_BY_INPLACE_INCLUDE
+
+#include <boost/numeric/mtl/mtl_fwd.hpp>
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/utility/tag.hpp>
+#include <boost/numeric/mtl/operation/assign_each_nonzero.hpp>
+#include <boost/numeric/mtl/operation/divide_by.hpp>
+
+
+namespace mtl {
+
+
+	/// Divide collection \p c (from right) by scalar factor \p alpha; \p c is altered
+	template <typename Factor, typename Coll>
+	void divide_by_inplace(Coll& c, const Factor& alpha, tag::scalar)
+	{
+	    // assign_each_nonzero(c, boost::lambda::_1 / alpha);
+	    assign_each_nonzero(c, tfunctor::divide_by<typename Collection<Coll>::value_type, Factor>(alpha));
+	}
+
+	/// Divide collection \p c (from right) by factor \p alpha; \p c is altered
+	template <typename Factor, typename Collection>
+	void divide_by_inplace(Collection& c, const Factor& alpha)
+	{
+	    divide_by_inplace(c, alpha, typename traits::category<Factor>::type());
+	}
+
+
+} // namespace mtl
+
+#endif // MTL_DIVIDE_BY_INPLACE_INCLUDE