AMDiS/lib/mtl4/boost/numeric/itl/pc/sub_matrix_pc.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 ITL_PC_SUB_MATRIX_PC_INCLUDE
+#define ITL_PC_SUB_MATRIX_PC_INCLUDE
+
+#include <boost/static_assert.hpp>
+#include <boost/mpl/if.hpp>
+
+#include <boost/numeric/mtl/utility/range_generator.hpp>
+#include <boost/numeric/mtl/utility/property_map.hpp>
+#include <boost/numeric/mtl/vector/dense_vector.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+#include <boost/numeric/itl/pc/solver.hpp>
+
+namespace itl { namespace pc {
+
+/// Class for applying \tparam Preconditioner only on a sub-matrix
+/** Other entries are just copied.
+    Optionally preconditioner can be referred from outside instead of storing it by setting
+    \tparam Store to true.
+**/
+template <typename Preconditioner, typename Matrix, bool Store= true>
+class sub_matrix_pc
+{
+    typedef mtl::vector::dense_vector<bool>                                               tag_type;
+    typedef typename boost::mpl::if_c<Store, Preconditioner, const Preconditioner&>::type pc_type;
+
+    struct matrix_container
+    {
+	matrix_container() : Ap(0) {}
+
+	matrix_container(const tag_type& tags, const Matrix& src)
+	{
+	    using std::size_t;
+	    using namespace mtl;
+
+	    mtl::vector::dense_vector<size_t> perm(size(tags));
+	    size_t n= 0;
+	    for (size_t i= 0; i < size(tags); ++i) {
+		perm[i]= n;
+		if (tags[i])
+		    ++n;
+	    }
+
+	    Ap= new Matrix(n, n);
+
+	    typename traits::row<Matrix>::type             row(src);
+	    typename traits::col<Matrix>::type             col(src);
+	    typename traits::const_value<Matrix>::type     value(src);
+	    typedef typename traits::range_generator<tag::major, Matrix>::type  cursor_type;
+
+	    matrix::inserter<Matrix> ins(*Ap, Ap->nnz() / Ap->dim1());
+
+	    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)
+		    if (tags[row(*icursor)] && tags[col(*icursor)])
+			ins(perm[row(*icursor)], perm[col(*icursor)]) << value(*icursor);
+	    }
+	}
+
+	~matrix_container() { delete Ap; }
+
+	Matrix* Ap;
+    };
+
+    size_t count_entries() const
+    {
+	using mtl::size;
+	size_t n= 0;
+	for (size_t i= 0; i < size(tags); ++i)
+	    if (tags[i]) ++n;
+	return n;
+    }
+
+  public:
+
+    sub_matrix_pc(const tag_type& tags, const Matrix& A)
+      : tags(tags), n(count_entries()), mc(tags, A), P(*mc.Ap)
+    {
+	BOOST_STATIC_ASSERT((Store));
+	delete mc.Ap;
+	mc.Ap= 0;
+    }
+
+#if 0 // to do later
+    sub_matrix_pc(const tag_type& tags, const Preconditioner& P)
+      : tags(tags), P(P)
+    {
+	// check sizes
+    }
+#endif
+
+  private:
+    template <typename VectorIn>
+    VectorIn& create_x0(VectorIn) const
+    {
+	static VectorIn  x0(n);
+	return x0;
+    }
+
+    template <typename VectorOut>
+    VectorOut& create_y0(VectorOut) const
+    {
+	static VectorOut  y0(n);
+	return y0;
+    }
+
+    template <typename VectorIn>
+    void restrict(const VectorIn& x, VectorIn& x0) const
+    {
+	for (size_t i= 0, j= 0; i < size(tags); ++i)
+	    if (tags[i])
+		x0[j++]= x[i];
+    }
+
+    template <typename VectorIn, typename VectorOut>
+    void prolongate(const VectorIn& x, const VectorOut& y0, VectorOut& y) const
+    {
+	for (size_t i= 0, j= 0; i < size(tags); ++i)
+	    if (tags[i])
+		y[i]= y0[j++];
+	    else
+		y[i]= x[i];
+    }
+
+  public:
+    /// Solve Px = y approximately on according sub-system; remaining entries are copied
+    template <typename VectorIn, typename VectorOut>
+    void solve(const VectorIn& x, VectorOut& y) const
+    {
+	mtl::vampir_trace<5056> tracer;
+	y.checked_change_resource(x);
+
+	VectorIn&  x0= create_x0(x);
+	VectorOut& y0= create_y0(y);
+
+	restrict(x, x0);
+	P.solve(x0, y0);
+	// y0= solve(P, x0); // doesn't compile yet for unknown reasons
+	prolongate(x, y0, y);
+    }
+
+    /// Solve Px = y approximately on according sub-system; remaining entries are copied
+    template <typename VectorIn, typename VectorOut>
+    void adjoint_solve(const VectorIn& x, VectorOut& y) const
+    {
+	mtl::vampir_trace<5057> tracer;
+	y.checked_change_resource(x);
+
+	VectorIn&  x0= create_x0(x);
+	VectorOut& y0= create_y0(y);
+
+	restrict(x, x0);
+	P.adjoint_solve(x0, y0);
+	// y0= adjoint_solve(P, x0); // doesn't compile yet for unknown reasons
+	prolongate(x, y0, y);
+    }
+
+  private:
+    tag_type          tags;
+    size_t            n;
+    matrix_container  mc;
+    pc_type           P;
+};
+
+template <typename Preconditioner, typename Matrix, bool Store, typename Vector>
+solver<sub_matrix_pc<Preconditioner, Matrix, Store>, Vector, false>
+inline solve(const sub_matrix_pc<Preconditioner, Matrix, Store>& P, const Vector& x)
+{
+    return solver<sub_matrix_pc<Preconditioner, Matrix, Store>, Vector, false>(P, x);
+}
+
+template <typename Preconditioner, typename Matrix, bool Store, typename Vector>
+solver<sub_matrix_pc<Preconditioner, Matrix, Store>, Vector, true>
+inline adjoint_solve(const sub_matrix_pc<Preconditioner, Matrix, Store>& P, const Vector& x)
+{
+    return solver<sub_matrix_pc<Preconditioner, Matrix, Store>, Vector, true>(P, x);
+}
+
+
+}} // namespace itl::pc
+
+#endif // ITL_PC_SUB_MATRIX_PC_INCLUDE

AMDiS/lib/mtl4/boost/numeric/itl/smoother/gauss_seidel.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 ITL_GAUSS_SEIDEL_INCLUDE
+#define ITL_GAUSS_SEIDEL_INCLUDE
+
+#include <boost/assert.hpp>
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/mtl/utility/is_row_major.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/matrix/compressed2D.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+namespace itl {
+
+/// Gauss-Seidel smoother
+/** Constructor takes references to a matrix and a right-hand side vector.
+    operator() is applied on a vector and changes it in place.
+    Matrix must be square, stored row-major and free of zero entries in the diagonal.
+    Vectors b and x must have the same number of rows as A.
+**/
+template <typename Matrix>
+class gauss_seidel
+{
+    typedef typename mtl::Collection<Matrix>::value_type Scalar;
+    typedef typename mtl::Collection<Matrix>::size_type  size_type;
+  public:
+    /// Construct with constant references to matrix and RHS vector
+    gauss_seidel(const Matrix& A) : A(A), dia_inv(num_rows(A))
+    {
+	BOOST_STATIC_ASSERT((mtl::traits::is_row_major<Matrix>::value)); // No CCS
+	assert(num_rows(A) == num_cols(A)); // Matrix must be square
+	for (size_type i= 0; i < num_rows(A); ++i) {
+	    Scalar a= A[i][i];
+	    MTL_THROW_IF(a == 0, mtl::missing_diagonal());
+	    dia_inv[i]= 1.0 / a;
+	}
+    }
+
+    /// Apply Gauss-Seidel on vector \p x, i.e. \p x is changed
+    template <typename Vector, typename RHSVector>
+    Vector& operator()(Vector& x, const RHSVector& b) const
+    {
+	mtl::vampir_trace<8551> tracer;
+	namespace tag= mtl::tag; using namespace mtl::traits;
+	using mtl::begin; using mtl::end;
+
+        typedef typename range_generator<tag::row, Matrix>::type       a_cur_type;
+        typedef typename range_generator<tag::nz, a_cur_type>::type    a_icur_type;
+	typename col<Matrix>::type                   col_a(A);
+	typename const_value<Matrix>::type           value_a(A);
+
+	typedef typename mtl::Collection<Vector>::value_type           value_type;
+
+	a_cur_type ac= begin<tag::row>(A), aend= end<tag::row>(A);
+	for (unsigned i= 0; ac != aend; ++ac, ++i) {
+	    value_type tmp= b[i];
+	    for (a_icur_type aic= begin<tag::nz>(ac), aiend= end<tag::nz>(ac); aic != aiend; ++aic)
+		if (col_a(*aic) != i)
+		    tmp-= value_a(*aic) * x[col_a(*aic)];
+	    x[i]= dia_inv[i] * tmp;
+	}
+ 	return x;
+    }
+
+   private:
+    const Matrix&    A;
+    mtl::vector::dense_vector<Scalar>  dia_inv;
+};
+
+
+template <typename Value, typename Parameters>
+class gauss_seidel<mtl::matrix::compressed2D<Value, Parameters> >
+{
+    typedef mtl::matrix::compressed2D<Value, Parameters> Matrix;
+    typedef typename mtl::Collection<Matrix>::value_type Scalar;
+    typedef typename mtl::Collection<Matrix>::size_type  size_type;
+  public:
+    /// Construct with constant references to matrix and RHS vector
+    gauss_seidel(const Matrix& A)
+      : A(A), dia_inv(num_rows(A)), dia_pos(num_rows(A))
+    {
+	BOOST_STATIC_ASSERT((mtl::traits::is_row_major<Matrix>::value)); // No CCS
+	assert(num_rows(A) == num_cols(A)); // Matrix must be square
+	for (size_type i= 0; i < num_rows(A); ++i) {
+	    mtl::utilities::maybe<size_type> pos = A.indexer(A, i, i);
+	    MTL_THROW_IF(!pos, mtl::missing_diagonal());
+	    dia_inv[i]= 1.0 / A.value_from_offset(pos);
+	    dia_pos[i]= pos;
+	}
+    }
+
+    /// Apply Gauss-Seidel on vector \p x, i.e. \p x is changed
+    template <typename Vector, typename RHSVector>
+    Vector& operator()(Vector& x, const RHSVector& b) const
+    {
+	mtl::vampir_trace<8551> tracer;
+	typedef typename mtl::Collection<Vector>::value_type           value_type;
+	typedef typename mtl::Collection<Matrix>::size_type            size_type;
+	const size_type nr= num_rows(A);
+	size_type cj1= A.ref_major()[0];
+	for (size_type i= 0; i < nr; ++i) {
+	    value_type tmp= b[i];
+	    size_type cj0= cj1, cjm= dia_pos[i];
+	    cj1= A.ref_major()[i+1];
+	    for (; cj0 < cjm; cj0++)
+		tmp-= A.data[cj0] * x[A.ref_minor()[cj0]];
+	    for (size_type j= cjm+1; j < cj1; j++)
+		tmp-= A.data[j] * x[A.ref_minor()[j]];
+	    x[i]= dia_inv[i] * tmp;
+	}
+ 	return x;
+    }
+
+
+  private:
+    const Matrix&    A;
+    mtl::vector::dense_vector<Scalar>     dia_inv;
+    mtl::vector::dense_vector<size_type>  dia_pos;
+};
+
+
+} // namespace itl
+
+#endif // ITL_GAUSS_SEIDEL_INCLUDE

