AMDiS/lib/mtl4/boost/numeric/mtl/operation/lazy_assign.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_LAZY_ASSIGN_INCLUDE
+#define MTL_LAZY_ASSIGN_INCLUDE
+
+namespace mtl {
+
+
+/// Helper class for lazy assignment semantics
+template <typename T, typename U, typename Assign>
+struct lazy_assign
+{
+    typedef Assign  assign_type;
+
+    lazy_assign(T& first, const U& second) : first(first), second(second) {}
+    void delay_assign() const {}
+
+    T&       first;
+    const U& second;
+
+};
+
+} // namespace mtl
+
+#endif // MTL_LAZY_ASSIGN_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/left_scale_inplace.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_LEFT_SCALE_INPLACE_INCLUDE
+#define MTL_LEFT_SCALE_INPLACE_INCLUDE
+
+#include <boost/numeric/mtl/mtl_fwd.hpp>
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/utility/category.hpp>
+#include <boost/numeric/mtl/utility/tag.hpp>
+#include <boost/numeric/mtl/utility/enable_if.hpp>
+#include <boost/numeric/mtl/operation/assign_each_nonzero.hpp>
+#include <boost/numeric/mtl/operation/mult.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+
+namespace mtl {
+
+    namespace impl {
+
+	template <typename Factor, typename Coll>
+	inline Coll& left_scale_inplace(const Factor& alpha, tag::scalar, Coll& c)
+	{
+	    vampir_trace<3014> tracer;
+	    assign_each_nonzero(c, tfunctor::scale<Factor, typename Collection<Coll>::value_type>(alpha));
+
+	    return c;
+	}
+
+	template <typename Factor, typename Collection>
+	inline Collection& left_scale_inplace(const Factor& alpha, tag::matrix, Collection& c)
+	{
+	    vampir_trace<3014> tracer;
+	    Collection tmp(alpha * c);
+	    swap(tmp, c);
+	    return c;
+	}
+    }
+
+    namespace matrix {
+
+	/// Scale matrix \p c from left with scalar or matrix factor \p alpha; \p c is altered
+	template <typename Factor, typename Matrix>
+	typename mtl::traits::enable_if_matrix<Matrix, Matrix&>::type
+	inline left_scale_inplace(const Factor& alpha, Matrix& A)
+	{
+	    return mtl::impl::left_scale_inplace(alpha, typename mtl::traits::category<Factor>::type(), A);
+	}
+    }
+
+    namespace vector {
+
+	/// Scale vector \p c from left with scalar or matrix factor \p alpha; \p c is altered
+	template <typename Factor, typename Vector>
+	typename mtl::traits::enable_if_vector<Vector, Vector&>::type
+	inline left_scale_inplace(const Factor& alpha, Vector& v)
+	{
+	    return mtl::impl::left_scale_inplace(alpha, typename mtl::traits::category<Factor>::type(), v);
+	}
+    }
+
+    using vector::left_scale_inplace;
+    using matrix::left_scale_inplace;
+
+} // namespace mtl
+
+#endif // MTL_LEFT_SCALE_INPLACE_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/look_at_each_nonzero.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_FOR_EACH_NONZERO_INCLUDE
+#define MTL_FOR_EACH_NONZERO_INCLUDE
+
+#include <utility>
+
+#include <boost/numeric/mtl/utility/tag.hpp>
+#include <boost/numeric/mtl/utility/enable_if.hpp>
+#include <boost/numeric/mtl/utility/range_generator.hpp>
+#include <boost/numeric/mtl/utility/property_map.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+
+namespace mtl {
+
+    namespace vector {
+
+	/// Perform \p f(v[i]) on each non-zero i in constant vector \p v; thus the must keep the result in its state
+	template <typename Vector, typename Functor>
+	inline void look_at_each_nonzero(const Vector& v, Functor& f)
+	{
+	    vampir_trace<2007> tracer;
+	    typedef typename mtl::traits::range_generator<tag::const_iter::nz, Vector>::type iterator;
+	    for (iterator i= begin<tag::iter::nz>(v), iend= end<tag::iter::nz>(v); i != iend; ++i)
+		f(*i);
+	}
+
+	/// Perform \p f(v[i], i) on each non-zero i in constant vector \p v; thus the must keep the result in its state
+	template <typename Vector, typename Functor>
+	typename mtl::traits::enable_if_vector<Vector>::type // to be called for vectors only
+	inline look_at_each_nonzero_pos(const Vector& v, Functor& f)
+	{
+	    vampir_trace<2008> tracer;
+	    typename mtl::traits::index<Vector>::type           index(v);
+	    typename mtl::traits::const_value<Vector>::type     value(v);
+
+	    typedef typename traits::range_generator<tag::nz, Vector>::type iterator;
+	    for (iterator i= begin<tag::nz>(v), iend= end<tag::nz>(v); i != iend; ++i)
+		f(value(*i), index(*i));
+	}
+
+    } // namespace vector
+
+    namespace matrix {
+
+	/// Perform a potentially mutating \p f(A[i][j]) on each non-zero entry in matrix \p A
+	template <typename Matrix, typename Functor>
+	inline void look_at_each_nonzero(const Matrix& A, Functor& f)
+	{
+	    vampir_trace<3016> tracer;
+	    typename mtl::traits::const_value<Matrix>::type     value(A);
+
+	    typedef typename mtl::traits::range_generator<tag::major, Matrix>::type     cursor_type;
+	    typedef typename mtl::traits::range_generator<tag::nz, cursor_type>::type   icursor_type;
+
+	    for (cursor_type cursor = begin<tag::major>(A), cend = end<tag::major>(A); cursor != cend; ++cursor)
+		for (icursor_type icursor = begin<tag::nz>(cursor), icend = end<tag::nz>(cursor);
+		     icursor != icend; ++icursor)
+		    f(value(*icursor));
+	}
+
+	/// Perform a potentially mutating \p f(A[i][j], make_pair(i, j)) on each non-zero entry in matrix \p A
+	template <typename Matrix, typename Functor>
+	typename mtl::traits::enable_if_matrix<Matrix>::type // to be called for matrices only
+	inline look_at_each_nonzero_pos(const Matrix& A, Functor& f)
+	{
+	    vampir_trace<3017> tracer;
+	    typename mtl::traits::row<Matrix>::type             row(A);
+	    typename mtl::traits::col<Matrix>::type             col(A);
+	    typename mtl::traits::const_value<Matrix>::type     value(A);
+
+	    typedef typename mtl::traits::range_generator<tag::major, Matrix>::type     cursor_type;
+	    typedef typename mtl::traits::range_generator<tag::nz, cursor_type>::type   icursor_type;
+
+	    for (cursor_type cursor = begin<tag::major>(A), cend = end<tag::major>(A); cursor != cend; ++cursor)
+		for (icursor_type icursor = begin<tag::nz>(cursor), icend = end<tag::nz>(cursor);
+		     icursor != icend; ++icursor)
+		    f(value(*icursor), std::make_pair(row(*icursor), col(*icursor)));
+	}
+
+    } // namespace matrix
+
+} // namespace mtl
+
+#endif // MTL_FOR_EACH_NONZERO_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/lower_trisolve.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_LOWER_TRISOLVE_INCLUDE
+#define MTL_LOWER_TRISOLVE_INCLUDE
+
+#include <boost/mpl/int.hpp>
+#include <boost/type_traits/is_same.hpp>
+#include <boost/numeric/mtl/utility/tag.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/mtl/utility/property_map.hpp>
+#include <boost/numeric/mtl/utility/range_generator.hpp>
+#include <boost/numeric/mtl/utility/category.hpp>
+#include <boost/numeric/mtl/utility/static_assert.hpp>
+#include <boost/numeric/mtl/concept/collection.hpp>
+
+#include <boost/numeric/linear_algebra/identity.hpp>
+#include <boost/numeric/linear_algebra/inverse.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+namespace mtl { namespace matrix {
+
+
+namespace detail {
+
+    /// Class that implements lower trisolver
+    /** DiaTag can be tag::regular_diagonal, tag::unit_diagonal, or tag::inverse_diagonal.
+	CompactStorage means that matrix contains only lower entries (strict lower when DiaTag == unit_diagonal). \sa \ref trisolve_object **/
+    template <typename Matrix, typename DiaTag, bool CompactStorage= false>
+    struct lower_trisolve_t
+    {
+	MTL_STATIC_ASSERT((boost::is_same<DiaTag, tag::regular_diagonal>::value
+			   || boost::is_same<DiaTag, tag::unit_diagonal>::value
+			   || boost::is_same<DiaTag, tag::inverse_diagonal>::value),
+			  "DiaTag must be either tag::regular_diagonal, tag::unit_diagonal, or tag::inverse_diagonal.");
+
+	typedef typename Collection<Matrix>::value_type         	       	    value_type;
+	typedef typename Collection<Matrix>::size_type          	       	    size_type;
+	typedef typename OrientedCollection<Matrix>::orientation	       	    my_orientation;
+	typedef typename mtl::traits::category<Matrix>::type    	       	    my_category;
+	typedef typename mtl::traits::range_generator<tag::major, Matrix>::type     a_cur_type; // row or col depending on Matrix
+	typedef typename mtl::traits::range_generator<tag::nz, a_cur_type>::type    a_icur_type;
+
+	/// Construction from matrix \p A
+	lower_trisolve_t(const Matrix& A) : A(A), value_a(A), col_a(A), row_a(A)
+	{    MTL_THROW_IF(num_rows(A) != num_cols(A), matrix_not_square());	}
+
+	template <typename M, typename D, bool C>
+	struct generic_version
+	  : boost::mpl::int_<(mtl::traits::is_row_major<M>::value ? 0 : 2)
+	                     + (boost::is_same<D, tag::unit_diagonal>::value ? 1 : 2)
+                            > {};
+
+	template <typename M, typename D, bool C>
+	struct version
+	  : generic_version<M, D, C> {};
+
+
+	template <typename Value, typename Para, typename D>
+	struct version<compressed2D<Value, Para>, D, true>
+	  : boost::mpl::if_<mtl::traits::is_row_major<Para>,
+			    typename boost::mpl::if_<boost::is_same<D, tag::unit_diagonal>,
+						     boost::mpl::int_<5>,
+						     boost::mpl::int_<6> >::type,
+			    generic_version<compressed2D<Value, Para>, D, true>
+	                   >::type {};
+
+	/// Solve \p w = A * \p v
+	template <typename VectorIn, typename VectorOut>
+	void operator()(const VectorIn& v, VectorOut& w) const
+	{   vampir_trace<5022> tracer; apply(v, w, version<Matrix, DiaTag, CompactStorage>()); }
+
+
+      private:
+	template <typename Value>
+	Value inline lower_trisolve_diavalue(const Value& v, tag::regular_diagonal) const
+	{   using math::reciprocal; return reciprocal(v); }
+
+	template <typename Value>
+	Value lower_trisolve_diavalue(const Value& v, tag::inverse_diagonal) const
+	{  return v;	}
+
+	template <typename Tag> int dia_inc(Tag) { return 0; }
+	int dia_inc(tag::unit_diagonal) { return 1; }
+
+	// Generic row-major unit_diagonal
+	template <typename VectorIn, typename VectorOut>
+	void apply(const VectorIn& v, VectorOut& w, boost::mpl::int_<1>) const
+	{
+	    using namespace tag;
+	    a_cur_type ac= begin<row>(A), aend= end<row>(A);
+	    for (size_type r= 0; ac != aend; ++r, ++ac) {
+		a_icur_type aic= begin<nz>(ac), aiend= CompactStorage ? end<nz>(ac) : lower_bound<nz>(ac, r);
+		typename Collection<VectorOut>::value_type rr= v[r];
+		for (; aic != aiend; ++aic) {
+		    MTL_DEBUG_THROW_IF(col_a(*aic) >= r, logic_error("Matrix entries must be sorted for this."));
+		    rr-= value_a(*aic) * w[col_a(*aic)];
+		}
+		w[r]= rr;
+	    }
+	}
+
+	// Generic row-major not unit_diagonal
+	template <typename VectorIn, typename VectorOut>
+	void apply(const VectorIn& v, VectorOut& w, boost::mpl::int_<2>) const
+	{
+	    using namespace tag;
+	    a_cur_type ac= begin<row>(A), aend= end<row>(A);
+	    for (size_type r= 0; ac != aend; ++r, ++ac) {
+		a_icur_type aic= begin<nz>(ac), aiend= CompactStorage ? end<nz>(ac) : lower_bound<nz>(ac, r+1);
+		MTL_THROW_IF(aic == aiend, missing_diagonal());
+		--aiend;
+		MTL_THROW_IF(col_a(*aiend) != r, missing_diagonal());
+
+		value_type dia= value_a(*aiend);
+		typename Collection<VectorOut>::value_type rr= v[r];
+
+		for (; aic != aiend; ++aic) {
+		    MTL_DEBUG_THROW_IF(col_a(*aic) >= r, logic_error("Matrix entries must be sorted for this."));
+		    rr-= value_a(*aic) * w[col_a(*aic)];
+		}
+		w[r]= rr * lower_trisolve_diavalue(dia, DiaTag());
+	    }
+	}
+
+	// Generic column-major unit_diagonal
+	template <typename VectorIn, typename VectorOut>
+	void apply(const VectorIn& v, VectorOut& w, boost::mpl::int_<3>) const
+	{
+	    using namespace tag;
+	    w= v;
+	    a_cur_type ac= begin<col>(A), aend= end<col>(A);
+	    for (size_type r= 0; ac != aend; ++r, ++ac) {
+		a_icur_type aic= CompactStorage ? begin<nz>(ac) : lower_bound<nz>(ac, r+1), aiend= end<nz>(ac);
+		typename Collection<VectorOut>::value_type rr= w[r];
+
+		for (; aic != aiend; ++aic) {
+		    MTL_DEBUG_THROW_IF(row_a(*aic) <= r, logic_error("Matrix entries must be sorted for this."));
+		    w[row_a(*aic)]-= value_a(*aic) * rr;
+		}
+	    }
+	}
+
+	// Generic column-major not unit_diagonal
+	template <typename VectorIn, typename VectorOut>
+	void apply(const VectorIn& v, VectorOut& w, boost::mpl::int_<4>) const
+	{
+	    using namespace tag;
+	    w= v;
+	    a_cur_type ac= begin<col>(A), aend= end<col>(A);
+	    for (size_type r= 0; ac != aend; ++r, ++ac) {
+		a_icur_type aic= CompactStorage ? begin<nz>(ac) : lower_bound<nz>(ac, r), aiend= end<nz>(ac);
+		MTL_DEBUG_THROW_IF(aic == aiend || row_a(*aic) != r, missing_diagonal());
+		typename Collection<VectorOut>::value_type rr= w[r]*= lower_trisolve_diavalue(value_a(*aic), DiaTag());
+
+		for (++aic; aic != aiend; ++aic) {
+		    MTL_DEBUG_THROW_IF(row_a(*aic) <= r, logic_error("Matrix entries must be sorted for this."));
+		    w[row_a(*aic)]-= value_a(*aic) * rr;
+		}
+	    }
+	}
+
+	// Tuning for IC_0 and similar using compressed2D row-major compact with implicit unit diagonal
+	template <typename VectorIn, typename VectorOut>
+	void apply(const VectorIn& v, VectorOut& w, boost::mpl::int_<5>) const
+	{
+	    vampir_trace<5048> tracer;
+	    if (num_rows(A) == 0) return;
+	    size_type j1= A.ref_major()[1];
+	    for (size_type r= 0, rend= num_rows(A); r != rend; ++r) {
+		size_type j0= j1;
+		j1= A.ref_major()[r+1];
+		typename Collection<VectorOut>::value_type rr= v[r];
+		for (; j0 != j1; ++j0) {
+		    MTL_DEBUG_THROW_IF(A.ref_minor()[j0] > r, logic_error("Matrix entries from U in lower triangular."));
+		    rr-= A.data[j0] * w[A.ref_minor()[j0]];
+		}
+		w[r]= rr;
+	    }
+	}
+
+	// Tuning for IC_0 and similar using compressed2D row-major compact with explicitly stored diagonal (possibly already inverted)
+	template <typename VectorIn, typename VectorOut>
+	void apply(const VectorIn& v, VectorOut& w, boost::mpl::int_<6>) const
+	{
+	    vampir_trace<5047> tracer;
+	    for (size_type r= 0, rend= num_rows(A); r != rend; ++r) {
+		size_type j0= A.ref_major()[r], j1= A.ref_major()[r+1];
+		MTL_THROW_IF(j0 == j1, missing_diagonal());
+		--j1;
+		MTL_THROW_IF(A.ref_minor()[j1] != r, missing_diagonal());
+		value_type dia= A.data[j1];
+		typename Collection<VectorOut>::value_type rr= v[r];
+		for (; j0 != j1; ++j0) {
+		    MTL_DEBUG_THROW_IF(A.ref_minor()[j0] > r, logic_error("Matrix entries from U in lower triangular."));
+		    rr-= A.data[j0] * w[A.ref_minor()[j0]];
+		}
+		w[r]= rr * lower_trisolve_diavalue(dia, DiaTag());
+	    }
+	}
+
+	const Matrix&                                    A;
+	typename mtl::traits::const_value<Matrix>::type  value_a;
+	typename mtl::traits::col<Matrix>::type          col_a;
+	typename mtl::traits::row<Matrix>::type          row_a;
+    };
+
+}  // detail
+
+
+/// Solves the lower triangular matrix A  with the rhs v and returns the solution vector
+template <typename Matrix, typename Vector>
+Vector inline lower_trisolve(const Matrix& A, const Vector& v)
+{
+    Vector w(resource(v));
+    detail::lower_trisolve_t<Matrix, tag::regular_diagonal> solver(A);
+    solver(v, w);
+    return w;
+}
+
+/// Solves the lower triangular matrix A  with the rhs v
+template <typename Matrix, typename VectorIn, typename VectorOut>
+inline void lower_trisolve(const Matrix& A, const VectorIn& v, VectorOut& w)
+{
+    detail::lower_trisolve_t<Matrix, tag::regular_diagonal> solver(A);
+    solver(v, w);
+}
+
+/// Solves the lower triangular matrix A (only one's in the diagonal) with the rhs v and returns the solution vector
+template <typename Matrix, typename Vector>
+Vector inline unit_lower_trisolve(const Matrix& A, const Vector& v)
+{
+    Vector w(resource(v));
+    detail::lower_trisolve_t<Matrix, tag::unit_diagonal> solver(A);
+    solver(v, w);
+    return w;
+}
+
+/// Solves the lower triangular matrix A (only one's in the diagonal) with the rhs v and returns the solution vector
+template <typename Matrix, typename VectorIn, typename VectorOut>
+inline void unit_lower_trisolve(const Matrix& A, const VectorIn& v, VectorOut& w)
+{
+    detail::lower_trisolve_t<Matrix, tag::unit_diagonal> solver(A);
+    solver(v, w);
+}
+
+/// Solves the lower triangular matrix A (inverse the diagonal) with the rhs v and returns the solution vector
+template <typename Matrix, typename Vector>
+Vector inline inverse_lower_trisolve(const Matrix& A, const Vector& v)
+{
+    Vector w(resource(v));
+    detail::lower_trisolve_t<Matrix, tag::inverse_diagonal> solver(A);
+    solver(v, w);
+    return w;
+}
+
+/// Solves the lower triangular matrix A (inverse the diagonal) with the rhs v and returns the solution vector
+template <typename Matrix, typename VectorIn, typename VectorOut>
+inline void inverse_lower_trisolve(const Matrix& A, const VectorIn& v, VectorOut& w)
+{
+    detail::lower_trisolve_t<Matrix, tag::inverse_diagonal> solver(A);
+    solver(v, w);
+}
+
+template <typename Matrix, typename Vector, typename DiaTag>
+Vector inline lower_trisolve(const Matrix& A, const Vector& v, DiaTag)
+{
+    Vector w(resource(v));
+    detail::lower_trisolve_t<Matrix, DiaTag> solver(A);
+    solver(v, w);
+    return w;
+}
+
+template <typename Matrix, typename VectorIn, typename VectorOut, typename DiaTag>
+inline void lower_trisolve(const Matrix& A, const VectorIn& v, VectorOut& w, DiaTag)
+{
+    detail::lower_trisolve_t<Matrix, DiaTag> solver(A);
+    solver(v, w);
+}
+
+}} // namespace mtl::matrix
+
+#endif // MTL_LOWER_TRISOLVE_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/lu.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_LU_INCLUDE
+#define MTL_MATRIX_LU_INCLUDE
+
+#include <cmath>
+#include <boost/numeric/linear_algebra/identity.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/mtl/utility/irange.hpp>
+#include <boost/numeric/mtl/utility/lu_matrix_type.hpp>
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/matrix/upper.hpp>
+#include <boost/numeric/mtl/matrix/lower.hpp>
+#include <boost/numeric/mtl/matrix/permutation.hpp>
+#include <boost/numeric/mtl/operation/adjoint.hpp>
+#include <boost/numeric/mtl/operation/lower_trisolve.hpp>
+#include <boost/numeric/mtl/operation/upper_trisolve.hpp>
+#include <boost/numeric/mtl/operation/max_pos.hpp>
+#include <boost/numeric/mtl/operation/swap_row.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+
+namespace mtl { namespace matrix {
+
+/// LU factorization in place (without pivoting and optimization so far)
+/** eps is tolerance for pivot element. If less or equal the matrix is considered singular.
+    eps is given as double right now, might be refactored to the magnitude type of the value type in the future. **/
+template <typename Matrix>
+void inline lu(Matrix& LU, double eps= 0)
+{
+    vampir_trace<5023> tracer;
+    using std::abs;
+    MTL_THROW_IF(num_rows(LU) != num_cols(LU), matrix_not_square());
+
+    for (std::size_t k= 0; k < num_rows(LU)-1; k++) {
+	if(abs(LU[k][k]) <= eps) throw matrix_singular();
+	irange r(k+1, imax); // Interval [k+1, n-1]
+	LU[r][k]/= LU[k][k];
+	LU[r][r]-= LU[r][k] * LU[k][r];
+    }
+}
+
+/// LU factorization in place (with pivoting and without optimization so far)
+/** eps is tolerance for pivot element. If less or equal the matrix is considered singular.
+    eps is given as double right now, might be refactored to the magnitude type of the value type in the future. **/
+template <typename Matrix, typename PermuationVector>
+void inline lu(Matrix& A, PermuationVector& P, double eps= 0)
+{
+    vampir_trace<5024> tracer;
+    using math::zero; using std::abs;
+    typedef typename Collection<Matrix>::size_type    size_type;
+    size_type ncols = num_cols(A), nrows = num_rows(A);
+
+    MTL_THROW_IF(ncols != nrows , matrix_not_square());
+    P.change_dim(nrows);
+
+    for (size_type i= 0; i < nrows; i++)
+        P[i]= i;
+
+    for (size_type i= 0; i < nrows-1; i++) {
+	irange r(i+1, imax), ir(i, i+1); // Intervals [i+1, n-1], [i, i]
+	size_type rmax= max_abs_pos(A[irange(i, imax)][ir]).first + i;
+	swap_row(A, i, rmax);
+	swap_row(P, i, rmax);
+
+	if(abs(A[i][i]) <= eps) throw matrix_singular(); // other gmres test doesn't work
+
+	A[r][i]/= A[i][i];              // Scale column i
+	A[r][r]-= A[r][i] * A[i][r]; 	 // Decrease bottom right block of matrix
+    }
+}
+
+
+/// LU factorization without factorization that returns the matrix
+/** eps is tolerance for pivot element. If less or equal the matrix is considered singular.
+    eps is given as double right now, might be refactored to the magnitude type of the value type in the future. **/
+template <typename Matrix>
+Matrix inline lu_f(const Matrix& A, double eps= 0)
+{
+    vampir_trace<5025> tracer;
+    Matrix LU(A);
+    lu(LU, eps);
+    return LU;
+}
+
+/// Solve Ax = b by LU factorization without pivoting; vector x is returned
+template <typename Matrix, typename Vector>
+Vector inline lu_solve_straight(const Matrix& A, const Vector& b, double eps= 0)
+{
+    vampir_trace<5026> tracer;
+    Matrix LU(A);
+    lu(LU, eps);
+    return upper_trisolve(upper(LU), unit_lower_trisolve(strict_lower(LU), b));
+}
+
+/// Solve LUx = b by with forward and backward-LU subsitution, lu(LU) was allready done
+template <typename Matrix, typename Vector>
+Vector inline lu_solve_apply(const Matrix& LU, const Vector& b)
+{
+    vampir_trace<5026> tracer;
+    return upper_trisolve(upper(LU), unit_lower_trisolve(strict_lower(LU), b));
+}
+
+/// Apply the factorization L*U with permutation P on vector b to solve Ax = b
+template <typename Matrix, typename PermVector, typename Vector>
+Vector inline lu_apply(const Matrix& LU, const PermVector& P, const Vector& b)
+{
+    vampir_trace<5027> tracer;
+    return upper_trisolve(upper(LU), unit_lower_trisolve(strict_lower(LU), Vector(matrix::permutation(P) * b)));
+}
+
+
+/// Solve Ax = b by LU factorization with column pivoting; vector x is returned
+template <typename Matrix, typename Vector>
+Vector inline lu_solve(const Matrix& A, const Vector& b, double eps= 0)
+{
+    vampir_trace<5028> tracer;
+    mtl::vector::dense_vector<std::size_t, vector::parameters<> > P(num_rows(A));
+    Matrix                    LU(A);
+
+    lu(LU, P, eps);
+    return lu_apply(LU, P, b);
+}
+
+
+/// Apply the factorization L*U with permutation P on vector b to solve adjoint(A)x = b
+/** That is \f$P^{-1}(LU)^H x = b\f$ --> \f$x= P^{-1}L^{-H} U^{-H} b\f$ where \f$P^{{-1}^{{-1}^H}} = P^{-1}\f$ **/
+template <typename Matrix, typename PermVector, typename Vector>
+Vector inline lu_adjoint_apply(const Matrix& LU, const PermVector& P, const Vector& b)
+{
+    vampir_trace<5029> tracer;
+    return Vector(trans(matrix::permutation(P)) * unit_upper_trisolve(adjoint(LU), lower_trisolve(adjoint(LU), b)));
+}
+
+
+/// Solve \f$adjoint(A)x = b\f$ by LU factorization with column pivoting; vector x is returned
+template <typename Matrix, typename Vector>
+Vector inline lu_adjoint_solve(const Matrix& A, const Vector& b, double eps= 0)
+{
+    vampir_trace<5030> tracer;
+    mtl::vector::dense_vector<std::size_t, vector::parameters<> > P(num_rows(A));
+    Matrix                    LU(A);
+
+    lu(LU, P, eps);
+    return lu_adjoint_apply(LU, P, b);
+}
+
+/// Class that keeps LU factorization (and permutation); using column pivoting
+template <typename Matrix>
+class lu_solver
+{
+    typedef typename mtl::traits::lu_matrix_type<Matrix>::type            matrix_type;
+    typedef mtl::vector::dense_vector<std::size_t, vector::parameters<> > permutation_type;
+  public:
+    /// Construct from matrix \p A and use optionally threshold \p eps in factorization
+    explicit lu_solver(const Matrix& A, double eps= 0)
+      : LU(A), P(num_rows(A))
+    {
+	lu(LU, P, eps);
+    }
+
+    /// Solve A*x = b with factorization from constructor
+    template <typename VectorIn, typename VectorOut>
+    void solve(const VectorIn& b, VectorOut& x) const
+    {
+	x= upper_trisolve(upper(LU), unit_lower_trisolve(strict_lower(LU), VectorIn(matrix::permutation(P) * b)));
+    }
+    /// Solve \f$adjoint(A)x = b\f$ using LU factorization
+    template <typename VectorIn, typename VectorOut>
+    void adjoint_solve(const VectorIn& b, VectorOut& x) const
+    {
+	x= trans(matrix::permutation(P)) * unit_upper_trisolve(adjoint(LU), lower_trisolve(adjoint(LU), b));
+    }
+
+  private:
+    matrix_type      LU;
+    permutation_type P;
+};
+
+
+}} // namespace mtl::matrix
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+// ### For illustration purposes
+#if 0
+
+namespace mtl { namespace matrix {
+
+template <typename Matrix>
+void inline lu(Matrix& LU)
+{
+    MTL_THROW_IF(num_rows(LU) != num_cols(LU), matrix_not_square());
+
+    typedef typename Collection<Matrix>::value_type   value_type;
+    typedef typename Collection<Matrix>::size_type    size_type;
+
+    size_type n= num_rows(LU);
+    for (size_type k= 0; k < num_rows(LU); k++) {
+	value_type pivot= LU[k][k];
+	for (size_type j= k+1; j < n; j++) {
+	    value_type alpha= LU[j][k]/= pivot;
+	    for (size_type i= k+1; i < n; i++)
+		LU[j][i]-= alpha * LU[k][i];
+	}
+    }
+}
+
+
+
+}} // namespace mtl::matrix
+
+#endif
+
+#endif // MTL_MATRIX_LU_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/make_sparse.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, www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also tools/license/license.mtl.txt in the distribution.
+
+#ifndef MTL_VECTOR_MAKE_SPARSE_INCLUDE
+#define MTL_VECTOR_MAKE_SPARSE_INCLUDE
+
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/operation/size.hpp>
+#include <boost/numeric/mtl/operation/update.hpp>
+#include <boost/numeric/mtl/matrix/compressed2D.hpp>
+#include <boost/numeric/mtl/matrix/inserter.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+
+namespace mtl { namespace vector {
+
+// Commands in Matlab
+// S = sparse(i,j,s,m,n,nzmax)
+// S = sparse(i,j,s,m,n)
+// S = sparse(i,j,s)
+// S = sparse(m,n)
+
+template <typename SizeVector, typename ValueVector>
+struct make_sparse_trait
+{
+    typedef typename Collection<SizeVector>::value_type  size_type;
+    typedef typename Collection<ValueVector>::value_type value_type;
+    typedef matrix::parameters<row_major, index::c_index, mtl::non_fixed::dimensions, false, size_type> paras;
+    typedef matrix::compressed2D<value_type, paras>   type;
+};
+
+/// Generates an \p m by \p n matrix from the vectors \p rows, \p cols, and \p values.
+/** A sparse matrix is created (compressed2D). The value type is the same as the element type
+    of the value vector and the size type the same as the entries of the vectors with the row indices.
+    Zero entries in \p values are ignored. Entries with same coordinates are added.
+    Same as <a href="http://www.mathworks.de/de/help/matlab/ref/sparse.html">Matlab's sparse</a> function besides that it is zero-indexed.
+ **/
+template <typename SizeVector1, typename SizeVector2, typename ValueVector>
+inline typename make_sparse_trait<SizeVector1, ValueVector>::type
+make_sparse(const SizeVector1& rows, const SizeVector2& cols, const ValueVector& values,
+	    std::size_t m, std::size_t n)
+{
+    MTL_THROW_IF(size(rows) != size(cols), incompatible_size());
+    MTL_THROW_IF(size(rows) != size(values), incompatible_size());
+
+    typedef make_sparse_trait<SizeVector1, ValueVector> traits;
+    typedef typename traits::type       matrix_type;
+    typedef typename traits::value_type value_type;
+    typedef typename traits::size_type  size_type;
+
+    size_type               ms(m), ns(n);  // shouldn't be needed :-!
+    matrix_type             A(ms, ns);
+    matrix::inserter<matrix_type, update_plus<value_type> > ins(A, size_type(size(rows) / m + 1));
+
+    for (std::size_t i= 0; i < size(rows); i++)
+	if (values[i] != value_type(0))
+	    ins[rows[i]][cols[i]] << values[i];
+    return A;
+}
+
+/// Generates an \p m by \p n matrix from the vectors \p rows, \p cols, and \p values.
+/** A sparse matrix is created (compressed2D). The value type is the same as the element type
+    of the value vector and the size type the same as the entries of the vectors with the row indices.
+    Zero entries in \p values are ignored. Entries with same coordinates are added. Last parameter
+    is ignored and can be omitted.
+    Same as <a href="http://www.mathworks.de/de/help/matlab/ref/sparse.html">Matlab's sparse</a> function besides that it is zero-indexed.
+ **/
+template <typename SizeVector1, typename SizeVector2, typename ValueVector>
+inline typename make_sparse_trait<SizeVector1, ValueVector>::type
+make_sparse(const SizeVector1& rows, const SizeVector2& cols, const ValueVector& values,
+	    std::size_t m, std::size_t n, std::size_t)
+{
+    return make_sparse(rows, cols, values, m, n);
+}
+
+/// Generates a matrix from the vectors \p rows, \p cols, and \p values.
+/** A sparse matrix is created (compressed2D).
+    The number of rows/columns is one plus the maximum of the entries of \p rows and \p cols.
+    The value type is the same as the element type
+    of the value vector and the size type the same as the entries of the vectors with the row indices.
+    Zero entries in \p values are ignored. Entries with same coordinates are added.
+    Same as <a href="http://www.mathworks.de/de/help/matlab/ref/sparse.html">Matlab's sparse</a> function besides that it is zero-indexed.
+ **/
+template <typename SizeVector1, typename SizeVector2, typename ValueVector>
+inline typename make_sparse_trait<SizeVector1, ValueVector>::type
+make_sparse(const SizeVector1& rows, const SizeVector2& cols, const ValueVector& values)
+{
+    return make_sparse(rows, cols, values, max(rows)+1, max(cols)+1);
+}
+
+/// Generates an empty \p m by \p n matrix.
+/** A sparse matrix is created (compressed2D<double>).
+    Same as <a href="http://www.mathworks.de/de/help/matlab/ref/sparse.html">Matlab's sparse</a> function besides that it is zero-indexed.
+ **/
+inline matrix::compressed2D<double> make_sparse(std::size_t m, std::size_t n)
+{
+    return matrix::compressed2D<double>(m, n);
+}
+
+
+
+
+
+} // namespace :vector
+
+    using vector::make_sparse;
+
+} // namespace mtl
+
+#endif // MTL_VECTOR_MAKE_SPARSE_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/make_tag_vector.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, www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also tools/license/license.mtl.txt in the distribution.
+
+#ifndef MTL_MAKE_TAG_VECTOR_INCLUDE
+#define MTL_MAKE_TAG_VECTOR_INCLUDE
+
+#include <boost/numeric/mtl/utility/irange.hpp>
+#include <boost/numeric/mtl/utility/srange.hpp>
+#include <boost/numeric/mtl/vector/dense_vector.hpp>
+
+namespace mtl {
+
+    template <typename Range>
+    inline vector::dense_vector<bool>
+    make_tag_vector(std::size_t n, const Range& r)
+    {
+	vector::dense_vector<bool> v(n);
+	for (std::size_t i= 0; i < n; i++)
+	    v[i]= r.in_range(i);
+	return v;
+    }
+
+} // namespace mtl
+
+#endif // MTL_MAKE_TAG_VECTOR_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/mat_cvec_times_expr.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_MAT_CVEC_TIMES_EXPR_INCLUDE
+#define MTL_MAT_CVEC_TIMES_EXPR_INCLUDE
+
+#include <boost/numeric/mtl/operation/bin_op_expr.hpp>
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/concept/std_concept.hpp>
+#include <boost/numeric/mtl/mtl_fwd.hpp>
+
+namespace mtl {
+
+template <typename E1, typename E2>
+struct mat_cvec_times_expr
+  : public bin_op_expr< E1, E2 >,
+    public mtl::vector::vec_expr< mat_cvec_times_expr<E1, E2> >
+{
+    typedef bin_op_expr< E1, E2 >         base;
+    typedef mat_cvec_times_expr<E1, E2>   self;
+
+    typedef typename Multiplicable<typename Collection<E1>::value_type,
+				   typename Collection<E2>::value_type>::result_type   value_type;
+    typedef typename Collection<E2>::size_type    size_type;
+
+    mat_cvec_times_expr( E1 const& matrix, E2 const& vector ) : base(matrix, vector)  {}
+};
+
+    namespace vector {
+	template <typename E1, typename E2>
+	std::size_t inline size(const mat_cvec_times_expr<E1, E2>& x)
+	{ return num_rows(x.first); }
+    }
+
+} // namespace mtl
+
+#endif // MTL_MAT_CVEC_TIMES_EXPR_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/mat_vec_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_MAT_VEC_MULT_INCLUDE
+#define MTL_MAT_VEC_MULT_INCLUDE
+
+#include <cassert>
+// #include <iostream>
+#include <boost/mpl/bool.hpp>
+
+#include <boost/numeric/mtl/config.hpp>
+#include <boost/numeric/mtl/mtl_fwd.hpp>
+#include <boost/numeric/mtl/utility/category.hpp>
+#include <boost/numeric/mtl/utility/property_map.hpp>
+#include <boost/numeric/mtl/utility/range_generator.hpp>
+#include <boost/numeric/mtl/utility/tag.hpp>
+#include <boost/numeric/mtl/utility/is_static.hpp>
+#include <boost/numeric/mtl/utility/tag.hpp>
+#include <boost/numeric/mtl/utility/enable_if.hpp>
+#include <boost/numeric/mtl/utility/multi_tmp.hpp>
+#include <boost/numeric/mtl/utility/static_assert.hpp>
+#include <boost/numeric/mtl/vector/dense_vector.hpp>
+#include <boost/numeric/mtl/utility/omp_size_type.hpp>
+#include <boost/numeric/mtl/operation/set_to_zero.hpp>
+#include <boost/numeric/mtl/operation/update.hpp>
+#include <boost/numeric/linear_algebra/identity.hpp>
+#include <boost/numeric/meta_math/loop.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+
+namespace mtl { namespace matrix {
+
+namespace impl {
+
+    template <std::size_t Index0, std::size_t Max0, std::size_t Index1, std::size_t Max1, typename Assign>
+    struct fully_unroll_mat_cvec_mult
+      : public meta_math::loop2<Index0, Max0, Index1, Max1>
+    {
+	typedef meta_math::loop2<Index0, Max0, Index1, Max1>                              base;
+	typedef fully_unroll_mat_cvec_mult<base::next_index0, Max0, base::next_index1, Max1, Assign>  next_t;
+
+	template <typename Matrix, typename VectorIn, typename VectorOut>
+	static inline void apply(const Matrix& A, const VectorIn& v, VectorOut& w)
+	{
+	    Assign::update(w[base::index0], A[base::index0][base::index1] * v[base::index1]);
+	    next_t::apply(A, v, w);
+	}
+    };
+
+    // need specialization here for not going back to column 0 but column 1
+    template <std::size_t Index0, std::size_t Max0, std::size_t Max1, typename Assign>
+    struct fully_unroll_mat_cvec_mult<Index0, Max0, Max1, Max1, Assign>
+      : public meta_math::loop2<Index0, Max0, Max1, Max1>
+    {
+	typedef meta_math::loop2<Index0, Max0, Max1, Max1>                              base;
+	typedef fully_unroll_mat_cvec_mult<base::next_index0, Max0, 2, Max1, Assign>  next_t;
+
+	template <typename Matrix, typename VectorIn, typename VectorOut>
+	static inline void apply(const Matrix& A, const VectorIn& v, VectorOut& w)
+	{
+	    Assign::update(w[base::index0], A[base::index0][base::index1] * v[base::index1]);
+	    next_t::apply(A, v, w);
+	}
+    };
+
+    template <std::size_t Max0, std::size_t Max1, typename Assign>
+    struct fully_unroll_mat_cvec_mult<Max0, Max0, Max1, Max1, Assign>
+      : public meta_math::loop2<Max0, Max0, Max1, Max1>
+    {
+	typedef meta_math::loop2<Max0, Max0, Max1, Max1>                              base;
+
+	template <typename Matrix, typename VectorIn, typename VectorOut>
+	static inline void apply(const Matrix& A, const VectorIn& v, VectorOut& w)
+	{
+	    Assign::update(w[base::index0], A[base::index0][base::index1] * v[base::index1]);
+	}
+    };
+
+    struct noop
+    {
+	template <typename Matrix, typename VectorIn, typename VectorOut>
+	static inline void apply(const Matrix&, const VectorIn&, VectorOut&) {}
+    };
+} // impl
+
+// Dense matrix vector multiplication with run-time matrix size
+template <typename Matrix, typename VectorIn, typename VectorOut, typename Assign>
+inline void dense_mat_cvec_mult(const Matrix& A, const VectorIn& v, VectorOut& w, Assign, boost::mpl::true_)
+{
+    vampir_trace<3017> tracer;
+    typedef typename static_num_rows<Matrix>::type size_type;
+    static const size_type rows_a= static_num_rows<Matrix>::value, cols_a= static_num_cols<Matrix>::value;
+
+    assert(rows_a > 0 && cols_a > 0);
+    // w= A[all][0] * v[0];  N.B.: 1D is unrolled by the compiler faster (at least on gcc)
+    for (size_type i= 0; i < rows_a; i++)
+	Assign::first_update(w[i], A[i][0] * v[0]);
+
+    // corresponds to w+= A[all][1:] * v[1:]; if necessary
+    typedef impl::fully_unroll_mat_cvec_mult<1, rows_a, 2, cols_a, Assign>  f2;
+    typedef typename boost::mpl::if_c<(cols_a > 1), f2, impl::noop>::type   f3;
+    f3::apply(A, v, w);
+}
+
+template <unsigned Index, unsigned Size>
+struct init_ptrs
+{
+    template <typename Matrix, typename Ptrs>
+    inline static void apply(const Matrix& A, Ptrs& ptrs)
+    {
+	ptrs.value= &A[Index][0];
+	init_ptrs<Index+1, Size>::apply(A, ptrs.sub);
+    }
+};
+
+template <unsigned Size>
+struct init_ptrs<Size, Size>
+{
+    template <typename Matrix, typename Ptrs>
+    inline static void apply(const Matrix&, Ptrs&) {}
+};
+
+template <unsigned Index, unsigned Size>
+struct square_cvec_mult_rows
+{
+    template <typename Tmps, typename Ptrs, typename ValueIn>
+    inline static void compute(Tmps& tmps, Ptrs& ptrs, ValueIn vi)
+    {
+	tmps.value+= *ptrs.value++ * vi;
+	square_cvec_mult_rows<Index+1, Size>::compute(tmps.sub, ptrs.sub, vi);
+    }
+
+    template <typename As, typename VectorOut, typename Tmps>
+    inline static void update(As, VectorOut& w, const Tmps& tmps)
+    {
+	As::first_update(w[Index], tmps.value);
+	square_cvec_mult_rows<Index+1, Size>::update(As(), w, tmps.sub);
+    }
+};
+
+template <unsigned Size>
+struct square_cvec_mult_rows<Size, Size>
+{
+    template <typename Tmps, typename Ptrs, typename ValueIn>
+    inline static void compute(Tmps&, Ptrs&, ValueIn) {}
+
+    template <typename As, typename VectorOut, typename Tmps>
+    inline static void update(As, VectorOut&, const Tmps&) {}
+};
+
+template <unsigned Index, unsigned Size>
+struct square_cvec_mult_cols
+{
+    template <typename Tmps, typename Ptrs, typename VPtr>
+    inline static void compute(Tmps& tmps, Ptrs& ptrs, VPtr vp)
+    {
+	square_cvec_mult_rows<0, Size>::compute(tmps, ptrs, *vp);
+	square_cvec_mult_cols<Index+1, Size>::compute(tmps, ptrs, ++vp);
+    }
+};
+
+template <unsigned Size>
+struct square_cvec_mult_cols<Size, Size>
+{
+    template <typename Tmps, typename Ptrs, typename VPtr>
+    inline static void compute(Tmps&, Ptrs&, VPtr) {}
+};
+
+
+// Dense matrix vector multiplication with run-time matrix size
+template <unsigned Size, typename MValue, typename MPara, typename ValueIn, typename ParaIn,
+	  typename VectorOut, typename Assign>
+inline void square_cvec_mult(const dense2D<MValue, MPara>& A, const mtl::vector::dense_vector<ValueIn, ParaIn>& v, VectorOut& w, Assign)
+{
+    // vampir_trace<3067> tracer;
+    MTL_STATIC_ASSERT((mtl::traits::is_row_major<MPara>::value), "Only row-major matrices supported in this function.");
+
+    typedef typename Collection<VectorOut>::value_type value_type;
+    multi_tmp<Size, value_type> tmps(math::zero(w[0]));
+
+    multi_tmp<Size, const MValue*>    ptrs;
+    init_ptrs<0, Size>::apply(A, ptrs);
+
+    const ValueIn* vp= &v[0];
+
+    square_cvec_mult_cols<0, Size>::compute(tmps, ptrs, vp); // outer loop over columns
+    square_cvec_mult_rows<0, Size>::update(Assign(), w, tmps);         // update rows
+}
+
+
+// Dense matrix vector multiplication with run-time matrix size
+template <typename MValue, typename MPara, typename ValueIn, typename ParaIn, typename VectorOut, typename Assign>
+typename boost::enable_if<mtl::traits::is_row_major<MPara> >::type
+inline dense_mat_cvec_mult(const dense2D<MValue, MPara>& A, const mtl::vector::dense_vector<ValueIn, ParaIn>& v, VectorOut& w, Assign, boost::mpl::false_)
+{
+    // vampir_trace<3066> tracer;
+
+    using math::zero;
+    if (mtl::vector::size(w) == 0) return;
+
+    typedef typename Collection<VectorOut>::value_type value_type;
+    // typedef ValueIn                                    value_in_type;
+    typedef typename MPara::size_type                  size_type;
+
+    const size_type  nr= num_rows(A), nc= num_cols(A);
+    if (nr == nc && nr <= 8)
+	switch (nr) {
+	  case 1: Assign::first_update(w[0], A[0][0] * v[0]); return;
+	  case 2: square_cvec_mult<2>(A, v, w, Assign()); return;
+	  case 3: square_cvec_mult<3>(A, v, w, Assign()); return;
+	  case 4: square_cvec_mult<4>(A, v, w, Assign()); return;
+	  case 5: square_cvec_mult<5>(A, v, w, Assign()); return;
+	  case 6: square_cvec_mult<6>(A, v, w, Assign()); return;
+	  case 7: square_cvec_mult<7>(A, v, w, Assign()); return;
+	  case 8: square_cvec_mult<8>(A, v, w, Assign()); return;
+	}
+
+
+    const size_type  nrb= nr / 4 * 4;
+    const value_type z(math::zero(w[0]));
+
+    for (size_type i= 0; i < nrb; i+= 4) {
+	value_type      tmp0(z), tmp1(z), tmp2(z), tmp3(z);
+	const MValue *p0= &A[i][0], *pe= p0 + nc, *p1= &A[i+1][0], *p2= &A[i+2][0], *p3= &A[i+3][0];
+	const ValueIn* vp= &v[0];
+	for (; p0 != pe; ) {
+	    const ValueIn vj= *vp++;
+	    tmp0+= *p0++ * vj;
+	    tmp1+= *p1++ * vj;
+	    tmp2+= *p2++ * vj;
+	    tmp3+= *p3++ * vj;
+	}
+	Assign::first_update(w[i], tmp0);
+	Assign::first_update(w[i+1], tmp1);
+	Assign::first_update(w[i+2], tmp2);
+	Assign::first_update(w[i+3], tmp3);
+    }
+
+    for (size_type i= nrb; i < nr; i++) {
+	value_type tmp= z;
+	const ValueIn* vp= &v[0];
+	for (const MValue *p0= &A[i][0], *pe= p0 + nc; p0 != pe; )
+	    tmp+= *p0++ * *vp++;
+	Assign::first_update(w[i], tmp);
+    }
+}
+
+
+// Dense matrix vector multiplication with run-time matrix size
+template <typename Matrix, typename VectorIn, typename VectorOut, typename Assign>
+inline void dense_mat_cvec_mult(const Matrix& A, const VectorIn& v, VectorOut& w, Assign, boost::mpl::false_)
+{
+    vampir_trace<3018> tracer;
+    // Naive implementation, will be moved to a functor and complemented with more efficient ones
+
+    using math::zero;
+    if (mtl::vector::size(w) == 0) return;
+    // std::cout << "Bin in richtiger Funktion\n";
+
+    // if (Assign::init_to_zero) set_to_zero(w); // replace update with first_update insteda
+
+    typedef typename Collection<VectorOut>::value_type value_type;
+    typedef typename Collection<VectorIn>::value_type  value_in_type;
+    typedef typename Collection<Matrix>::size_type     size_type;
+
+    const value_type z(math::zero(w[0]));
+    const size_type nr= num_rows(A), nrb= nr / 4 * 4, nc= num_cols(A);
+
+    for (size_type i= 0; i < nrb; i+= 4) {
+	value_type      tmp0(z), tmp1(z), tmp2(z), tmp3(z);
+	for (size_type j= 0; j < nc; j++) {
+	    const value_in_type vj= v[j];
+	    tmp0+= A[i][j] * vj;
+	    tmp1+= A[i+1][j] * vj;
+	    tmp2+= A[i+2][j] * vj;
+	    tmp3+= A[i+3][j] * vj;
+	}
+	Assign::first_update(w[i], tmp0);
+	Assign::first_update(w[i+1], tmp1);
+	Assign::first_update(w[i+2], tmp2);
+	Assign::first_update(w[i+3], tmp3);
+    }
+
+    for (size_type i= nrb; i < nr; i++) {
+	value_type tmp= z;
+	for (size_type j= 0; j < nc; j++)
+	    tmp+= A[i][j] * v[j];
+	Assign::first_update(w[i], tmp);
+    }
+}
+
+// Dense matrix vector multiplication
+template <typename Matrix, typename VectorIn, typename VectorOut, typename Assign>
+inline void mat_cvec_mult(const Matrix& A, const VectorIn& v, VectorOut& w, Assign, tag::flat<tag::dense>)
+{
+# ifdef MTL_NOT_UNROLL_FSIZE_MAT_VEC_MULT
+    boost::mpl::false_        selector;
+# else
+	mtl::traits::is_static<Matrix> selector;
+# endif
+    dense_mat_cvec_mult(A, v, w, Assign(), selector);
+}
+
+// Element structure vector multiplication
+template <typename Matrix, typename VectorIn, typename VectorOut, typename Assign>
+inline void mat_cvec_mult(const Matrix& A, const VectorIn& v, VectorOut& w, Assign, tag::flat<tag::element_structure>)
+{
+    vampir_trace<3048> tracer;
+    if (mtl::vector::size(w) == 0) return;
+
+    typedef typename Collection<VectorOut>::value_type value_type;
+    typedef typename Collection<VectorIn>::value_type  value_in_type;
+
+    value_in_type varray[1024];
+    value_type    warray[1024];
+
+    if (Assign::init_to_zero) set_to_zero(w);
+    for(int elmi= 0; elmi < A.m_total_elements; elmi++){
+	const typename Matrix::element_type& elementi= A.m_elements[elmi];
+	const typename Matrix::element_type::index_type& indices= elementi.get_indices();
+	unsigned int n(size(indices));
+
+	if (n <= 1024) {
+	    VectorIn vtmp(n, varray);
+	    for (unsigned int i= 0; i < n; i++)
+		vtmp[i]= v[indices[i]];
+	    VectorOut wtmp(n, warray);
+	    wtmp= elementi.get_values() * vtmp;
+	    for (unsigned int i= 0; i < n; i++)
+		Assign::update(w[indices[i]], wtmp[i]);
+	} else {
+	    VectorIn vtmp(n);
+	    for (unsigned int i= 0; i < n; i++)
+		vtmp[i]= v[indices[i]];
+	    VectorOut wtmp(elementi.get_values() * vtmp);
+	    for (unsigned int i= 0; i < n; i++)
+		Assign::update(w[indices[i]], wtmp[i]);
+	}
+    }
+}
+
+// Multi-vector vector multiplication (tag::multi_vector is derived from dense)
+template <typename Matrix, typename VectorIn, typename VectorOut, typename Assign>
+inline void mat_cvec_mult(const Matrix& A, const VectorIn& v, VectorOut& w, Assign, tag::flat<tag::multi_vector>)
+{
+    vampir_trace<3019> tracer;
+    if (Assign::init_to_zero) set_to_zero(w);
+    for (unsigned i= 0; i < num_cols(A); i++)
+	Assign::update(w, A.vector(i) * v[i]);
+}
+
+// Transposed multi-vector vector multiplication (tag::transposed_multi_vector is derived from dense)
+template <typename TransposedMatrix, typename VectorIn, typename VectorOut, typename Assign>
+inline void mat_cvec_mult(const TransposedMatrix& A, const VectorIn& v, VectorOut& w, Assign, tag::flat<tag::transposed_multi_vector>)
+{
+    vampir_trace<3020> tracer;
+    typename TransposedMatrix::const_ref_type B= A.ref; // Referred matrix
+
+    if (Assign::init_to_zero) set_to_zero(w);
+    for (unsigned i= 0; i < num_cols(B); i++)
+	Assign::update(w[i], dot_real(B.vector(i), v));
+}
+
+// Hermitian multi-vector vector multiplication (tag::hermitian_multi_vector is derived from dense)
+template <typename HermitianMatrix, typename VectorIn, typename VectorOut, typename Assign>
+inline void mat_cvec_mult(const HermitianMatrix& A, const VectorIn& v, VectorOut& w, Assign, tag::flat<tag::hermitian_multi_vector>)
+{
+    vampir_trace<3021> tracer;
+    typename HermitianMatrix::const_ref_type B= A.const_ref(); // Referred matrix
+
+    if (Assign::init_to_zero) set_to_zero(w);
+    for (unsigned i= 0; i < num_cols(B); i++)
+	Assign::update(w[i], dot(B.vector(i), v));
+}
+
+
+
+// Sparse row-major matrix vector multiplication
+template <typename Matrix, typename VectorIn, typename VectorOut, typename Assign>
+inline void smat_cvec_mult(const Matrix& A, const VectorIn& v, VectorOut& w, Assign, tag::row_major)
+{
+    vampir_trace<3022> tracer;
+    using namespace tag;
+    using mtl::traits::range_generator;
+    using math::zero;
+    using mtl::vector::set_to_zero;
+
+    typedef typename range_generator<row, Matrix>::type       a_cur_type;
+    typedef typename range_generator<nz, a_cur_type>::type    a_icur_type;
+    typename mtl::traits::col<Matrix>::type                   col_a(A);
+    typename mtl::traits::const_value<Matrix>::type           value_a(A);
+
+    if (Assign::init_to_zero) set_to_zero(w);
+
+    typedef typename Collection<VectorOut>::value_type        value_type;
+    a_cur_type ac= begin<row>(A), aend= end<row>(A);
+    for (unsigned i= 0; ac != aend; ++ac, ++i) {
+	value_type tmp= zero(w[i]);
+	for (a_icur_type aic= begin<nz>(ac), aiend= end<nz>(ac); aic != aiend; ++aic)
+	    tmp+= value_a(*aic) * v[col_a(*aic)];
+	Assign::update(w[i], tmp);
+    }
+}
+
+// Row-major compressed2D with very few entries (i.e. Very Sparse MATrix) times vector
+template <typename MValue, typename MPara, typename VectorIn, typename VectorOut, typename Assign>
+inline void vsmat_cvec_mult(const compressed2D<MValue, MPara>& A, const VectorIn& v, VectorOut& w, Assign, tag::row_major)
+{
+    vampir_trace<3064> tracer;
+    using math::zero;
+
+    typedef compressed2D<MValue, MPara>                       Matrix;
+    typedef typename Collection<Matrix>::size_type            size_type;
+    typedef typename Collection<VectorOut>::value_type        value_type;
+
+    if (size(w) == 0) return;
+    const value_type z(math::zero(w[0]));
+
+    // std::cout << "very sparse: nnz = " << A.nnz() << ", num_rows = " << num_rows(A) << '\n';
+
+    size_type nr= num_rows(A);
+    for (size_type i1= 0, i2= std::min<size_type>(1024, nr); i1 < i2; i1= i2, i2= std::min<size_type>(i2 + 1024, nr)) {
+	// std::cout << "range = " << i1 << " .. " << i2 << "\n";
+	if (A.ref_major()[i1] < A.ref_major()[i2])
+	    for (size_type i= i1; i < i2; ++i) {
+		const size_type cj0= A.ref_major()[i], cj1= A.ref_major()[i+1];
+		value_type      tmp0(z);
+		for (size_type j0= cj0; j0 != cj1; ++j0)
+		    tmp0+= A.data[j0] * v[A.ref_minor()[j0]];
+		Assign::first_update(w[i], tmp0);
+	    }
+    }
+}
+
+#ifdef MTL_CRS_CVEC_MULT_TUNING
+template <unsigned Index, unsigned BSize, typename SizeType>
+struct crs_cvec_mult_block
+{
+    template <typename Matrix, typename VectorIn, typename CBlock, typename TBlock>
+    void operator()(const Matrix& A, const VectorIn& v, const CBlock& cj, TBlock& tmp) const
+    {
+	for (SizeType j= cj.value; j != cj.sub.value; ++j) // cj is one index larger
+	    tmp.value+= A.data[j] * v[A.ref_minor()[j]];
+	sub(A, v, cj.sub, tmp.sub);
+    }
+
+    template <typename VectorOut, typename TBlock, typename Assign>
+    void first_update(VectorOut& w, SizeType i, const TBlock& tmp, Assign as) const
+    {
+	Assign::first_update(w[i + Index], tmp.value);
+	sub.first_update(w, i, tmp.sub, as);
+    }
+
+    crs_cvec_mult_block<Index+1, BSize, SizeType> sub;
+};
+
+
+template <unsigned BSize, typename SizeType>
+struct crs_cvec_mult_block<BSize, BSize, SizeType>
+{
+    template <typename Matrix, typename VectorIn, typename CBlock, typename TBlock>
+    void operator()(const Matrix& A, const VectorIn& v, const CBlock& cj, TBlock& tmp) const
+    {
+	for (SizeType j= cj.value; j != cj.sub.value; ++j)// cj is one index larger
+	    tmp.value+= A.data[j] * v[A.ref_minor()[j]];
+    }
+
+    template <typename VectorOut, typename TBlock, typename Assign>
+    void first_update(VectorOut& w, SizeType i, const TBlock& tmp, Assign) const
+    {
+	Assign::first_update(w[i + BSize], tmp.value);
+    }
+};
+
+
+// Row-major compressed2D vector multiplication
+template <unsigned BSize, typename MValue, typename MPara, typename VectorIn, typename VectorOut, typename Assign>
+inline void smat_cvec_mult(const compressed2D<MValue, MPara>& A, const VectorIn& v, VectorOut& w, Assign as, tag::row_major)
+{
+    vampir_trace<3049> tracer;
+    using math::zero;
+
+    if (A.nnz() < num_rows(A)) {
+	vsmat_cvec_mult(A, v, w, as, tag::row_major());
+	return;
+    }
+
+    typedef compressed2D<MValue, MPara>                       Matrix;
+    typedef typename Collection<VectorOut>::value_type        value_type;
+    typedef typename mtl::traits::omp_size_type<typename Collection<Matrix>::size_type>::type size_type;
+
+    if (size(w) == 0) return;
+    const value_type z(math::zero(w[0]));
+
+    size_type nr= num_rows(A), nrb= nr / BSize * BSize;
+
+    #ifdef MTL_WITH_OPENMP
+    #   pragma omp parallel
+    #endif
+    {
+    	#ifdef MTL_WITH_OPENMP
+	    vampir_trace<8004> tracer;
+    	#   pragma omp for
+    	#endif
+	for (size_type i= 0; i < nrb; i+= BSize) {
+	    multi_constant_from_array<0, BSize+1, size_type> cj(A.ref_major(), i);
+	    multi_tmp<BSize, value_type>                     tmp(z);
+	    crs_cvec_mult_block<0, BSize-1, size_type>       block;
+
+	    block(A, v, cj, tmp);
+	    block.first_update(w, i, tmp, as);
+	}
+    }
+
+    for (size_type i= nrb; i < nr; ++i) {
+	const size_type cj0= A.ref_major()[i], cj1= A.ref_major()[i+1];
+	value_type      tmp0(z);
+	for (size_type j0= cj0; j0 != cj1; ++j0)
+	    tmp0+= A.data[j0] * v[A.ref_minor()[j0]];
+	Assign::first_update(w[i], tmp0);
+    }
+}
+
+template <typename MValue, typename MPara, typename VectorIn, typename VectorOut, typename Assign>
+typename mtl::traits::enable_if_scalar<typename Collection<VectorOut>::value_type>::type
+inline smat_cvec_mult(const compressed2D<MValue, MPara>& A, const VectorIn& v, VectorOut& w, Assign, tag::row_major)
+{
+    smat_cvec_mult<crs_cvec_mult_block_size>(A, v, w, Assign(), tag::row_major());
+}
+#endif
+
+
+#if !defined(MTL_CRS_CVEC_MULT_NO_ACCEL) && !defined(MTL_CRS_CVEC_MULT_TUNING)
+
+template <typename MValue, typename MPara, typename VectorIn, typename VectorOut, typename Assign>
+typename mtl::traits::enable_if_scalar<typename Collection<VectorOut>::value_type>::type
+inline adapt_crs_cvec_mult(const compressed2D<MValue, MPara>& A, const VectorIn& v, VectorOut& w, Assign)
+{
+    vampir_trace<3065> tracer;
+    using math::zero;
+    assert(!Assign::init_to_zero);
+
+    typedef compressed2D<MValue, MPara>                       Matrix;
+    typedef typename Collection<Matrix>::size_type            size_type;
+    typedef typename Collection<VectorOut>::value_type        value_type;
+
+    const value_type z(math::zero(w[0]));
+    size_type nr= num_rows(A), nrb= nr / 4 * 4, nrb2= nr / 64 * 64;
+
+    for (size_type i1= 0; i1 < nrb2; i1+= 64)
+	if (A.ref_major()[i1] != A.ref_major()[i1 + 64])
+	    for (size_type i2= i1, i2e= i1+64; i2 < i2e; i2+= 16)
+		if (A.ref_major()[i2] != A.ref_major()[i2 + 16])
+		    for (size_type i= i2, i3e= i2+16; i < i3e; i+= 4)
+			if (A.ref_major()[i] != A.ref_major()[i + 4]) {
+			    const size_type cj0= A.ref_major()[i], cj1= A.ref_major()[i+1], cj2= A.ref_major()[i+2],
+				cj3= A.ref_major()[i+3], cj4= A.ref_major()[i+4];
+			    value_type      tmp0(z), tmp1(z), tmp2(z), tmp3(z);
+			    for (size_type j0= cj0; j0 != cj1; ++j0)
+				tmp0+= A.data[j0] * v[A.ref_minor()[j0]];
+			    for (size_type j1= cj1; j1 != cj2; ++j1)
+				tmp1+= A.data[j1] * v[A.ref_minor()[j1]];
+			    for (size_type j2= cj2; j2 != cj3; ++j2)
+				tmp2+= A.data[j2] * v[A.ref_minor()[j2]];
+			    for (size_type j3= cj3; j3 != cj4; ++j3)
+				tmp3+= A.data[j3] * v[A.ref_minor()[j3]];
+
+			    Assign::first_update(w[i], tmp0);
+			    Assign::first_update(w[i+1], tmp1);
+			    Assign::first_update(w[i+2], tmp2);
+			    Assign::first_update(w[i+3], tmp3);
+			}
+
+    for (size_type i= nrb2; i < nrb; i+= 4)
+	if (A.ref_major()[i] != A.ref_major()[i + 4])  {
+	    const size_type cj0= A.ref_major()[i], cj1= A.ref_major()[i+1], cj2= A.ref_major()[i+2],
+		cj3= A.ref_major()[i+3], cj4= A.ref_major()[i+4];
+	    value_type      tmp0(z), tmp1(z), tmp2(z), tmp3(z);
+	    for (size_type j0= cj0; j0 != cj1; ++j0)
+		tmp0+= A.data[j0] * v[A.ref_minor()[j0]];
+	    for (size_type j1= cj1; j1 != cj2; ++j1)
+		tmp1+= A.data[j1] * v[A.ref_minor()[j1]];
+	    for (size_type j2= cj2; j2 != cj3; ++j2)
+		tmp2+= A.data[j2] * v[A.ref_minor()[j2]];
+	    for (size_type j3= cj3; j3 != cj4; ++j3)
+		tmp3+= A.data[j3] * v[A.ref_minor()[j3]];
+
+	    Assign::first_update(w[i], tmp0);
+	    Assign::first_update(w[i+1], tmp1);
+	    Assign::first_update(w[i+2], tmp2);
+	    Assign::first_update(w[i+3], tmp3);
+	}
+
+    for (size_type i= nrb; i < nr; ++i) {
+	const size_type cj0= A.ref_major()[i], cj1= A.ref_major()[i+1];
+	value_type      tmp0(z);
+	for (size_type j0= cj0; j0 != cj1; ++j0)
+	    tmp0+= A.data[j0] * v[A.ref_minor()[j0]];
+	Assign::first_update(w[i], tmp0);
+    }
+}
+
+// Row-major compressed2D vector multiplication
+template <typename MValue, typename MPara, typename VectorIn, typename VectorOut, typename Assign>
+typename mtl::traits::enable_if_scalar<typename Collection<VectorOut>::value_type>::type
+inline smat_cvec_mult(const compressed2D<MValue, MPara>& A, const VectorIn& v, VectorOut& w, Assign as, tag::row_major)
+{
+    vampir_trace<3049> tracer;
+    // vampir_trace<5056> tttracer;
+    using math::zero;
+
+    if (A.nnz() < num_rows(A) && !as.init_to_zero) {
+	vsmat_cvec_mult(A, v, w, as, tag::row_major());
+	return;
+    }
+
+    typedef compressed2D<MValue, MPara>                       Matrix;
+    typedef typename Collection<VectorOut>::value_type        value_type;
+    typedef typename mtl::traits::omp_size_type<typename Collection<Matrix>::size_type>::type size_type;
+
+    if (size(w) == 0) return;
+    const value_type z(math::zero(w[0]));
+
+    size_type nr= num_rows(A), nrb= nr / 4 * 4;
+    if (nr > 10) {
+	size_type nh= nr / 2, nq= nr / 4, nt= nr - nq;
+	if (!as.init_to_zero &&
+	    (A.ref_major()[1] == A.ref_major()[0]
+	     || A.ref_major()[nq] == A.ref_major()[nq+1]
+	     || A.ref_major()[nh] == A.ref_major()[nh+1]
+	     || A.ref_major()[nt] == A.ref_major()[nt+1]
+	     || A.ref_major()[nr-1] == A.ref_major()[nr])) {
+	    adapt_crs_cvec_mult(A, v, w, as);
+	    return;
+	}
+    }
+
+    #ifdef MTL_WITH_OPENMP
+    #   pragma omp parallel
+    #endif
+    {
+    	#ifdef MTL_WITH_OPENMP
+	    vampir_trace<8004> tracer;
+    	#   pragma omp for
+    	#endif
+	    for (size_type i= 0; i < nrb; i+= 4) {
+		const size_type cj0= A.ref_major()[i], cj1= A.ref_major()[i+1], cj2= A.ref_major()[i+2],
+		    cj3= A.ref_major()[i+3], cj4= A.ref_major()[i+4];
+		value_type      tmp0(z), tmp1(z), tmp2(z), tmp3(z);
+		for (size_type j0= cj0; j0 != cj1; ++j0)
+		    tmp0+= A.data[j0] * v[A.ref_minor()[j0]];
+		for (size_type j1= cj1; j1 != cj2; ++j1)
+		    tmp1+= A.data[j1] * v[A.ref_minor()[j1]];
+		for (size_type j2= cj2; j2 != cj3; ++j2)
+		    tmp2+= A.data[j2] * v[A.ref_minor()[j2]];
+		for (size_type j3= cj3; j3 != cj4; ++j3)
+		    tmp3+= A.data[j3] * v[A.ref_minor()[j3]];
+
+		Assign::first_update(w[i], tmp0);
+		Assign::first_update(w[i+1], tmp1);
+		Assign::first_update(w[i+2], tmp2);
+		Assign::first_update(w[i+3], tmp3);
+	    }
+    }
+
+    for (size_type i= nrb; i < nr; ++i) {
+	const size_type cj0= A.ref_major()[i], cj1= A.ref_major()[i+1];
+	value_type      tmp0(z);
+	for (size_type j0= cj0; j0 != cj1; ++j0)
+	    tmp0+= A.data[j0] * v[A.ref_minor()[j0]];
+	Assign::first_update(w[i], tmp0);
+    }
+}
+#endif
+
+// Row-major ell_matrix vector multiplication
+template <typename MValue, typename MPara, typename VectorIn, typename VectorOut, typename Assign>
+typename mtl::traits::enable_if_scalar<typename Collection<VectorOut>::value_type>::type
+inline smat_cvec_mult(const ell_matrix<MValue, MPara>& A, const VectorIn& v, VectorOut& w, Assign, tag::row_major)
+{
+    typedef typename MPara::size_type size_type;
+
+    const size_type stride= A.stride(), slots= A.slots();
+    for (size_type r= 0; r < A.dim1(); ++r) {
+	MValue s(0);
+	for (size_type k= r, i= 0; i < slots; ++i, k+= stride)
+	    s+= A.ref_data()[k] * v[A.ref_minor()[k]];
+	Assign::first_update(w[r], s);
+    }
+ }
+
+
+// Row-major sparse_banded vector multiplication
+template <typename MValue, typename MPara, typename VectorIn, typename VectorOut, typename Assign>
+typename mtl::traits::enable_if_scalar<typename Collection<VectorOut>::value_type>::type
+inline smat_cvec_mult(const sparse_banded<MValue, MPara>& A, const VectorIn& v, VectorOut& w, Assign, tag::row_major)
+{
+    vampir_trace<3069> tracer;
+    typedef sparse_banded<MValue, MPara>                      Matrix;
+    typedef typename Collection<VectorOut>::value_type        value_type;
+    typedef typename Matrix::band_size_type                   band_size_type;
+    typedef typename MPara::size_type                         size_type;
+    typedef mtl::vector::dense_vector<band_size_type, vector::parameters<> > vector_type;
+
+    if (size(w) == 0) return;
+    const value_type z(math::zero(w[0]));
+
+    size_type nr= num_rows(A), nc= num_cols(A), nb= A.ref_bands().size();
+    if (nb == size_type(0) && Assign::init_to_zero) {
+	set_to_zero(w);
+	return;
+    }
+
+    vector_type bands(A.ref_bands()), begin_rows(max(0, -bands)), end_rows(min(nr, nc - bands));
+    assert(end_rows[nb-1] > 0);
+
+    // std::cout << "bands = " << bands << ", begin_rows = " << begin_rows << ", end_rows = " << end_rows << "\n";
+    size_type begin_pos= 0, end_pos= nb - 1;
+
+    // find lowest diagonal in row 0
+    while (begin_pos < nb && begin_rows[begin_pos] > 0) begin_pos++;
+    // if at the end, the first rows are empty
+    if (begin_pos == nb && Assign::init_to_zero) {
+	w[irange(begin_rows[--begin_pos])]= z;
+	// std::cout << "w[0.." << begin_rows[begin_pos] << "] <- 0\n";
+    }
+
+    band_size_type from= begin_rows[begin_pos];
+    // find first entry with same value
+    while (begin_pos > 0 && begin_rows[begin_pos - 1] == from) {
+	assert(from == 0); // should only happen when multiple bands start in row 0
+	begin_pos--;
+    }
+    for (bool active= true; active; ) {
+	// search backwards for the next-largest entry
+	band_size_type to= begin_pos > 0 && begin_rows[begin_pos - 1] <= end_rows[end_pos] ? begin_rows[begin_pos - 1] : end_rows[end_pos];
+
+	// std::cout << "rows " << from << ".." << to << ": with bands ";
+	// for (size_type i= begin_pos; i <= end_pos; i++)
+	//     std::cout << bands[i] << (i < end_pos ? ", " : "\n");
+
+	const MValue* Aps= A.ref_data() + (from * nb + begin_pos);
+
+	const band_size_type blocked_to= ((to - from) & -4) + from;
+	assert((blocked_to - from) % 4 == 0 && blocked_to >= band_size_type(from) && blocked_to <= band_size_type(to));
+	for (band_size_type r= from; r < blocked_to; r+= 4) {
+	    value_type     tmp0(z), tmp1(z), tmp2(z), tmp3(z);
+	    const MValue   *Ap0= Aps, *Ap1= Aps + nb, *Ap2= Ap1 + nb, *Ap3= Ap2 + nb;
+	    for (size_type b= begin_pos; b <= end_pos; ++b, ++Ap0, ++Ap1, ++Ap2, ++Ap3) {
+		tmp0+= *Ap0 * v[r + bands[b]];
+		tmp1+= *Ap1 * v[r + bands[b] + 1];
+		tmp2+= *Ap2 * v[r + bands[b] + 2];
+		tmp3+= *Ap3 * v[r + bands[b] + 3];
+	    }
+	    Assign::first_update(w[r], tmp0);
+	    Assign::first_update(w[r+1], tmp1);
+	    Assign::first_update(w[r+2], tmp2);
+	    Assign::first_update(w[r+3], tmp3);
+	    Aps+= 4 * nb;
+	}
+
+	for (band_size_type r= blocked_to; r < band_size_type(to); r++) {
+	    value_type     tmp(z);
+	    const MValue*  Ap= Aps;
+	    for (size_type b= begin_pos; b <= end_pos; ++b, ++Ap)
+		tmp+= *Ap * v[r + bands[b]];
+	    Assign::first_update(w[r], tmp);
+	    Aps+= nb;
+	}
+
+	if (begin_pos > 0) {
+	    if (begin_rows[begin_pos-1] == to)
+		begin_pos--;
+	    if (end_rows[end_pos] == to)
+		end_pos--;
+	} else { // begin == 0 -> decrement end_pos or finish
+	    if (end_rows[0] == to) {
+		active= false;
+		assert(end_rows[0] = end_rows[end_pos]);
+	    } else {
+		assert(end_pos > 0);
+		end_pos--;
+	    }
+	}
+	assert(begin_pos <= end_pos);
+	from= to;
+    }
+
+    if (size_type(end_rows[0]) < nr  && Assign::init_to_zero) {
+	w[irange(end_rows[0], nr)]= z;
+	// std::cout << "w[" << end_rows[0] << ".." << nr << "] <- 0\n";
+    }
+}
+
+// Sparse column-major matrix vector multiplication
+template <typename Matrix, typename VectorIn, typename VectorOut, typename Assign>
+inline void smat_cvec_mult(const Matrix& A, const VectorIn& v, VectorOut& w, Assign, tag::col_major)
+{
+    vampir_trace<3023> tracer;
+    using namespace tag; namespace traits = mtl::traits;
+    using traits::range_generator;
+    using mtl::vector::set_to_zero;
+    typedef typename range_generator<col, Matrix>::type       a_cur_type;
+    typedef typename range_generator<nz, a_cur_type>::type    a_icur_type;
+
+    typename traits::row<Matrix>::type                        row_a(A);
+    typename traits::const_value<Matrix>::type                value_a(A);
+
+    if (Assign::init_to_zero) set_to_zero(w);
+
+    unsigned rv= 0; // traverse all rows of v
+    for (a_cur_type ac= begin<col>(A), aend= end<col>(A); ac != aend; ++ac, ++rv) {
+	typename Collection<VectorIn>::value_type    vv= v[rv];
+	for (a_icur_type aic= begin<nz>(ac), aiend= end<nz>(ac); aic != aiend; ++aic)
+	    Assign::update(w[row_a(*aic)], value_a(*aic) * vv);
+    }
+}
+
+// Sparse matrix vector multiplication
+template <typename Matrix, typename VectorIn, typename VectorOut, typename Assign>
+inline void mat_cvec_mult(const Matrix& A, const VectorIn& v, VectorOut& w, Assign, tag::flat<tag::sparse>)
+{
+    smat_cvec_mult(A, v, w, Assign(), typename OrientedCollection<Matrix>::orientation());
+}
+
+
+
+}} // namespace mtl::matrix
+
+
+
+
+#endif // MTL_MAT_VEC_MULT_INCLUDE
+

