AMDiS/lib/mtl4/boost/numeric/itl/pc/solver.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_PC_SOLVER_INCLUDE
+#define ITL_PC_SOLVER_INCLUDE
+
+#include <boost/mpl/bool.hpp>
+#include <boost/numeric/mtl/vector/assigner.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+namespace itl { namespace pc {
+
+/// Helper class for delayed (i.e. copy-free) evaluation of preconditioners
+template <typename PC, typename Vector, bool adjoint= false>
+struct solver
+  : mtl::vector::assigner<solver<PC, Vector> >
+{
+    typedef PC  pc_type;
+
+    /// Constructor taking preconditioner and source vector
+    solver(const pc_type& P, const Vector& x) : P(P), x(x) {}
+
+    /// Assign result to vector \p y, if possible without copying
+    template <typename VectorOut>
+    void assign_to(VectorOut& y) const
+    {
+	mtl::vampir_trace<5055> tracer;
+	assign_to(y, boost::mpl::bool_<adjoint>());
+    }
+
+  protected:
+
+    template <typename VectorOut>
+    void assign_to(VectorOut& y, boost::mpl::false_) const
+    {	P.solve(x, y);    }
+
+    template <typename VectorOut>
+    void assign_to(VectorOut& y, boost::mpl::true_) const
+    {	P.adjoint_solve(x, y);    }
+
+    const pc_type&        P;
+    const Vector&         x;
+};
+
+}} // namespace itl::pc
+
+#endif // ITL_PC_SOLVER_INCLUDE