AMDiS/lib/mtl4/boost/numeric/itl/smoother/jacobi.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 ITL_JACOBI_INCLUDE
+#define ITL_JACOBI_INCLUDE
+
+#include <boost/assert.hpp>
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/operations.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/mtl/utility/is_row_major.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/matrix/compressed2D.hpp>
+
+namespace itl {
+
+/// Jacobi smoother
+/** Constructor takes references to a matrix and a right-hand side vector.
+    operator() is applied on a vector and changes it in place.
+    Matrix must be square, stored row-major and free of zero entries in the diagonal.
+    Vectors b and x must have the same number of rows as A.
+**/
+template <typename Matrix>
+class jacobi
+{
+    typedef typename mtl::Collection<Matrix>::value_type Scalar;
+    typedef typename mtl::Collection<Matrix>::size_type  size_type;
+  public:
+    /// Construct with constant references to matrix and RHS vector
+    jacobi(const Matrix& A) : A(A), dia_inv(num_rows(A))
+    {
+	BOOST_STATIC_ASSERT((mtl::traits::is_row_major<Matrix>::value)); // No CCS
+	assert(num_rows(A) == num_cols(A)); // Matrix must be square
+	for (size_type i= 0; i < num_rows(A); ++i) {
+	    Scalar a= A[i][i];
+	    MTL_THROW_IF(a == 0, mtl::missing_diagonal());
+	    dia_inv[i]= 1.0 / a;
+	}
+    }
+
+    /// Apply Jacobi on vector \p x, i.e. \p x is changed
+    template <typename Vector, typename RHSVector>
+    Vector& operator()(Vector& x, const RHSVector& b) const
+    {
+	namespace tag= mtl::tag; using namespace mtl::traits;
+	using mtl::begin; using mtl::end; using std::swap;
+
+        typedef typename range_generator<tag::row, Matrix>::type       a_cur_type;
+        typedef typename range_generator<tag::nz, a_cur_type>::type    a_icur_type;
+	typename col<Matrix>::type                   col_a(A);
+	typename const_value<Matrix>::type           value_a(A);
+
+	typedef typename mtl::Collection<Vector>::value_type           value_type;
+
+	static Vector x0(resource(x));
+	a_cur_type ac= begin<tag::row>(A), aend= end<tag::row>(A);
+	for (unsigned i= 0; ac != aend; ++ac, ++i) {
+	    value_type tmp= b[i];
+	    for (a_icur_type aic= begin<tag::nz>(ac), aiend= end<tag::nz>(ac); aic != aiend; ++aic)
+		if (col_a(*aic) != i)
+		    tmp-= value_a(*aic) * x[col_a(*aic)];
+	    x0[i]= dia_inv[i] * tmp;
+	}
+	swap(x0, x);
+ 	return x;
+    }
+
+   private:
+    const Matrix&    A;
+    mtl::vector::dense_vector<Scalar>  dia_inv;
+};
+
+
+template <typename Value, typename Parameters>
+class jacobi<mtl::matrix::compressed2D<Value, Parameters> >
+{
+    typedef mtl::matrix::compressed2D<Value, Parameters> Matrix;
+    typedef typename mtl::Collection<Matrix>::value_type Scalar;
+    typedef typename mtl::Collection<Matrix>::size_type  size_type;
+  public:
+    /// Construct with constant references to matrix and RHS vector
+    jacobi(const Matrix& A)
+      : A(A), dia_inv(num_rows(A)), dia_pos(num_rows(A))
+    {
+	BOOST_STATIC_ASSERT((mtl::traits::is_row_major<Matrix>::value)); // No CCS
+	assert(num_rows(A) == num_cols(A)); // Matrix must be square
+	for (size_type i= 0; i < num_rows(A); ++i) {
+	    mtl::utilities::maybe<size_type> pos = A.indexer(A, i, i);
+	    MTL_THROW_IF(!pos, mtl::missing_diagonal());
+	    dia_inv[i]= 1.0 / A.value_from_offset(pos);
+	    dia_pos[i]= pos;
+	}
+    }
+
+    /// Apply Jacobi on vector \p x, i.e. \p x is changed
+    template <typename Vector, typename RHSVector>
+    Vector& operator()(Vector& x, const RHSVector& b) const
+    {
+	typedef typename mtl::Collection<Vector>::value_type           value_type;
+	typedef typename mtl::Collection<Matrix>::size_type            size_type;
+	const size_type nr= num_rows(A);
+	static Vector x0(resource(x));
+	size_type cj1= A.ref_major()[0];
+	for (size_type i= 0; i < nr; ++i) {
+	    value_type tmp= b[i];
+	    size_type cj0= cj1, cjm= dia_pos[i];
+	    cj1= A.ref_major()[i+1];
+	    for (; cj0 < cjm; cj0++)
+		tmp-= A.data[cj0] * x[A.ref_minor()[cj0]];
+	    for (size_type j= cjm+1; j < cj1; j++)
+		tmp-= A.data[j] * x[A.ref_minor()[j]];
+	    x0[i]= dia_inv[i] * tmp;
+	}
+	swap(x0, x);
+ 	return x;
+    }
+
+
+  private:
+    const Matrix&    A;
+    mtl::vector::dense_vector<Scalar>     dia_inv;
+    mtl::vector::dense_vector<size_type>  dia_pos;
+};
+
+
+} // namespace itl
+
+#endif // ITL_JACOBI_INCLUDE

AMDiS/lib/mtl4/boost/numeric/itl/smoother/jor.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 ITL_JOR_INCLUDE
+#define ITL_JOR_INCLUDE
+
+#include <boost/assert.hpp>
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/mtl/utility/is_row_major.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/matrix/compressed2D.hpp>
+
+#include "./relaxation_parameter.hpp"
+
+namespace itl {
+
+/// Jacobi smoother with relaxation
+/** Constructor takes references to a matrix and a right-hand side vector.
+    operator() is applied on a vector and changes it in place.
+    Matrix must be square, stored row-major and free of zero entries in the diagonal.
+    Vectors b and x must have the same number of rows as A.
+**/
+template <typename Matrix, typename Omega = default_omega>
+class jor
+{
+    typedef typename mtl::Collection<Matrix>::value_type Scalar;
+    typedef typename mtl::Collection<Matrix>::size_type  size_type;
+  public:
+    /// Construct with constant references to matrix and RHS vector
+    jor(const Matrix& A) : A(A), dia_inv(num_rows(A))
+    {
+	BOOST_STATIC_ASSERT((mtl::traits::is_row_major<Matrix>::value)); // No CCS
+	assert(num_rows(A) == num_cols(A)); // Matrix must be square
+	for (size_type i= 0; i < num_rows(A); ++i) {
+	    Scalar a= A[i][i];
+	    MTL_THROW_IF(a == 0, mtl::missing_diagonal());
+	    dia_inv[i]= 1.0 / a;
+	}
+    }
+
+    /// Apply JOR on vector \p x, i.e. \p x is changed
+    template <typename Vector, typename RHSVector>
+    Vector& operator()(Vector& x, const RHSVector& b) const
+    {
+	namespace tag= mtl::tag; using namespace mtl::traits;
+	using mtl::begin; using mtl::end; using std::swap;
+
+        typedef typename range_generator<tag::row, Matrix>::type       a_cur_type;
+        typedef typename range_generator<tag::nz, a_cur_type>::type    a_icur_type;
+	typename col<Matrix>::type                   col_a(A);
+	typename const_value<Matrix>::type           value_a(A);
+
+	typedef typename mtl::Collection<Vector>::value_type           value_type;
+
+	static Vector x0(resource(x));
+	a_cur_type ac= begin<tag::row>(A), aend= end<tag::row>(A);
+	for (unsigned i= 0; ac != aend; ++ac, ++i) {
+	    value_type tmp= b[i];
+	    for (a_icur_type aic= begin<tag::nz>(ac), aiend= end<tag::nz>(ac); aic != aiend; ++aic)
+		if (col_a(*aic) != i)
+		    tmp-= value_a(*aic) * x[col_a(*aic)];
+	    x0[i]= Omega::value * (dia_inv[i] * tmp) + (1. - Omega::value)*x0[i];
+	}
+	swap(x0, x);
+ 	return x;
+    }
+
+   private:
+    const Matrix&    A;
+    mtl::vector::dense_vector<Scalar>  dia_inv;
+};
+
+
+template <typename Value, typename Parameters, typename Omega>
+class jor<mtl::matrix::compressed2D<Value, Parameters> , Omega>
+{
+    typedef mtl::matrix::compressed2D<Value, Parameters> Matrix;
+    typedef typename mtl::Collection<Matrix>::value_type Scalar;
+    typedef typename mtl::Collection<Matrix>::size_type  size_type;
+  public:
+    /// Construct with constant references to matrix and RHS vector
+    jor(const Matrix& A)
+      : A(A), dia_inv(num_rows(A)), dia_pos(num_rows(A))
+    {
+	BOOST_STATIC_ASSERT((mtl::traits::is_row_major<Matrix>::value)); // No CCS
+	assert(num_rows(A) == num_cols(A)); // Matrix must be square
+	for (size_type i= 0; i < num_rows(A); ++i) {
+	    mtl::utilities::maybe<size_type> pos = A.indexer(A, i, i);
+	    MTL_THROW_IF(!pos, mtl::missing_diagonal());
+	    dia_inv[i]= 1.0 / A.value_from_offset(pos);
+	    dia_pos[i]= pos;
+	}
+    }
+
+    /// Apply JOR on vector \p x, i.e. \p x is changed
+    template <typename Vector, typename RHSVector>
+    Vector& operator()(Vector& x, const RHSVector& b) const
+    {
+	typedef typename mtl::Collection<Vector>::value_type           value_type;
+	typedef typename mtl::Collection<Matrix>::size_type            size_type;
+	const size_type nr= num_rows(A);
+	static Vector x0(resource(x));
+	size_type cj1= A.ref_major()[0];
+	for (size_type i= 0; i < nr; ++i) {
+	    value_type tmp= b[i];
+	    size_type cj0= cj1, cjm= dia_pos[i];
+	    cj1= A.ref_major()[i+1];
+	    for (; cj0 < cjm; cj0++)
+		tmp-= A.data[cj0] * x[A.ref_minor()[cj0]];
+	    for (size_type j= cjm+1; j < cj1; j++)
+		tmp-= A.data[j] * x[A.ref_minor()[j]];
+	    x0[i] = Omega::value * (dia_inv[i] * tmp) + (1. - Omega::value)*x0[i];
+	}
+	swap(x0, x);
+ 	return x;
+    }
+
+
+  private:
+    const Matrix&    A;
+    mtl::vector::dense_vector<Scalar>     dia_inv;
+    mtl::vector::dense_vector<size_type>  dia_pos;
+};
+
+
+} // namespace itl
+
+#endif // ITL_JOR_INCLUDE