AMDiS/lib/mtl4/boost/numeric/mtl/operation/matrix_bracket.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_MATRIX_BRACKETS_INCLUDE
+#define MTL_MATRIX_BRACKETS_INCLUDE
+
+#include <algorithm>
+#include <boost/utility/enable_if.hpp>
+#include <boost/type_traits/is_integral.hpp>
+#include <boost/type_traits/is_same.hpp>
+#include <boost/type_traits/remove_reference.hpp>
+#include <boost/mpl/bool.hpp>
+#include <boost/numeric/mtl/mtl_fwd.hpp>
+#include <boost/numeric/mtl/utility/irange.hpp>
+#include <boost/numeric/mtl/utility/iset.hpp>
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/operation/column_in_matrix.hpp>
+#include <boost/numeric/mtl/operation/row_in_matrix.hpp>
+
+namespace mtl { namespace operations {
+
+    template <typename Matrix, typename Ref, typename ValueRef>
+    struct bracket_proxy
+    {
+	typedef typename Matrix::value_type      value_type;
+	typedef typename Matrix::size_type       size_type;
+	typedef RowInMatrix<typename boost::remove_reference<Ref>::type> row_traits;
+
+	explicit bracket_proxy(Ref matrix, size_type row) : matrix(matrix), row(row) {}
+
+	ValueRef operator[] (size_type col) { return matrix(row, col);	}
+
+	template <typename T> struct my_traits { static const bool value= boost::is_same<T, mtl::irange>::value && row_traits::exists; };
+
+	template <typename T>
+	typename boost::lazy_enable_if_c<my_traits<T>::value, row_traits>::type
+	operator[] (const T& col_range)
+	{
+	    return row_traits::apply(matrix, row, col_range);
+	}
+      protected:
+	Ref         matrix;
+	size_type   row;
+    };
+
+
+
+    template <typename Matrix, typename Ref, typename ValueRef>
+    struct range_bracket_proxy
+    {
+	typedef typename Matrix::size_type       size_type;
+	typedef ColumnInMatrix<typename boost::remove_reference<Ref>::type> col_traits;
+
+	explicit range_bracket_proxy(Ref matrix, const irange& row_range) : matrix(matrix), row_range(row_range) {}
+
+	ValueRef operator[] (const irange& col_range)
+	{
+	    return sub_matrix(matrix, row_range.start(), row_range.finish(),
+			      col_range.start(), col_range.finish());
+	}
+
+	template <typename T> struct my_traits { static const bool value = boost::is_integral<T>::value && col_traits::exists; };
+
+	template <typename T>
+	typename boost::lazy_enable_if_c<my_traits<T>::value, col_traits>::type
+	operator[] (T col)  { return col_traits::apply(matrix, row_range, col); }
+
+      protected:
+	Ref         matrix;
+	irange      row_range;
+    };
+
+    template <typename Matrix, typename Ref, typename ValueRef>
+    struct set_bracket_proxy
+    {
+	set_bracket_proxy(Ref matrix, const iset& row_set) : matrix(matrix), row_set(row_set) {}
+
+	mtl::matrix::indirect<Matrix> operator[] (const iset& col_set)
+	{
+	    return mtl::matrix::indirect<Matrix>(matrix, row_set, col_set);
+	}
+
+      protected:
+	Ref         matrix;
+	iset        row_set;
+     };
+
+} // namespace operations
+
+} // NAMESPACE mtl
+
+#endif // MTL_MATRIX_BRACKETS_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/max.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_MAX_INCLUDE
+#define MTL_MAX_INCLUDE
+
+#include <iostream>
+#include <cmath>
+
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/utility/tag.hpp>
+#include <boost/numeric/mtl/utility/category.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 {
+
+    namespace impl {
+
+	// Do we really need this for matrices?
+	// Then we need a different dispatching
+
+	template <unsigned long Unroll, typename Vector>
+	typename Collection<Vector>::value_type
+	inline max(const Vector& vector, tag::vector)
+	{
+	    typedef typename Collection<Vector>::value_type result_type;
+	    return vector::reduction<Unroll, vector::max_functor, result_type>::apply(vector);
+	}
+
+    } // namespace impl
+
+///Returns maximal entry of %vector v
+template <unsigned long Unroll, typename Value>
+typename Collection<Value>::value_type
+inline max(const Value& value)
+{
+    vampir_trace<2010> tracer;
+    return impl::max<Unroll>(value, typename traits::category<Value>::type());
+}
+
+template <typename Value>
+typename Collection<Value>::value_type
+inline max(const Value& value)
+{
+    return max<8>(value);
+}
+
+} // namespace vector
+
+using vector::max;
+
+} // namespace mtl::vector
+
+#endif // MTL_MAX_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/max_abs_pos.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_MATRIX_MAX_ABS_POS_INCLUDE
+#define MTL_MATRIX_MAX_ABS_POS_INCLUDE
+
+#include <utility>
+#include <cmath>
+#include <boost/numeric/linear_algebra/identity.hpp>
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/concept/magnitude.hpp>
+#include <boost/numeric/mtl/utility/range_generator.hpp>
+#include <boost/numeric/mtl/utility/tag.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+namespace mtl {
+
+namespace matrix {
+///Returns pair (row, col) from absolut maximal entry of %matrix A
+    template <typename Matrix>
+    typename mtl::traits::enable_if_matrix<Matrix, std::pair<typename Collection<Matrix>::size_type, typename Collection<Matrix>::size_type> >::type
+    inline max_abs_pos(const Matrix& A)
+    {
+	vampir_trace<3024> tracer;
+	namespace traits = mtl::traits;
+	using std::abs;
+	typedef typename Collection<Matrix>::value_type   value_type;
+	typedef typename Collection<Matrix>::size_type    size_type;
+
+	typename RealMagnitude<value_type>::type max(abs(A[0][0]));
+	size_type r= 0, c= 0;
+
+	typename traits::row<Matrix>::type             row(A);
+	typename traits::col<Matrix>::type             col(A);
+	typename traits::const_value<Matrix>::type     value(A);
+	typedef typename traits::range_generator<tag::major, Matrix>::type  cursor_type;
+
+	for (cursor_type cursor = begin<tag::major>(A), cend = end<tag::major>(A); cursor != cend; ++cursor) {
+	    typedef typename traits::range_generator<tag::nz, cursor_type>::type icursor_type;
+	    for (icursor_type icursor = begin<tag::nz>(cursor), icend = end<tag::nz>(cursor); icursor != icend; ++icursor)
+		if (abs(value(*icursor)) > max) {
+		    max= abs(value(*icursor));
+		    r= row(*icursor);
+		    c= col(*icursor);
+		}
+	}
+
+	return std::make_pair(r, c);
+    }
+
+} // namespace matrix
+
+namespace vector {
+///Returns position from absolut maximal entry of %vector v
+    template <typename Vector>
+    typename mtl::traits::enable_if_vector<Vector, typename Collection<Vector>::size_type>::type
+    inline max_abs_pos(const Vector& v)
+    {
+	vampir_trace<2011> tracer;
+	using std::abs;
+	typedef typename Collection<Vector>::size_type    size_type;
+	typedef typename Collection<Vector>::value_type   value_type;
+
+	size_type i= 0;
+	size_type max_col= size(v);
+	value_type max(abs(v[0]));
+
+	for(size_type j= 1; j < max_col; j++)
+	    if(abs(v[j]) > max) {
+		max = abs(v[j]);
+		i= j;
+	    }
+	return i;
+    }
+
+} // namespace vector
+
+} // namespace mtl
+
+#endif // MTL_MATRIX_MAX_ABS_POS_INCLUDE
+