AMDiS/lib/mtl4/boost/numeric/itl/pc/sorting.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.
+//
+// Algorithm inspired by Nick Vannieuwenhoven, written by Cornelius Steinhardt
+
+/*
+ *
+ *  Created on: Nov 3, 2009
+ *      Author: heazk
+ */
+
+#ifndef MTL_SORTING_INCLUDE
+#define MTL_SORTING_INCLUDE
+
+#include <string.h>
+#include <stdlib.h>
+
+#define BITS_PER_BYTE 8
+
+////////////////////////////////////////////////////////////////////////////////
+// RADIX-SORT
+////////////////////////////////////////////////////////////////////////////////
+
+template<
+	class Type
+>
+void radix_sort(
+	Type* values,
+	const int size,
+	Type* tmp = 0
+) {
+	if(size < 2) {
+		return;
+	}
+
+	Type *const orig_ptr = values;
+	if(tmp == 0) {
+		tmp = new Type[size];
+	}
+
+#if VERBOSE_MODE > 10
+	std::cout << "unsorted: ";
+	for(int i = 0; i < size; ++i) {
+		std::cout << values[i] << ", ";
+	}
+	std::cout << std::endl;
+#endif
+
+	// Determine the maximum number of bits needed.
+	int max = INT_MIN;
+	for(int i = 0; i < size; ++i) {
+		max = (values[i] <= max) ? max : values[i];
+	}
+	int maxBits = 1;
+	for(; max; max >>= 1, ++maxBits){};
+	for(max = 0; max < maxBits; max += BITS_PER_BYTE){};
+	maxBits = max;
+
+	const int groupBits = BITS_PER_BYTE;
+	const int intBits   = maxBits;
+	const int groupSize = 1 << groupBits;
+
+	const int nbIterations = intBits / groupBits;
+	const int mask = groupSize-1;
+
+	// Counting and prefix arrays.
+	int* count = new int[groupSize];
+
+	for (
+			int shift = 0, it = 0;
+			it < nbIterations;
+			++it, shift += groupBits
+	) {
+		// Reset count array.
+		memset(count, 0, groupSize*sizeof(int));
+
+		// Counting elements in the it-th group.
+		for (int i = 0; i < size; ++i) {
+			count[ (values[i] >> shift) & mask ]++;
+		}
+
+		// Calculate prefixes.
+		int currIncrement = -1;
+		int prevCount = count[0];
+		count[0] = 0;
+		for (int i = 1; i < groupSize; ++i) {
+			currIncrement = prevCount;
+			prevCount = count[i];
+			count[i] = count[i-1] + currIncrement;
+		}
+
+		// Sort elements in the it-th group.
+		for (int i = 0; i < size; ++i) {
+			int& index = count[( values[i] >> shift) & mask];
+			tmp[ index ] = values[i];
+			++index;
+		}
+
+		// Swap pointers.
+		Type* tmp_values = values;
+		values = tmp;
+		tmp = tmp_values;
+	}
+
+	// Make sure the sorted keys and values reside in the input array.
+	if( values != orig_ptr ) {
+		memcpy( orig_ptr, values, sizeof(Type)*size );
+		tmp = values;
+	}
+
+#if VERBOSE_MODE > 10
+	std::cout << "sorted: ";
+	for(int i = 0; i < size; ++i) {
+		std::cout << values[i] << ", ";
+	}
+	std::cout << std::endl;
+#endif
+
+	if(count) {
+		delete[] count;
+	}
+	if(tmp) {
+		delete[] tmp;
+	}
+}
+
+template<
+	class Key, class Val
+>
+void radix_sort(
+		Key* keys, Val* vals,
+		const int size,
+		Key* tmpKeys = 0, Val* tmpVals = 0
+) {
+	if(size < 2) {
+		return;
+	}
+
+	// The pointers to the original data.
+	Key *const origKey = keys;
+	Val *const origVal = vals;
+
+	// Helper arrays.
+	bool ownKeys = false;
+	bool ownVals = false;
+	if(!tmpKeys) {
+		tmpKeys = new Key[size];
+		ownKeys = true;
+	}
+	if(!tmpVals) {
+		tmpVals = new Val[size];
+		ownVals = true;
+	}
+
+	// Determine the maximum number of bits needed.
+	int max = -(INT_MAX-1);
+	for(int i = 0; i < size; ++i) {
+		max = keys[i] <= max ? max : keys[i];
+	}
+	int maxBits = 1;
+	for(; max; max >>= 1, ++maxBits){};
+	for(max = 0; max < maxBits; max += BITS_PER_BYTE){};
+	maxBits = max;
+
+	const int groupBits = BITS_PER_BYTE;
+	const int intBits   = maxBits;
+	const int groupSize = 1 << groupBits;
+
+	const int nbIterations = intBits / groupBits;
+	const int mask = groupSize-1;
+
+	// Counting and prefix arrays.
+	int* count = new int[groupSize];
+	// Allocation of memory could fail!
+
+	for (
+			int shift = 0, it = 0;
+			it < nbIterations;
+			++it, shift += groupBits
+	) {
+		// Reset count array.
+		for(int i = 0; i < groupSize; ++i) {
+			count[i] = 0;
+		}
+
+		// Counting elements in the it-th group.
+		for (int i = 0; i < size; ++i) {
+			count[ (keys[i] >> shift) & mask ]++;
+		}
+
+		// Calculate prefixes.
+		int currIncrement = -1;
+		int prevCount = count[0];
+		count[0] = 0;
+		for (int i = 1; i < groupSize; ++i) {
+			currIncrement = prevCount;
+			prevCount = count[i];
+			count[i] = count[i-1] + currIncrement;
+		}
+
+		// Sort elements in the it-th group.
+		for (int i = 0; i < size; ++i) {
+			int& index = count[( keys[i] >> shift) & mask];
+			tmpKeys[ index ] = keys[i];
+			tmpVals[ index ] = vals[i];
+			++index;
+		}
+
+		// Swap pointers.
+		Key* tmpKeyPtr = tmpKeys;
+		Val* tmpValuePtr = tmpVals;
+		tmpKeys = keys;
+		tmpVals = vals;
+		keys = tmpKeyPtr;
+		vals = tmpValuePtr;
+	}
+
+	// Make sure the sorted keys and values reside in the input array.
+	if(keys != origKey) {
+		memcpy( origKey, keys, sizeof(Key)*size );
+		tmpKeys = keys;
+	}
+	if(vals != origVal) {
+		memcpy( origVal, vals, sizeof(Val)*size );
+		tmpVals = vals;
+	}
+
+	if(count) 				delete[] count;
+	if(ownKeys && tmpKeys) 	delete[] tmpKeys;
+	if(ownVals && tmpVals)	delete[] tmpVals;
+	tmpKeys = 0;
+	tmpVals = 0;
+	count = 0;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+// QUICK-SORT
+////////////////////////////////////////////////////////////////////////////////
+
+template<typename Key, typename Val>
+inline void swap(Key& key0, Val& val0, Key& key1, Val& val1) {
+	Key tmp_key = key0;
+	Val tmp_val = val0;
+	key0 = key1;
+	val0 = val1;
+	key1 = tmp_key;
+	val1 = tmp_val;
+}
+
+template<typename Key, typename Value>
+inline void sort_along(Key* keys, Value* values, int size) {
+
+
+	if(size < 2) {
+		return;
+	}
+	if( size == 2 ) {
+		if(keys[0] > keys[1]) {
+			swap(
+				keys[0], values[0],
+				keys[1], values[1]
+			);
+		}
+		return;
+	}
+
+	int piv = std::rand() % size;
+
+
+	swap(
+		keys[0],   values[0],
+		keys[piv], values[piv]
+	);
+	int begin = 1;
+	int end = size-1;
+	while(begin < end) {
+		for(; (begin < end) && (keys[begin] <= keys[0]); ++begin){};
+		for(; (begin <= end) && (keys[end] > keys[0]); --end){};
+
+		if(begin < end) {
+			swap(
+				keys[begin], values[begin],
+				keys[end],   values[end]
+			);
+		}
+	}
+	piv = end;
+	swap(
+		keys[0],   values[0],
+		keys[piv], values[piv]
+	);
+
+
+	sort_along( keys, values, piv );
+	sort_along( keys+(piv+1), values+(piv+1), size-piv-1);
+}
+
+#endif // MTL_SORTING_INCLUDE

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