AMDiS/lib/mtl4/boost/numeric/itl/smoother/relaxation_parameter.hpp

+/*
+ * Marcel Schiffel, 13.10.11
+ *
+ * definition of default relaxation parameter for iterative solvers
+ */
+
+
+#ifndef MTL_ITL_RELAXATION_PARAMETER_HPP
+#define MTL_ITL_RELAXATION_PARAMETER_HPP
+
+
+namespace itl {
+
+
+/**
+ * default relaxation parameter for iterative solvers
+ */
+struct default_omega
+{
+    static const double value;
+};
+
+const double default_omega::value= 2./3.;
+
+} /* namespace itl */
+
+#endif
+

AMDiS/lib/mtl4/boost/numeric/itl/smoother/repeated.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 ITL_REPEATED_INCLUDE
+#define ITL_REPEATED_INCLUDE
+
+namespace itl {
+
+/// Repeat the smoother
+template <typename Smoother, std::size_t N= 1>
+class repeated
+{
+  public:
+
+    typedef Smoother smoother_type;
+
+    /// Construct with \p smoother
+    repeated(const Smoother& smoother) : smoother(smoother) {}
+
+    /// Construct with \p smoother and number of repetitions \p n
+    template <typename Matrix>
+    repeated(const Matrix& A) : smoother(A) {}
+
+    /// Apply smoother n times on vector \p x, i.e. \p x is changed
+    template <typename Vector, typename RHSVector>
+    Vector& operator()(Vector& x, const RHSVector& b) const
+    {
+	for (std::size_t i= 0; i < N; i++)
+	    smoother(x, b);
+ 	return x;
+    }
+
+  private:
+    Smoother    smoother;
+    std::size_t n;
+};
+
+
+} // namespace itl
+
+#endif // ITL_REPEATED_INCLUDE

AMDiS/lib/mtl4/boost/numeric/itl/smoother/sor.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 ITL_SOR_INCLUDE
+#define ITL_SOR_INCLUDE
+
+#include <boost/assert.hpp>
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/mtl/utility/is_row_major.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/matrix/compressed2D.hpp>
+
+#include "./relaxation_parameter.hpp"
+
+namespace itl {
+
+/// Gauss-Seidel smoother with relaxation
+/** Constructor takes references to a matrix and a right-hand side vector.
+    operator() is applied on a vector and changes it in place.
+    Matrix must be square, stored row-major and free of zero entries in the diagonal.
+    Vectors b and x must have the same number of rows as A.
+**/
+template <typename Matrix, typename Omega = default_omega>
+class sor
+{
+    typedef typename mtl::Collection<Matrix>::value_type Scalar;
+    typedef typename mtl::Collection<Matrix>::size_type  size_type;
+
+  public:
+    /// Construct with constant references to matrix and RHS vector
+    sor(const Matrix& A) : A(A), dia_inv(num_rows(A))
+    {
+	BOOST_STATIC_ASSERT((mtl::traits::is_row_major<Matrix>::value)); // No CCS
+	assert(num_rows(A) == num_cols(A)); // Matrix must be square
+	for (size_type i= 0; i < num_rows(A); ++i) {
+	    Scalar a= A[i][i];
+	    MTL_THROW_IF(a == 0, mtl::missing_diagonal());
+	    dia_inv[i]= 1.0 / a;
+	}
+    }
+
+    /// Apply SOR on vector \p x, i.e. \p x is changed
+    template <typename Vector, typename RHSVector>
+    Vector& operator()(Vector& x, const RHSVector& b) const
+    {
+	namespace tag= mtl::tag; using namespace mtl::traits;
+	using mtl::begin; using mtl::end;
+
+        typedef typename range_generator<tag::row, Matrix>::type       a_cur_type;
+        typedef typename range_generator<tag::nz, a_cur_type>::type    a_icur_type;
+	typename col<Matrix>::type                   col_a(A);
+	typename const_value<Matrix>::type           value_a(A);
+
+	typedef typename mtl::Collection<Vector>::value_type           value_type;
+
+	a_cur_type ac= begin<tag::row>(A), aend= end<tag::row>(A);
+	for (unsigned i= 0; ac != aend; ++ac, ++i) {
+	    value_type tmp= b[i];
+	    for (a_icur_type aic= begin<tag::nz>(ac), aiend= end<tag::nz>(ac); aic != aiend; ++aic)
+		if (col_a(*aic) != i)
+		    tmp-= value_a(*aic) * x[col_a(*aic)];
+	    x[i] = Omega::value * (dia_inv[i] * tmp) + (1. - Omega::value)*x[i];
+	}
+ 	return x;
+    }
+
+   private:
+    const Matrix&    A;
+    mtl::vector::dense_vector<Scalar>  dia_inv;
+};
+
+
+template <typename Value, typename Parameters, typename Omega>
+class sor<mtl::matrix::compressed2D<Value, Parameters> , Omega>
+{
+    typedef mtl::matrix::compressed2D<Value, Parameters> Matrix;
+    typedef typename mtl::Collection<Matrix>::value_type Scalar;
+    typedef typename mtl::Collection<Matrix>::size_type  size_type;
+
+  public:
+    /// Construct with constant references to matrix and RHS vector
+    sor(const Matrix& A)
+      : A(A), dia_inv(num_rows(A)), dia_pos(num_rows(A))
+    {
+	BOOST_STATIC_ASSERT((mtl::traits::is_row_major<Matrix>::value)); // No CCS
+	assert(num_rows(A) == num_cols(A)); // Matrix must be square
+	for (size_type i= 0; i < num_rows(A); ++i) {
+	    mtl::utilities::maybe<size_type> pos = A.indexer(A, i, i);
+	    MTL_THROW_IF(!pos, mtl::missing_diagonal());
+	    dia_inv[i]= 1.0 / A.value_from_offset(pos);
+	    dia_pos[i]= pos;
+	}
+    }
+
+    /// Apply SOR on vector \p x, i.e. \p x is changed
+    template <typename Vector, typename RHSVector>
+    Vector& operator()(Vector& x, const RHSVector& b) const
+    {
+	typedef typename mtl::Collection<Vector>::value_type           value_type;
+	typedef typename mtl::Collection<Matrix>::size_type            size_type;
+	const size_type nr= num_rows(A);
+	size_type cj1= A.ref_major()[0];
+	for (size_type i= 0; i < nr; ++i) {
+	    value_type tmp= b[i];
+	    size_type cj0= cj1, cjm= dia_pos[i];
+	    cj1= A.ref_major()[i+1];
+	    for (; cj0 < cjm; cj0++)
+		tmp-= A.data[cj0] * x[A.ref_minor()[cj0]];
+	    for (size_type j= cjm+1; j < cj1; j++)
+		tmp-= A.data[j] * x[A.ref_minor()[j]];
+	    x[i]= Omega::value * (dia_inv[i] * tmp) + (1. - Omega::value)*x[i];
+	}
+ 	return x;
+    }
+
+
+  private:
+    const Matrix&    A;
+    mtl::vector::dense_vector<Scalar>     dia_inv;
+    mtl::vector::dense_vector<size_type>  dia_pos;
+};
+
+
+} // namespace itl
+
+#endif // ITL_GAUSS_SEIDEL_INCLUDE

AMDiS/lib/mtl4/boost/numeric/itl/stepper/armijo.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, Cornelius Steinhardt and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef ITL_ARMIJO_INCLUDE
+#define ITL_ARMIJO_INCLUDE
+
+namespace itl {
+
+/// Step size control by Armijo
+/**
+ **/
+template <typename Value= double>
+class armijo
+{
+  public:
+    typedef Value   value_type;
+
+    // Defaults from Prof. Fischer's lecture
+    armijo(Value delta= 0.5, Value gamma= 0.5, Value beta1= 0.25, Value beta2= 0.5)
+      : delta(delta), gamma(gamma), beta1(beta1), beta((beta1 + beta2) / 2.0) {}
+
+    ///
+    template <typename Vector, typename F, typename Grad>
+    typename mtl::Collection<Vector>::value_type
+    operator() (const Vector& x, const Vector& d, F f, Grad grad_f) const
+    {
+	// Star's step size
+	typename mtl::Collection<Vector>::value_type alpha= -gamma * dot(grad_f(x), d) / dot(d, d);
+	Vector     x_k(x + alpha * d);
+
+	while (f(x_k) > f(x) + (beta1 * alpha) * dot(grad_f(x), d)) {
+	    alpha*= beta;
+	    x_k= x+ alpha * d;
+	}
+	return alpha;
+    }
+  private:
+    Value delta, gamma, beta1, beta;
+};
+
+
+} // namespace itl
+
+#endif // ITL_ARMIJO_INCLUDE