AMDiS/lib/mtl4/boost/numeric/mtl/operation/max_of_sums.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_MAX_OF_SUMS_INCLUDE
+#define MTL_MAX_OF_SUMS_INCLUDE
+
+#include <boost/numeric/mtl/concept/magnitude.hpp>
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/linear_algebra/operators.hpp>
+#include <boost/numeric/linear_algebra/identity.hpp>
+#include <boost/numeric/mtl/utility/tag.hpp>
+#include <boost/numeric/mtl/utility/range_generator.hpp>
+#include <boost/numeric/mtl/vector/dense_vector.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+#include <numeric>
+#include <cmath>
+
+
+namespace mtl { namespace impl {
+
+// We need property map of the minor index
+template <typename Matrix, typename MinorIndex>
+typename RealMagnitude<typename Collection<Matrix>::value_type>::type
+inline max_of_sums(const Matrix& matrix, bool aligned, MinorIndex minor_index, unsigned dim2)
+{
+    vampir_trace<2012> tracer;
+    using std::max; using std::abs; using math::zero;
+
+    typedef typename Collection<Matrix>::value_type   value_type;
+    typedef typename RealMagnitude<value_type>::type  real_type;
+    real_type ref, my_zero= zero(ref);
+
+    // If matrix is empty then the result is the identity from the default-constructed value
+    if (num_rows(matrix) == 0 || num_cols(matrix) == 0)
+	return my_zero;
+
+    typedef typename traits::range_generator<tag::major, Matrix>::type     cursor_type;
+    typedef typename traits::range_generator<tag::nz, cursor_type>::type   icursor_type;
+    typename traits::const_value<Matrix>::type                             value(matrix);
+
+    if (aligned) {
+	real_type maxv= my_zero;
+	for (cursor_type cursor = begin<tag::major>(matrix), cend = end<tag::major>(matrix); cursor != cend; ++cursor) {
+	    real_type sum= my_zero;
+	    for (icursor_type icursor = begin<tag::nz>(cursor), icend = end<tag::nz>(cursor); icursor != icend; ++icursor)
+		sum+= abs(value(*icursor));
+	    maxv= max(maxv, sum);
+	}
+	return maxv;
+    }
+
+    // If matrix has other orientation, we compute all sums in a vector
+    vector::dense_vector<real_type>   sums(dim2, my_zero);
+    for (cursor_type cursor = begin<tag::major>(matrix), cend = end<tag::major>(matrix); cursor != cend; ++cursor)
+	for (icursor_type icursor = begin<tag::nz>(cursor), icend = end<tag::nz>(cursor); icursor != icend; ++icursor)
+	    sums[minor_index(*icursor)]+= abs(value(*icursor));
+    // replace by mtl::accumulate<8>
+    return std::accumulate(sums.begin(), sums.end(), my_zero, math::max<real_type>());
+}
+
+
+}} // namespace mtl::impl
+
+#endif // MTL_MAX_OF_SUMS_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/max_pos.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_MATRIX_MAX_POS
+#define MTL_MATRIX_MAX_POS
+
+#include <cmath>
+#include <utility>
+#include <boost/numeric/linear_algebra/identity.hpp>
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/utility/pos_type.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/mtl/operation/look_at_each_nonzero.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+#if 0
+#include <boost/numeric/mtl/utility/range_generator.hpp>
+#include <boost/numeric/mtl/utility/tag.hpp>
+#include <boost/numeric/mtl/utility/property_map.hpp>
+#endif
+
+
+
+namespace mtl {
+
+    namespace vector {
+
+	template <typename Vector>
+	struct max_pos_functor
+	{
+	    typedef typename Collection<Vector>::value_type       value_type;
+	    typedef typename mtl::traits::pos_type<Vector>::type  pos_type;
+	    typedef std::pair<value_type, pos_type>               result_type;
+
+	    // initialize with max value and max position
+	    max_pos_functor()
+	    {
+		value.first= math::identity(math::max<value_type>(), value_type()); // minimal value for comparison
+		value.second= math::identity(math::min<pos_type>(), pos_type());    // maximal position to check if changed
+	    }
+
+	    void operator()(const value_type& x, const pos_type& p)
+	    {
+		if (x > value.first)
+		    value= std::make_pair(x, p);
+	    }
+
+	    bool unchanged() const { return value.second == math::identity(math::min<pos_type>(), pos_type()); }
+
+	    result_type  value;
+	};
+	///Returns position of maximal entry of %vector v
+	template <typename Vector>
+	typename max_pos_functor<Vector>::pos_type
+	inline max_pos(const Vector& v)
+	{
+	    vampir_trace<2013> tracer;
+	    max_pos_functor<Vector> f;
+	    look_at_each_nonzero_pos(v, f);
+
+	    MTL_DEBUG_THROW_IF(f.unchanged(), runtime_error("max_pos cannot be applied on empty container"));
+	    return f.value.second;
+	}
+
+    } // namespace vector
+
+    namespace matrix {
+	///Returns pair (row, col) from maximal entry of %matrix A
+	using mtl::vector::max_pos;
+    }
+
+
+#if 0
+namespace matrix {
+///Returns pair (row, col) from maximal entry of %matrix A
+    template <typename Matrix>
+    std::pair<typename Collection<Matrix>::size_type, typename Collection<Matrix>::size_type>
+    inline max_pos(const Matrix& A)
+    {
+	namespace traits = mtl::traits;
+	typedef typename Collection<Matrix>::value_type   value_type;
+	typedef typename Collection<Matrix>::size_type    size_type;
+
+	value_type max(A[0][0]);
+	size_type r= 0, c= 0;
+
+	typename traits::row<Matrix>::type             row(A);
+	typename traits::col<Matrix>::type             col(A);
+	typename traits::const_value<Matrix>::type     value(A);
+	typedef typename traits::range_generator<tag::major, Matrix>::type  cursor_type;
+
+	for (cursor_type cursor = begin<tag::major>(A), cend = end<tag::major>(A); cursor != cend; ++cursor) {
+	    typedef typename traits::range_generator<tag::nz, cursor_type>::type icursor_type;
+	    for (icursor_type icursor = begin<tag::nz>(cursor), icend = end<tag::nz>(cursor); icursor != icend; ++icursor)
+		if (value(*icursor) > max) {
+		    max= value(*icursor);
+		    r= row(*icursor);
+		    c= col(*icursor);
+		}
+	}
+
+	return std::make_pair(r, c);
+}
+
+} // namespace matrix
+
+namespace vector {
+///Returns position from maximal entry of %vector v
+    template <typename Vector>
+    typename Collection<Vector>::size_type
+    inline max_pos(const Vector& v)
+    {
+	typedef typename Collection<Vector>::size_type    size_type;
+	typedef typename Collection<Vector>::value_type   value_type;
+
+	size_type i= 0;
+
+	size_type max_col= size(v);
+	value_type max(v[0]);
+
+	for(size_type j= 1;i < max_col; j++)
+	    if(v[j] > max) {
+		max = v[j];
+		i= j;
+	    }
+	return i;
+    }
+
+} // namespace vector
+
+#endif
+
+
+} // namespace mtl
+
+#endif // MTL_MATRIX_MAX_POS
+