AMDiS/lib/mtl4/boost/numeric/itl/stepper/wolf.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, Cornelius Steinhardt and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef ITL_WOLF_INCLUDE
+#define ITL_WOLF_INCLUDE
+
+namespace itl {
+
+/// Step size control by Wolf
+/**
+ **/
+template <typename Value= double>
+class wolf
+{
+  public:
+    typedef Value   value_type;
+
+    // Defaults from Prof. Fischer's lecture
+    wolf(Value delta= 0.5, Value gamma= 0.5, Value beta1= 0.25, Value beta2= 0.5)
+      : delta(delta), gamma(gamma), beta1(beta1), beta2(beta2) {}
+
+    ///
+    template <typename Vector, typename F, typename Grad>
+    typename mtl::Collection<Vector>::value_type
+    operator() (const Vector& x, const Vector& d, F f, Grad grad_f) const
+    {
+	// Star's step size
+	typename mtl::Collection<Vector>::value_type alpha= -gamma * dot(grad_f(x), d) / dot(d, d);
+	Vector     x_k(x + alpha * d);
+
+	Value beta= (beta1 + beta2) / 2;
+	while (f(x_k) > f(x) + (beta1 * alpha) * dot(grad_f(x), d)
+	       && dot(grad_f(x_k), d) < beta2 * dot(grad_f(x), d)) {
+	    alpha*= beta;
+	    x_k= x+ alpha * d;
+	}
+	return alpha;
+    }
+  private:
+    Value delta, gamma, beta1, beta2;
+};
+
+
+} // namespace itl
+
+#endif // ITL_WOLF_INCLUDE

AMDiS/lib/mtl4/boost/numeric/itl/updater/bfgs.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 ITL_BFGS_INCLUDE
+#define ITL_BFGS_INCLUDE
+
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/itl/utility/exception.hpp>
+
+namespace itl {
+
+/// Update of Hessian matrix for e.g. Quasi-Newton by Broyden, Fletcher, Goldfarb, and Shanno
+struct bfgs
+{
+    /// \f$ H_{k+1}=B_{k+1}^{-1}=(I-\frac{y_k\cdot s_k^T}{y_k^T\cdot s_k})^T\cdot H_k \cdot (I-\frac{y_k\cdot s_k^T}{y_k^T\cdot s_k}) + \frac{s_k\cdot s_k^T}{y_k^T\cdot s_k}\f$
+    template <typename Matrix, typename Vector>
+    void operator() (Matrix& H, const Vector& y, const Vector& s)
+    {
+	typedef typename mtl::Collection<Vector>::value_type value_type;
+	assert(num_rows(H) == num_cols(H));
+
+	value_type gamma= 1 / dot(y,s);
+	MTL_THROW_IF(gamma == 0.0, unexpected_orthogonality());
+	Matrix     A(math::one(H) - gamma * s * trans(y)),
+	           H2(A * H * trans(A) + gamma * s * trans(s));
+	swap(H2, H); // faster than H= H2
+    }
+};
+
+
+
+} // namespace itl
+
+#endif // ITL_BFGS_INCLUDE

AMDiS/lib/mtl4/boost/numeric/itl/updater/broyden.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 ITL_BROYDEN_INCLUDE
+#define ITL_BROYDEN_INCLUDE
+
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/itl/utility/exception.hpp>
+
+namespace itl {
+
+/// Update of Hessian matrix for e.g. Quasi-Newton by Broyden formula
+struct broyden
+{
+    /// \f$ H_{k+1}=B_{k+1}^{-1}=H_k+\frac{(s_k-H_k\cdot y_k)\cdot y_k^T\cdot H_k}{y_k^T\cdot H_k\cdot s_k} \f$
+    template <typename Matrix, typename Vector>
+    void operator() (Matrix& H, const Vector& y, const Vector& s)
+    {
+	typedef typename mtl::Collection<Vector>::value_type value_type;
+	assert(num_rows(H) == num_cols(H));
+
+	Vector     h(H * y), d(s - h);
+	value_type gamma= 1 / dot(y, h);
+	MTL_THROW_IF(gamma == 0.0, unexpected_orthogonality());
+	Matrix     A(gamma * d * trans(y)),
+	           H2(H + A * H);
+	swap(H2, H); // faster than H= H2
+   }
+};
+
+
+
+} // namespace itl
+
+#endif // ITL_BROYDEN_INCLUDE

AMDiS/lib/mtl4/boost/numeric/itl/updater/dfp.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 ITL_DFP_INCLUDE
+#define ITL_DFP_INCLUDE
+
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/itl/utility/exception.hpp>
+
+namespace itl {
+
+/// Update of Hessian matrix for e.g. Quasi-Newton by Davidon, Fletcher and Powell formula
+struct dfp
+{
+    /// \f$ H_{k+1}=B_{k+1}^{-1}=H_k+\frac{s_k\cdot s_k^T}{y_k^T\cdot s_k}- \frac{H_k\cdot y_k\cdot y_k^T\cdot H_k^T}{y_k^TH_k\cdot y_k}\f$
+    template <typename Matrix, typename Vector>
+    void operator() (Matrix& H, const Vector& y, const Vector& s)
+    {
+	typedef typename mtl::Collection<Vector>::value_type value_type;
+	assert(num_rows(H) == num_cols(H));
+
+	Vector h(H*y);
+	value_type gamma= 1 / dot(y,s), alpha= 1 / dot(y,h);
+	MTL_THROW_IF(gamma == 0.0, unexpected_orthogonality());
+	MTL_THROW_IF(alpha == 0.0, unexpected_orthogonality());
+	Matrix     A(alpha * y * trans(y)),
+	           H2(H - H * A * H  + gamma * s * trans(s));
+	swap(H2, H); // faster than H= H2
+   }
+};
+
+
+
+} // namespace itl
+
+#endif // ITL_DFP_INCLUDE

AMDiS/lib/mtl4/boost/numeric/itl/updater/psb.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 ITL_PSB_INCLUDE
+#define ITL_PSB_INCLUDE
+
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/itl/utility/exception.hpp>
+
+namespace itl {
+
+/// Update of Hessian matrix for e.g. Quasi-Newton by Powell's symmetric Broyden formula
+struct psb
+{
+    /// \f$ H_{k+1}=B_{k+1}^{-1}=H_k+\frac{(y_k-H_k\cdot s_k)s_k^T+s_k(y_k-H_k\cdot s_k)^T}{s_k^T\cdot s_k}-\frac{(y_k-H_k\cdot s_k)^T\cdot s_k}{(s_k^ts_k)^2}s_k\cdot s_k^T \f$
+    template <typename Matrix, typename Vector>
+    void operator() (Matrix& H, const Vector& y, const Vector& s)
+    {
+	typedef typename mtl::Collection<Vector>::value_type value_type;
+	assert(num_rows(H) == num_cols(H));
+	Vector     a(s - H * y);
+	value_type gamma= 1 / dot (y, y);
+        MTL_THROW_IF(gamma == 0.0, unexpected_orthogonality());
+
+        H+= gamma * a * trans(y) + gamma * y * trans(a) - dot(a, y) * gamma * gamma * y * trans(y);
+   }
+};
+
+
+
+} // namespace itl
+
+#endif // ITL_PSB_INCLUDE
+

AMDiS/lib/mtl4/boost/numeric/itl/updater/sr1.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 ITL_SR1_INCLUDE
+#define ITL_SR1_INCLUDE
+
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/itl/utility/exception.hpp>
+
+namespace itl {
+
+/// Update of Hessian matrix for e.g. Quasi-Newton by SR1 formulas
+struct sr1
+{
+    /// \f$ H_{k+1}=B_{k+1}^{-1}=H_k+\frac{(s_k-H_k\cdot y_k)(s_k-H_k\cdot y_k)^T}{(s_k-H_k\cdot y_k)^T\cdot y_k} \f$
+    template <typename Matrix, typename Vector>
+    void operator() (Matrix& H, const Vector& y, const Vector& s)
+    {
+	assert(num_rows(H) == num_cols(H));
+	Vector     d(s - H * y);
+	MTL_THROW_IF(dot(d, y) == 0.0, unexpected_orthogonality());
+	H+= 1 / dot(d, y) * d * trans(d);
+   }
+};
+
+
+
+} // namespace itl
+
+#endif // ITL_SR1_INCLUDE

AMDiS/lib/mtl4/boost/numeric/itl/utility/exception.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 ITL_EXCEPTION_INCLUDE
+#define ITL_EXCEPTION_INCLUDE
+
+#include <boost/numeric/mtl/utility/exception.hpp>
+
+namespace itl {
+
+/// Exception for iterative solvers that exhausted the search space, i.e. search direction(s) parallel to already visited Krylov subspace
+/** Either your matrix is too badly conditioned or your termination criterion is too strict for the used arithmetic (maybe use Gnu Multi-precision library). **/
+struct search_space_exhaustion
+  : mtl::runtime_error
+{
+    /// Error can be specified more precisely in constructor if desired
+    explicit search_space_exhaustion(const char *s= "Iterative solvers that exhausted the search space, i.e. search direction(s) parallel to already visited Krylov subspace")
+      : mtl::runtime_error(s) {}
+};
+
+/// Vectors unexpectedly become orthogonal, special case of search_space_exhaustion.
+struct unexpected_orthogonality : search_space_exhaustion {};
+
+} // namespace itl
+
+#endif // ITL_EXCEPTION_INCLUDE

AMDiS/lib/mtl4/boost/numeric/linear_algebra/accumulate.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 MATH_ACCUMULATE_INCLUDE
+#define MATH_ACCUMULATE_INCLUDE
+
+#include <concepts>
+
+#include <boost/numeric/linear_algebra/identity.hpp>
+#include <boost/numeric/linear_algebra/new_concepts.hpp>
+
+namespace math {
+
+// Dispatching between simple and unrolled version
+template <std::InputIterator Iter, std::CopyConstructible Value, typename Op>
+     requires std::Callable2<Op, Value, Value>
+           && std::CopyAssignable<Value, std::Callable2<Op, Value, Value>::result_type>
+Value inline accumulate(Iter first, Iter last, Value init, Op op)
+{
+    // std::cout << "Simple accumulate\n";
+    for (; first != last; ++first)
+	init= op(init, Value(*first));
+    return init;
+}
+
+
+template <std::RandomAccessIterator Iter, std::CopyConstructible Value, typename Op>
+    requires std::Callable2<Op, Value, Value>
+          && std::CopyAssignable<Value, std::Callable2<Op, Value, Value>::result_type>
+          && Commutative<Op, Value>
+          && Monoid<Op, Value>
+          && std::Convertible<Monoid<Op, Value>::identity_result_type, Value>
+Value inline accumulate(Iter first, Iter last, Value init, Op op)
+{
+    // std::cout << "Unrolled accumulate\n";
+    typedef typename std::RandomAccessIterator<Iter>::difference_type difference_type;
+    Value            t0= identity(op, init), t1= identity(op, init),
+	             t2= identity(op, init), t3= init;
+    difference_type  size= last - first, bsize= size >> 2 << 2, i;
+
+    for (i= 0; i < bsize; i+= 4) {
+	t0= op(t0, Value(first[i]));
+	t1= op(t1, Value(first[i+1]));
+	t2= op(t2, Value(first[i+2]));
+	t3= op(t3, Value(first[i+3]));
+    }
+    for (; i < size; i++)
+	t0= op(t0, Value(first[i]));
+
+    t0= op(t0, t1), t2= op(t2, t3), t0= op(t0, t2);
+    return t0;
+}
+
+} // namespace math
+
+#endif // MATH_ACCUMULATE_INCLUDE

AMDiS/lib/mtl4/boost/numeric/linear_algebra/algebraic_concepts.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 LA_ALGEBRAIC_CONCEPTS_DOC_INCLUDE
+#define LA_ALGEBRAIC_CONCEPTS_DOC_INCLUDE
+
+#ifdef __GXX_CONCEPTS__
+#  include <concepts>
+#else
+#  include <boost/numeric/linear_algebra/pseudo_concept.hpp>
+#endif
+
+#include <boost/numeric/linear_algebra/identity.hpp>
+#include <boost/numeric/linear_algebra/inverse.hpp>
+
+/// Namespace for purely algebraic concepts
+namespace algebra {
+
+/** @addtogroup Concepts
+ *  @{
+ */
+
+#ifndef __GXX_CONCEPTS__
+    //! Concept Commutative
+    /*!
+        \param Operation A functor implementing a binary operation
+        \param Element The type upon the binary operation is defined
+
+        \par Notation:
+        <table summary="notation">
+          <tr>
+            <td>op</td>
+    	<td>Object of type Operation</td>
+          </tr>
+          <tr>
+            <td>x, y</td>
+    	<td>Objects of type Element</td>
+          </tr>
+        </table>
+        \invariant
+        <table summary="invariants">
+          <tr>
+            <td>Commutativity</td>
+    	<td>op(x, y) == op(y, x)</td>
+          </tr>
+        </table>
+     */
+    template <typename Operation, typename Element>
+    struct Commutative
+    {};
+#else
+    concept Commutative<typename Operation, typename Element>
+    {
+	axiom Commutativity(Operation op, Element x, Element y)
+	{
+	    op(x, y) == op(y, x);
+	}
+    };
+#endif
+
+
+#ifndef __GXX_CONCEPTS__
+    //! Concept Associative
+    /*!
+        \param Operation A functor implementing a binary operation
+        \param Element The type upon the binary operation is defined
+
+        \par Notation:
+        <table summary="notation">
+          <tr>
+            <td>op</td>
+    	<td>Object of type Operation</td>
+          </tr>
+          <tr>
+            <td>x, y, z</td>
+    	<td>Objects of type Element</td>
+          </tr>
+        </table>
+        \invariant
+        <table summary="invariants">
+          <tr>
+            <td>Associativity</td>
+    	<td>op(x, op(y, z)) == op(op(x, y), z)</td>
+          </tr>
+        </table>
+     */
+    template <typename Operation, typename Element>
+    struct Associative
+    {};
+#else
+    concept Associative<typename Operation, typename Element>
+    {
+        axiom Associativity(Operation op, Element x, Element y, Element z)
+        {
+	    op(x, op(y, z)) == op(op(x, y), z);
+        }
+    };
+#endif
+
+
+#ifndef __GXX_CONCEPTS__
+    //! Concept SemiGroup
+    /*!
+        \param Operation A functor implementing a binary operation
+        \param Element The type upon the binary operation is defined
+
+        \note
+        -# The algebraic concept SemiGroup only requires associativity and is identical with the concept Associative.
+     */
+    template <typename Operation, typename Element>
+    struct SemiGroup
+        : Associative<Operation, Element>
+    {};
+#else
+    auto concept SemiGroup<typename Operation, typename Element>
+      : Associative<Operation, Element>
+    {};
+#endif
+
+
+#ifndef __GXX_CONCEPTS__
+    //! Concept Monoid
+    /*!
+        \param Operation A functor implementing a binary operation
+        \param Element The type upon the binary operation is defined
+
+        \par Refinement of:
+	- SemiGroup
+        \par Notation:
+        <table summary="notation">
+          <tr>
+            <td>op</td>
+    	<td>Object of type Operation</td>
+          </tr>
+          <tr>
+            <td>x</td>
+    	<td>Object of type Element</td>
+          </tr>
+        </table>
+        \invariant
+        <table summary="invariants">
+          <tr>
+            <td>Neutrality from right</td>
+	    <td>op( x, identity(op, x) ) == x</td>
+          </tr>
+          <tr>
+            <td>Neutrality from left</td>
+	    <td>op( identity(op, x), x ) == x</td>
+          </tr>
+        </table>
+     */
+    template <typename Operation, typename Element>
+    struct Monoid
+      : SemiGroup<Operation, Element>
+    {
+	/// Associated type; if not defined in concept_map automatically detected as result of identity
+        typedef associated_type identity_result_type;
+        identity_result_type identity(Operation, Element); ///< Identity element of Operation
+    };
+#else
+    concept Monoid<typename Operation, typename Element>
+      : SemiGroup<Operation, Element>
+    {
+        typename identity_result_type;
+        identity_result_type identity(Operation, Element);
+
+        axiom Neutrality(Operation op, Element x)
+        {
+	    op( x, identity(op, x) ) == x;
+	    op( identity(op, x), x ) == x;
+        }
+    };
+#endif
+
+#ifdef __GXX_CONCEPTS__
+    auto concept Inversion<typename Operation, typename Element>
+    {
+        typename inverse_result_type;
+        inverse_result_type inverse(Operation, Element);
+
+    };
+#else
+    //! Concept Inversion
+    /*!
+        \param Operation A functor implementing a binary operation
+        \param Element The type upon the binary operation is defined
+
+	\par Associated Types:
+	- inverse_result_type
+	\par Valid Expressions:
+	- inverse(op, x);
+     */
+    template <typename Operation, typename Element>
+    struct Inversion
+    {
+	/// Associated type; if not defined in concept_map automatically detected as result of inverse
+        typedef associated_type inverse_result_type;
+
+	/// Returns inverse of \p x regarding operation \p op
+        inverse_result_type inverse(Operation op, Element x);
+    };
+#endif
+
+
+#ifdef __GXX_CONCEPTS__
+    concept Group<typename Operation, typename Element>
+      : Monoid<Operation, Element>, Inversion<Operation, Element>
+    {
+        axiom Inversion(Operation op, Element x)
+        {
+	    op( x, inverse(op, x) ) == identity(op, x);
+	    op( inverse(op, x), x ) == identity(op, x);
+        }
+    };
+#else
+    //! Concept Group
+    /*!
+        \param Operation A functor implementing a binary operation
+        \param Element The type upon the binary operation is defined
+
+        \par Refinement of:
+	- Monoid
+	- Inversion
+        \par Notation:
+        <table summary="notation">
+          <tr>
+            <td>op</td>
+    	<td>Object of type Operation</td>
+          </tr>
+          <tr>
+            <td>x</td>
+    	<td>Object of type Element</td>
+          </tr>
+        </table>
+        \invariant
+        <table summary="invariants">
+          <tr>
+            <td>Inverse from right</td>
+    	<td>op( x, inverse(op, x) ) == identity(op, x)</td>
+          </tr>
+          <tr>
+            <td>Inverse from left</td>
+    	<td>op( inverse(op, x), x ) == identity(op, x)</td>
+          </tr>
+        </table>
+     */
+    template <typename Operation, typename Element>
+    struct Group
+      : Monoid<Operation, Element>,
+	Inversion<Operation, Element>
+    {};
+#endif
+
+
+#ifdef __GXX_CONCEPTS__
+    auto concept AbelianGroup<typename Operation, typename Element>
+      : Group<Operation, Element>, Commutative<Operation, Element>
+    {};
+#else
+    //! Concept AbelianGroup
+    /*!
+        \param Operation A functor implementing a binary operation
+        \param Element The type upon the binary operation is defined
+
+        \par Refinement of:
+	- Group
+	- Commutative
+     */
+    template <typename Operation, typename Element>
+    struct AbelianGroup
+      : Group<Operation, Element>,
+	Commutative<Operation, Element>
+    {};
+#endif
+
+
+#ifdef __GXX_CONCEPTS__
+    concept Distributive<typename AddOp, typename MultOp, typename Element>
+    {
+        axiom Distributivity(AddOp add, MultOp mult, Element x, Element y, Element z)
+        {
+	    // From left
+	    mult(x, add(y, z)) == add(mult(x, y), mult(x, z));
+	    // z right
+	    mult(add(x, y), z) == add(mult(x, z), mult(y, z));
+        }
+    };
+#else
+    //! Concept Distributive
+    /*!
+        \param AddOp A functor implementing a binary operation representing addition
+        \param MultOp A functor implementing a binary operation representing multiplication
+        \param Element The type upon the binary operation is defined
+
+        \par Notation:
+        <table summary="notation">
+          <tr>
+            <td>add</td>
+	    <td>Object of type AddOp</td>
+          </tr>
+          <tr>
+            <td>mult</td>
+	    <td>Object of type Multop</td>
+          </tr>
+          <tr>
+            <td>x, y, z</td>
+	    <td>Objects of type Element</td>
+          </tr>
+        </table>
+        \invariant
+        <table summary="invariants">
+          <tr>
+            <td>Distributivity from left</td>
+	    <td>mult(x, add(y, z)) == add(mult(x, y), mult(x, z))</td>
+          </tr>
+          <tr>
+            <td>Distributivity from right</td>
+	    <td>mult(add(x, y), z) == add(mult(x, z), mult(y, z))</td>
+          </tr>
+        </table>
+     */
+    template <typename AddOp, typename MultOp, typename Element>
+    struct Distributive
+    {};
+#endif
+
+
+#ifdef __GXX_CONCEPTS__
+    auto concept Ring<typename AddOp, typename MultOp, typename Element>
+      : AbelianGroup<AddOp, Element>,
+        SemiGroup<MultOp, Element>,
+        Distributive<AddOp, MultOp, Element>
+    {};
+#else
+    //! Concept Ring
+    /*!
+        \param AddOp A functor implementing a binary operation representing addition
+        \param MultOp A functor implementing a binary operation representing multiplication
+        \param Element The type upon the binary operation is defined
+
+        \par Refinement of:
+	- AbelianGroup <MultOp, Element>
+	- SemiGroup <MultOp, Element>
+        - Distributive <AddOp, MultOp, Element>
+     */
+    template <typename AddOp, typename MultOp, typename Element>
+    struct Ring
+      : AbelianGroup<AddOp, Element>,
+        SemiGroup<MultOp, Element>,
+        Distributive<AddOp, MultOp, Element>
+    {};
+#endif
+
+
+#ifdef __GXX_CONCEPTS__
+    auto concept RingWithIdentity<typename AddOp, typename MultOp, typename Element>
+      : Ring<AddOp, MultOp, Element>,
+        Monoid<MultOp, Element>
+    {};
+#else
+    //! Concept RingWithIdentity
+    /*!
+        \param AddOp A functor implementing a binary operation representing addition
+        \param MultOp A functor implementing a binary operation representing multiplication
+        \param Element The type upon the binary operation is defined
+
+        \par Refinement of:
+	- Ring <AddOp, MultOp, Element>
+        - Monoid <MultOp, Element>
+     */
+    template <typename AddOp, typename MultOp, typename Element>
+    struct RingWithIdentity
+      : Ring<AddOp, MultOp, Element>,
+        Monoid<MultOp, Element>
+    {};
+#endif
+
+
+#ifdef __GXX_CONCEPTS__
+    concept DivisionRing<typename AddOp, typename MultOp, typename Element>
+      : RingWithIdentity<AddOp, MultOp, Element>,
+        Inversion<MultOp, Element>
+    {
+        // 0 != 1, otherwise trivial
+        axiom ZeroIsDifferentFromOne(AddOp add, MultOp mult, Element x)
+        {
+	    identity(add, x) != identity(mult, x);
+        }
+
+        // Non-zero divisibility from left and from right
+        axiom NonZeroDivisibility(AddOp add, MultOp mult, Element x)
+        {
+	    if (x != identity(add, x))
+		mult(inverse(mult, x), x) == identity(mult, x);
+	    if (x != identity(add, x))
+		mult(x, inverse(mult, x)) == identity(mult, x);
+        }
+    };
+#else
+    //! Concept DivisionRing
+    /*!
+        \param AddOp A functor implementing a binary operation representing addition
+        \param MultOp A functor implementing a binary operation representing multiplication
+        \param Element The type upon the binary operation is defined
+
+        \par Refinement of:
+	- RingWithIdentity <AddOp, MultOp, Element>
+        - Inversion <MultOp, Element>
+
+        \par Notation:
+        <table summary="notation">
+          <tr>
+            <td>add</td>
+	    <td>Object of type AddOp</td>
+          </tr>
+          <tr>
+            <td>mult</td>
+	    <td>Object of type Multop</td>
+          </tr>
+          <tr>
+            <td>x, y, z</td>
+	    <td>Objects of type Element</td>
+          </tr>
+        </table>
+
+        \invariant
+        <table summary="invariants">
+          <tr>
+            <td>Non-zero divisibility from left</td>
+	    <td>mult(inverse(mult, x), x) == identity(mult, x)</td>
+	    <td>if x != identity(add, x)</td>
+          </tr>
+          <tr>
+            <td>Non-zero divisibility from right</td>
+	    <td>mult(x, inverse(mult, x)) == identity(mult, x)</td>
+	    <td>if x != identity(add, x)</td>
+          </tr>
+          <tr>
+            <td>Zero is different from one</td>
+	    <td>identity(add, x) != identity(mult, x)</td>
+	    <td></td>
+          </tr>
+	</table>
+
+	\note
+	-# Zero and one can be theoretically identical in a DivisionRing.  However,
+	   this implies that there is only one element x in the Ring with x + x = x and
+	   x * x = x (which is actually even a Field).
+	   Because this structure has no practical value we exclude it from
+	   consideration.
+     */
+    template <typename AddOp, typename MultOp, typename Element>
+    struct DivisionRing
+      : RingWithIdentity<AddOp, MultOp, Element>,
+        Inversion<MultOp, Element>
+    {};
+#endif
+
+
+#ifdef __GXX_CONCEPTS__
+    // SkewField is defined as synonym for DivisionRing
+    auto concept SkewField<typename AddOp, typename MultOp, typename Element>
+      : DivisionRing<AddOp, MultOp, Element>
+    {};
+#else
+    //! Concept SkewField
+    /*!
+        \param AddOp A functor implementing a binary operation representing addition
+        \param MultOp A functor implementing a binary operation representing multiplication
+        \param Element The type upon the binary operation is defined
+
+        \par Refinement of:
+	- DivisionRing <AddOp, MultOp, Element>
+	\note
+	- Because the refinement of DivisionRing to SkewField is automatic the two concepts
+	  are identical.
+     */
+    template <typename AddOp, typename MultOp, typename Element>
+    struct SkewField
+      : DivisionRing<AddOp, MultOp, Element>
+    {};
+#endif
+
+
+#ifdef __GXX_CONCEPTS__
+    auto concept Field<typename AddOp, typename MultOp, typename Element>
+      : DivisionRing<AddOp, MultOp, Element>,
+        Commutative<MultOp, Element>
+    {};
+#else
+    //! Concept Field
+    /*!
+        \param AddOp A functor implementing a binary operation representing addition
+        \param MultOp A functor implementing a binary operation representing multiplication
+        \param Element The type upon the binary operation is defined
+
+        \par Refinement of:
+	- DivisionRing <AddOp, MultOp, Element>
+	- Commutative <MultOp, Element>
+     */
+    template <typename AddOp, typename MultOp, typename Element>
+    struct Field
+      : DivisionRing<AddOp, MultOp, Element>,
+        Commutative<MultOp, Element>
+    {};
+#endif
+
+/*@}*/ // end of group Concepts
+
+} // algebra
+
+#endif // LA_ALGEBRAIC_CONCEPTS_DOC_INCLUDE