AMDiS/lib/mtl4/boost/numeric/mtl/operation/merge_complex_matrix.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_MATRIX_MERGE_COMPLEX_MATRIX_INCLUDE
+#define MTL_MATRIX_MERGE_COMPLEX_MATRIX_INCLUDE
+
+namespace mtl { namespace matrix {
+
+
+
+
+}} // namespace mtl::matrix
+
+#endif // MTL_MATRIX_MERGE_COMPLEX_MATRIX_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/merge_complex_vector.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_MERGE_COMPLEX_VECTOR_INCLUDE
+#define MTL_VECTOR_MERGE_COMPLEX_VECTOR_INCLUDE
+
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+namespace mtl { namespace vector {
+
+/// Merge two real-valued vectors into one complex-valued vector.
+/** Elements of the complex vector must be constructible from two real elements.
+    Complex vector is resized if its size is 0 otherwise the vectors must have
+    the same length. **/
+template <typename VectorReal, typename VectorImaginary, typename VectorComplex>
+inline void merge_complex_vector(const VectorReal& r, const VectorImaginary& i, VectorComplex& c)
+{
+    vampir_trace<2014> tracer;
+    typedef typename Collection<VectorComplex>::value_type value_type;
+
+    MTL_THROW_IF(size(r) != size(i), incompatible_size());
+    c.checked_change_dim(size(r));
+
+    for (std::size_t j= 0; j < size(r); ++j)
+	c[j]= value_type(r[j], i[j]);
+}
+
+
+
+
+}} // namespace mtl::vector
+
+#endif // MTL_VECTOR_MERGE_COMPLEX_VECTOR_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/min.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_MIN_INCLUDE
+#define MTL_MIN_INCLUDE
+
+#include <iostream>
+#include <cmath>
+
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/utility/tag.hpp>
+#include <boost/numeric/mtl/utility/category.hpp>
+#include <boost/numeric/mtl/vector/reduction.hpp>
+#include <boost/numeric/mtl/vector/reduction_functors.hpp>
+
+
+namespace mtl { namespace vector {
+
+    namespace impl {
+
+	// Do we really need this for matrices?
+	// Then we need a different dispatching
+
+	template <unsigned long Unroll, typename Vector>
+	typename Collection<Vector>::value_type
+	inline min(const Vector& vector, tag::vector)
+	{
+	    typedef typename Collection<Vector>::value_type result_type;
+	    return vector::reduction<Unroll, vector::min_functor, result_type>::apply(vector);
+	}
+
+    } // namespace impl
+
+///Returns minimal value of %vector v
+template <unsigned long Unroll, typename Value>
+typename Collection<Value>::value_type
+inline min(const Value& value)
+{
+    return impl::min<Unroll>(value, typename traits::category<Value>::type());
+}
+
+template <typename Value>
+typename Collection<Value>::value_type
+inline min(const Value& value)
+{
+    return min<8>(value);
+}
+
+} // namespace vector
+
+using vector::min;
+
+} // namespace mtl
+
+#endif // MTL_MIN_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/min_pos.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_MIN_POS_INCLUDE
+#define MTL_MIN_POS_INCLUDE
+
+#include <utility>
+
+#include <boost/numeric/linear_algebra/identity.hpp>
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/utility/pos_type.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/mtl/operation/look_at_each_nonzero.hpp>
+
+
+namespace mtl {
+
+    namespace vector {
+
+	template <typename Vector>
+	struct min_pos_functor
+	{
+	    typedef typename Collection<Vector>::value_type       value_type;
+	    typedef typename mtl::traits::pos_type<Vector>::type  pos_type;
+	    typedef std::pair<value_type, pos_type>               result_type;
+
+	    // initialize with max value and max position
+	    min_pos_functor() : value(math::identity(math::min<result_type>(), result_type())) {}
+
+	    void operator()(const value_type& x, const pos_type& p)
+	    {
+		if (x < value.first)
+		    value= std::make_pair(x, p);
+ 	    }
+
+	    bool unchanged() const { return value.second == math::identity(math::min<pos_type>(), pos_type()); }
+
+	    result_type  value;
+	};
+	///Returns position of minimal entry of %vector v
+	template <typename Vector>
+	typename min_pos_functor<Vector>::pos_type
+	inline min_pos(const Vector& v)
+	{
+	    min_pos_functor<Vector> f;
+	    look_at_each_nonzero_pos(v, f);
+
+	    MTL_DEBUG_THROW_IF(f.unchanged(), runtime_error("min_pos cannot be applied on empty container"));
+	    return f.value.second;
+	}
+
+    } // namespace vector
+
+    namespace matrix {
+
+	using mtl::vector::min_pos;
+    }
+
+
+} // namespace mtl
+
+#endif // MTL_MIN_POS_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/minimal_increase.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_MINIMAL_INCREASE_INCLUDE
+#define MTL_MINIMAL_INCREASE_INCLUDE
+
+#include <limits>
+
+namespace mtl {
+
+/// Increase x minimally: if x == 0 take minimal value, if x > 0 multiply by (1+2eps) otherwise divide by
+template <typename T>
+T inline minimal_increase(const T& x)
+{
+    const T factor= T(1) + T(2) * std::numeric_limits<T>::epsilon();
+    if (x == T(0))
+	return std::numeric_limits<T>::denorm_min();
+    else
+	return x > T(0) ? x * factor : x / factor;
+}
+
+
+
+} // namespace mtl
+
+#endif // MTL_MINIMAL_INCREASE_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/misc.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_MISC_INCLUDE
+#define MTL_MISC_INCLUDE
+
+namespace mtl {
+
+    template <typename T>
+    T square(const T& x)
+    { return x * x; }
+
+}// namespace mtl
+
+
+#endif // MTL_MISC_INCLUDE
+