AMDiS/lib/mtl4/boost/numeric/linear_algebra/concept_maps.hpp

+// $COPYRIGHT$
+
+#ifndef MATH_CONCEPT_MAPS_INCLUDE
+#define MATH_CONCEPT_MAPS_INCLUDE
+
+#include <boost/numeric/linear_algebra/intrinsic_concept_maps.hpp>
+#include <boost/numeric/linear_algebra/new_concepts.hpp>
+
+
+namespace math {
+
+#if 0 // Why this doesn't work???
+    template <typename T>
+        requires IntrinsicUnsignedIntegral<T>
+    concept_map UnsignedIntegral<T> {}
+
+    template <typename T>
+        requires IntrinsicSignedIntegral<T>
+    concept_map SignedIntegral<T> {}
+#endif
+
+    concept_map UnsignedIntegral<unsigned int> {}
+    concept_map AdditiveMonoid<unsigned int> {}
+
+    concept_map SignedIntegral<int> {}
+    concept_map AdditiveSemiGroup<int> {}
+
+    concept_map Commutative< add<int>, int > {}
+    concept_map Monoid< add<int>, int > {}
+    concept_map AbelianGroup< add<int>, int > {}
+
+    concept_map Commutative< mult<int>, int > {}
+    concept_map Monoid< mult<int>, int > {}
+    concept_map PIMonoid< mult<int>, int > {}
+
+    concept_map Commutative< min<int>, int > {}
+    concept_map Monoid< min<int>, int > {}
+
+    concept_map Commutative< max<int>, int > {}
+    //concept_map Monoid< max<int>, int > {}
+
+
+    concept_map Commutative< add<float>, float > {}
+    concept_map Monoid< add<float>, float > {}
+    concept_map AbelianGroup< add<float>, float > {}
+
+    concept_map Commutative< mult<float>, float > {}
+    concept_map Monoid< mult<float>, float > {}
+    concept_map PIMonoid< mult<float>, float > {}
+
+    concept_map AdditiveMonoid< float > {}
+    // concept_map AdditiveAbelianGroup< float > {}
+
+    // concept_map OperatorField<float> {}
+
+    concept_map Commutative< min<float>, float > {}
+    concept_map Monoid< min<float>, float > {}
+
+    concept_map Commutative< max<float>, float > {}
+    concept_map Monoid< max<float>, float > {}
+
+
+
+    concept_map Commutative< add<double>, double > {}
+    concept_map Monoid< add<double>, double > {}
+    concept_map AbelianGroup< add<double>, double > {}
+
+    concept_map Commutative< mult<double>, double > {}
+    concept_map Monoid< mult<double>, double > {}
+    concept_map PIMonoid< mult<double>, double > {}
+
+    concept_map AdditiveMonoid< double > {}
+    // concept_map AdditiveAbelianGroup< double > {}
+
+    // concept_map OperatorField<double> {}
+
+    concept_map Commutative< min<double>, double > {}
+    concept_map Monoid< min<double>, double > {}
+
+    concept_map Commutative< max<double>, double > {}
+    concept_map Monoid< max<double>, double > {}
+
+
+
+} // namespace math
+
+#endif // MATH_CONCEPT_MAPS_INCLUDE

AMDiS/lib/mtl4/boost/numeric/linear_algebra/concepts.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_NEW_ALGEBRAIC_CONCEPTS_INCLUDE
+#define MTL_NEW_ALGEBRAIC_CONCEPTS_INCLUDE
+
+#include <concepts>
+#include <boost/numeric/linear_algebra/intrinsic_concept_maps.hpp>
+#include <boost/numeric/linear_algebra/operators.hpp>
+
+
+namespace math {
+
+concept Commutative<typename Operation, typename Element>
+  : std::Callable2<Op, Element, Element>
+{
+    axiom Commutativity(Operation op, Element x, Element y)
+    {
+	op(x, y) == op(y, x);
+    }
+};
+
+
+concept SemiGroup<typename Operation, typename Element>
+  : std::Callable2<Op, Element, Element>
+{
+    axiom Associativity(Operation op, Element x, Element y, Element z)
+    {
+	op(x, op(y, z)) == op(op(x, y), z);
+    }
+};
+
+
+concept Monoid<typename Operation, typename Element>
+  : SemiGroup<Operation, Element>
+{
+    typename identity_result_type;
+    identity_result_type identity(Operation, Element);
+
+    axiom Neutrality(Operation op, Element x)
+    {
+	op( x, identity(op, x) ) == x;
+	op( identity(op, x), x ) == x;
+    }
+};
+
+
+auto concept Inversion<typename Operation, typename Element>
+{
+    typename result_type;
+    result_type inverse(Operation, Element);
+
+};
+
+
+concept PIMonoid<typename Operation, typename Element>
+  : Monoid<Operation, Element>,
+    Inversion<Operation, Element>
+{
+    bool is_invertible(Operation, Element);
+
+    requires std::Convertible<Inversion<Operation, Element>::result_type, Element>;
+
+    axiom Invertibility(Operation op, Element x)
+    {
+	// Only for invertible elements:
+	if (is_invertible(op, x))
+	    op( x, inverse(op, x) ) == identity(op, x);
+	if ( is_invertible(op, x) )
+	    op( inverse(op, x), x ) == identity(op, x);
+    }
+}
+
+#if 0
+    // Alternative approach to convert the result of inversion to Element
+    // Unfortunately, this doesn't compile
+    template <typename Operation, typename Element>
+        requires PIMonoid<Operation, Element>
+    concept_map PIMonoid<Operation, Inversion<Operation, Element>::result_type> {}
+#endif
+
+
+concept Group<typename Operation, typename Element>
+  : PIMonoid<Operation, Element>
+{
+    bool is_invertible(Operation, Element) { return true; }
+
+    // Just in case somebody redefines is_invertible
+    axiom AlwaysInvertible(Operation op, Element x)
+    {
+	is_invertible(op, x);
+    }
+
+    // In fact this is implied by AlwaysInvertible and inherited Invertibility axiom
+    // Maybe remove
+    axiom GlobalInvertibility(Operation op, Element x)
+    {
+	op( x, inverse(op, x) ) == identity(op, x);
+	op( inverse(op, x), x ) == identity(op, x);
+    }
+};
+
+
+auto concept AbelianGroup<typename Operation, typename Element>
+  : Group<Operation, Element>, Commutative<Operation, Element>
+{};
+
+
+// =======================
+// Operator-based concepts
+// =======================
+
+
+concept Additive<typename Element>
+  : std::HasPlus<Element>
+{
+    typename plus_assign_result_type;
+    plus_assign_result_type operator+=(Element& x, Element y)
+    {
+	x= x + y; return x;
+    }
+
+    requires std::Convertible<plus_assign_result_type, Element&>;
+
+    // Do we need the opposite conversion too?
+    // This line produces a compiler error
+    // requires std::Convertible<add<Element>::result_type,
+    //                           std::HasPlus<Element>::result_type>;
+
+    axiom Consistency(add<Element> op, Element x, Element y)
+    {
+	op(x, y) == x + y;
+	op(x, y) == (x += y, x);
+    }
+}
+
+
+auto concept AdditiveCommutative<typename Element>
+  : Additive<Element>,
+    Commutative< add<Element>, Element >
+{}
+
+
+auto concept AdditiveSemiGroup<typename Element>
+  : Additive<Element>,
+    SemiGroup< add<Element>, Element >
+{}
+
+
+concept AdditiveMonoid<typename Element>
+  : AdditiveSemiGroup<Element>,
+    Monoid< add<Element>, Element >
+{
+    Element zero(Element x)
+    {
+	return identity(add<Element>(), x);
+    }
+
+    // If we don't use the default definition
+    axiom IdentityConsistency (add<Element> op, Element x)
+    {
+	zero(x) == identity(op, x);
+    }
+};
+
+#ifndef CONCEPTS_WITHOUT_OVERLOADED_REQUIREMENTS  // Uncompilable due to error in compiler
+concept AdditivePIMonoid<typename Element>
+  : std::HasMinus<Element>, AdditiveMonoid<Element>,
+    PIMonoid< add<Element>, Element >
+{
+    typename minus_assign_result_type;
+    minus_assign_result_type operator-=(Element& x, Element y)
+    {
+	x= x - y; return x;
+    }
+
+    requires std::Convertible<minus_assign_result_type, Element&>;
+
+    typename unary_result_type;
+    unary_result_type operator-(Element x)
+    {
+	return zero(x) - x;
+    }
+
+    axiom InverseConsistency(add<Element> op, Element x, Element y)
+    {
+	// consistency between additive and functor concept
+	if ( is_invertible(op, y) )
+	    op(x, inverse(op, y)) == x - y;
+	if ( is_invertible(op, y) )
+	    op(x, y) == (x -= y, x);
+
+	// consistency of unary inversion
+	if ( is_invertible(op, y) )
+	    inverse(op, y) == -y;
+
+	// consistency between unary and binary -
+	if ( is_invertible(op, x) )
+	    identity(op, x) - x == -x;
+    }
+}
+
+
+auto concept AdditiveGroup<typename Element>
+  : AdditivePIMonoid<Element>,
+    Group< add<Element>, Element >
+{};
+
+
+auto concept AdditiveAbelianGroup<typename Element>
+  : AdditiveGroup<Element>,
+    Commutative< add<Element>, Element >
+{}
+
+#endif
+
+
+concept Multiplicative<typename Element>
+  : std::HasMultiply<Element>
+{
+    typename times_assign_result_type;
+    times_assign_result_type operator*=(Element& x, Element y)
+    {
+	x= x * y; return x;
+    }
+
+    requires std::Convertible<times_assign_result_type, Element&>;
+
+    // Do we need the opposite conversion too?
+    // This line produces a compiler error
+    // requires std::Convertible<mult<Element>::result_type,
+    //                           std::HasMultiply<Element>::result_type>;
+
+    axiom Consistency(mult<Element> op, Element x, Element y)
+    {
+	op(x, y) == x * y;
+	op(x, y) == (x *= y, x);
+    }
+}
+
+
+auto concept MultiplicativeCommutative<typename Element>
+  : Multiplicative<Element>,
+    Commutative< mult<Element>, Element >
+{}
+
+
+auto concept MultiplicativeSemiGroup<typename Element>
+  : Multiplicative<Element>,
+    SemiGroup< mult<Element>, Element >
+{}
+
+
+concept MultiplicativeMonoid<typename Element>
+  : MultiplicativeSemiGroup<Element>,
+    Monoid< mult<Element>, Element >
+{
+    Element one(Element x)
+    {
+	return identity(mult<Element>(), x);
+    }
+
+    // If we don't use the default definition
+    axiom IdentityConsistency (math::mult<Element> op, Element x)
+    {
+	one(x) == identity(op, x);
+    }
+};
+
+#ifndef CONCEPTS_WITHOUT_OVERLOADED_REQUIREMENTS  // Uncompilable due to error in compiler
+concept MultiplicativePIMonoid<typename Element>
+  : std::HasDivide<Element>, MultiplicativeMonoid<Element>,
+    PIMonoid< mult<Element>, Element >
+{
+    typename divide_assign_result_type;
+    divide_assign_result_type operator/=(Element& x, Element y)
+    {
+	x= x / y; return x;
+    }
+
+    requires std::Convertible<divide_assign_result_type, Element&>;
+
+    axiom InverseConsistency(mult<Element> op, Element x, Element y)
+    {
+	// consistency between multiplicative and functor concept
+	if ( is_invertible(op, y) )
+	    op(x, inverse(op, y)) == x / y;
+	if ( is_invertible(op, y) )
+	    op(x, y) == (x /= y, x);
+    }
+}
+
+
+auto concept MultiplicativeGroup<typename Element>
+  : MultiplicativePIMonoid<Element>,
+    Group< mult<Element>, Element >
+{};
+
+
+auto concept MultiplicativeAbelianGroup<typename Element>
+  : MultiplicativeGroup<Element>,
+    Commutative< mult<Element>, Element >
+{}
+
+
+// ==========================
+// Concepts with 2 operations
+// ==========================
+
+
+
+concept Distributive<typename AddOp, typename MultOp, typename Element>
+{
+    axiom Distributivity(AddOp add, MultOp mult, Element x, Element y, Element z)
+    {
+	// From left
+	mult(x, add(y, z)) == add(mult(x, y), mult(x, z));
+	// from right
+	mult(add(x, y), z) == add(mult(x, z), mult(y, z));
+    }
+}
+
+
+auto concept Ring<typename AddOp, typename MultOp, typename Element>
+  : AbelianGroup<AddOp, Element>,
+    SemiGroup<MultOp, Element>,
+    Distributive<AddOp, MultOp, Element>
+{}
+
+
+auto concept RingWithIdentity<typename AddOp, typename MultOp, typename Element>
+  : Ring<AddOp, MultOp, Element>,
+    Monoid<MultOp, Element>
+{}
+
+
+concept DivisionRing<typename AddOp, typename MultOp, typename Element>
+  : RingWithIdentity<AddOp, MultOp, Element>,
+    Inversion<MultOp, Element>
+{
+    // 0 != 1, otherwise trivial
+    axiom ZeroIsDifferentFromOne(AddOp add, MultOp mult, Element x)
+    {
+	identity(add, x) != identity(mult, x);
+    }
+
+    // Non-zero divisibility from left and from right
+    axiom NonZeroDivisibility(AddOp add, MultOp mult, Element x)
+    {
+	if (x != identity(add, x))
+	    mult(inverse(mult, x), x) == identity(mult, x);
+	if (x != identity(add, x))
+	    mult(x, inverse(mult, x)) == identity(mult, x);
+    }
+}
+
+
+auto concept Field<typename AddOp, typename MultOp, typename Element>
+  : DivisionRing<AddOp, MultOp, Element>,
+    Commutative<MultOp, Element>
+{}
+
+
+auto concept OperatorRing<typename Element>
+  : AdditiveAbelianGroup<Element>,
+    MultiplicativeSemiGroup<Element>,
+    Ring<add<Element>, mult<Element>, Element>
+{}
+
+
+auto concept OperatorRingWithIdentity<typename Element>
+  : OperatorRing<Element>,
+    MultiplicativeMonoid<Element>,
+    RingWithIdentity<add<Element>, mult<Element>, Element>
+{}
+
+
+auto concept OperatorDivisionRing<typename Element>
+  : OperatorRingWithIdentity<Element>,
+    MultiplicativePIMonoid<Element>,
+    DivisionRing<add<Element>, mult<Element>, Element>
+{}
+
+
+auto concept OperatorField<typename Element>
+  : OperatorDivisionRing<Element>,
+    Field<add<Element>, mult<Element>, Element>
+{}
+
+
+
+#endif
+
+concept IntrinsicType<typename T> {}
+
+concept IntrinsicArithmetic<typename T> : IntrinsicType<T> {}
+
+concept IntrinsicIntegral<typename T> : IntrinsicArithmetic<T> {}
+
+concept IntrinsicSignedIntegral<typename T>
+  : std::SignedIntegralLike<T>,
+    IntrinsicIntegral<T>
+{}
+
+concept IntrinsicUnsignedIntegral<typename T>
+  : std::UnsignedIntegralLike<T>,
+    IntrinsicIntegral<T>
+{}
+
+concept IntrinsicFloatingPoint<typename T>
+  : std::FloatingPointLike<T>,
+    IntrinsicArithmetic<T>
+{}
+
+
+// Intrinsic types are chategorized:
+
+concept_map IntrinsicSignedIntegral<char> {}
+concept_map IntrinsicSignedIntegral<signed char> {}
+concept_map IntrinsicUnsignedIntegral<unsigned char> {}
+concept_map IntrinsicSignedIntegral<short> {}
+concept_map IntrinsicUnsignedIntegral<unsigned short> {}
+concept_map IntrinsicSignedIntegral<int> {}
+concept_map IntrinsicUnsignedIntegral<unsigned int> {}
+concept_map IntrinsicSignedIntegral<long> {}
+concept_map IntrinsicUnsignedIntegral<unsigned long> {}
+concept_map IntrinsicSignedIntegral<long long> {}
+concept_map IntrinsicUnsignedIntegral<unsigned long long> {}
+
+concept_map IntrinsicFloatingPoint<float> {}
+concept_map IntrinsicFloatingPoint<double> {}
+
+
+#ifndef CONCEPTS_WITHOUT_OVERLOADED_REQUIREMENTS
+
+// ====================
+// Default Concept Maps
+// ====================
+
+// ==============
+// Arithmetic
+// ==============
+
+// ----------------
+// Signed integrals
+// ----------------
+
+template <typename T>
+  requires IntrinsicSignedIntegral<T>
+concept_map OperatorRingWithIdentity<T> {}
+
+template <typename T>
+  requires IntrinsicSignedIntegral<T>
+concept_map MultiplicativeCommutative<T> {}
+
+// ------------------
+// Unsigned integrals
+// ------------------
+
+
+template <typename T>
+  requires IntrinsicUnsignedIntegral<T>
+concept_map AdditiveCommutative<T> {}
+
+template <typename T>
+  requires IntrinsicUnsignedIntegral<T>
+concept_map AdditiveMonoid<T> {}
+
+template <typename T>
+  requires IntrinsicUnsignedIntegral<T>
+concept_map MultiplicativeCommutative<T> {}
+
+template <typename T>
+  requires IntrinsicUnsignedIntegral<T>
+concept_map MultiplicativeMonoid<T> {}
+
+// ---------------
+// Floationg Point
+// ---------------
+
+
+template <typename T>
+  requires IntrinsicFloatingPoint<T>
+concept_map Field<T> {}
+
+template <typename T>
+  requires IntrinsicFloatingPoint<T>
+concept_map Field< std::complex<T> > {}
+
+
+// ===========
+// Min and Max
+// ===========
+
+
+template <typename T>
+  requires IntrinsicArithmetic<T>
+concept_map Commutative< max<T>, T > {}
+
+template <typename T>
+  requires IntrinsicArithmetic<T>
+concept_map Monoid< max<T>, T > {}
+
+template <typename T>
+  requires IntrinsicArithmetic<T>
+concept_map Commutative< min<T>, T > {}
+
+template <typename T>
+  requires IntrinsicArithmetic<T>
+concept_map Monoid< min<T>, T > {}
+
+
+// ==========
+// And and Or
+// ==========
+
+template <typename T>
+  requires Intrinsic<T> && std::HasLogicalAnd<T>
+concept_map Commutative< std::logical_and<T>, T > {}
+
+template <typename T>
+  requires Intrinsic<T> && std::HasLogicalAnd<T>
+concept_map Monoid< std::logical_and<T>, T > {}
+
+template <typename T>
+  requires Intrinsic<T> && std::HasLogicalOr<T>
+concept_map Commutative< std::logical_or<T>, T > {}
+
+template <typename T>
+  requires Intrinsic<T> && std::HasLogicalOr<T>
+concept_map Monoid< std::logical_or<T>, T > {}
+
+template <typename T>
+  requires Intrinsic<T> && std::HasLogicalAnd<T> && std::HasLogicalOr<T>
+concept_map Distributive<std::logical_and<T>, std::logical_or<T>, T> {}
+
+template <typename T>
+  requires Intrinsic<T> && std::HasLogicalAnd<T> && std::HasLogicalOr<T>
+concept_map Distributive<std::logical_or<T>, std::logical_and<T>, T> {}
+
+
+// ==================
+// Bitwise operations
+// ==================
+
+// not yet defined
+
+template <typename T>
+  requires IntrinsicIntegral<T>
+concept_map Commutative< bit_and<T>, T > {}
+
+template <typename T>
+  requires IntrinsicIntegral<T>
+concept_map Monoid< bit_and<T>, T > {}
+
+template <typename T>
+  requires IntrinsicIntegral<T>
+concept_map Commutative< bit_or<T>, T > {}
+
+template <typename T>
+  requires IntrinsicIntegral<T>
+concept_map Monoid< bit_or<T>, T > {}
+
+template <typename T>
+  requires IntrinsicIntegral<T>
+concept_map Distributive<bit_and<T>, bit_or<T>, T> {}
+
+template <typename T>
+  requires IntrinsicIntegral<T>
+concept_map Distributive<bit_or<T>, bit_and<T>, T> {}
+
+template <typename T>
+  requires IntrinsicIntegral<T>
+concept_map Commutative< bit_xor<T>, T > {}
+
+template <typename T>
+  requires IntrinsicIntegral<T>
+concept_map SemiGroup< bit_xor<T>, T > {}
+
+// ====================
+// String concatenation
+// ====================
+
+concept_map AdditiveMonoid<std::string> {}
+
+
+#endif
+
+
+
+} // namespace math
+
+#endif // MTL_NEW_ALGEBRAIC_CONCEPTS_INCLUDE

AMDiS/lib/mtl4/boost/numeric/linear_algebra/ets_concepts.hpp

+// Copyright 2006. Peter Gottschling, Matthias Troyer, Rolf Bonderer
+
+#ifndef LA_ETS_CONCEPTS_INCLUDE
+#define LA_ETS_CONCEPTS_INCLUDE
+
+#include <boost/numeric/linear_algebra/concepts.hpp>
+
+#ifdef  __GXX_CONCEPTS__
+
+//"Namespace ETS = Expression Templates Support":
+// This file contains concepts that support the concepts defined in the namespace math.
+// They are required since the math-concepts do not support ExpressionTemplates so far.
+// The list of valid expressions is in fact infinite, so we just name some of them.
+// Once ExpressionTemplates are supported by the math-concepts, the ets-concepts are no longer required.
+
+namespace ets {
+
+  auto concept Field<typename Element>
+  {
+    requires std::Assignable<Element, math::AdditivePartiallyInvertibleMonoid<Element>::unary_result_type>; //"x=-y" valid
+  };
+
+  auto concept VectorSpace<typename Vector, typename Scalar>
+  {
+    // valid expression: "vector2 += scalar*vector1"
+    typename res_type_1;
+    res_type_1 operator+=(Vector&, math::Multiplicable<Scalar, Vector>::result_type);
+
+    // valid expression: "vector2 -= scalar*vector1"
+    typename res_type_2;
+    res_type_2 operator-=(Vector&, math::Multiplicable<Scalar, Vector>::result_type);
+
+    //valid expression: "vector *= -scalar"
+    typename res_type_3;
+    res_type_3 operator*=(Vector&, const math::AdditivePartiallyInvertibleMonoid<Scalar>::unary_result_type&);
+
+    //valid expression: "vector3 = vector1 + scalar*vector2"
+    requires math::Addable<Vector, math::Multiplicable<Scalar, Vector>::result_type>; //"vector1+scalar*vector2" valid
+    requires std::Assignable<Vector, math::Addable<Vector, math::Multiplicable<Scalar, Vector>::result_type>::result_type>; //"vector3 = vector1 + scalar*vector2" valid
+
+  };
+
+} //namespace ets
+
+#endif //__GXX_CONCEPTS__
+
+#endif //LA_ETS_CONCEPTS_INCLUDE