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

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

+// Copyright 2006. Peter Gottschling, Matthias Troyer, Rolf Bonderer
+// 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_IDENTITY_INCLUDE
+#define MATH_IDENTITY_INCLUDE
+
+#include <boost/math/tools/config.hpp>
+#include <limits>
+#include <string>
+#include <functional>
+
+#include <boost/numeric/linear_algebra/operators.hpp>
+
+namespace math {
+
+template <typename Operation, typename Element>
+struct identity_t {};
+
+// TBD: Do we the case that the return type is different? Using std::unary_function?
+
+// Additive identity of Element type is by default a converted 0
+// However, for vectors one needs to know the dimension
+// (and in parallel eventually also the distribution).
+// Therefore, an element is passed as reference.
+// It is strongly recommended to specialize this functor
+// for better efficiency.
+template <typename Element>
+struct identity_t< add<Element>, Element >
+  : public std::binary_function<add<Element>, Element, Element>
+{
+    Element operator() (const add<Element>&, const Element& ref) const
+    {
+	Element tmp(ref);
+	tmp= 0;
+	return tmp;
+    }
+};
+
+template <>
+struct identity_t< add<std::string>, std::string >
+  : public std::binary_function<add<std::string>, std::string, std::string>
+{
+    std::string operator() (const add<std::string>&, const std::string&) const
+    {
+	return std::string();
+    }
+};
+
+// Multiplicative identity of Element type is by default a converted 1
+// Same comments as above.
+// In contrast to additive identity, this default more likely to be wrong (e.g. matrices with all 1s)
+template <typename Element>
+struct identity_t< mult<Element>, Element >
+  : public std::binary_function<mult<Element>, Element, Element>
+{
+    Element operator() (const mult<Element>&, const Element& ref) const
+    {
+	Element tmp(ref);
+	tmp= 1;
+	return tmp;
+    }
+};
+
+
+// Identity of max is minimal representable value, for standard types defined in numeric_limits
+template <typename Element>
+struct identity_t< max<Element>, Element >
+  : public std::binary_function<max<Element>, Element, Element>
+{
+    Element operator() (const max<Element>&, const Element& ) const
+    {
+	using std::numeric_limits;
+	return numeric_limits<Element>::min();
+    }
+};
+
+template <>
+struct identity_t< max<float>, float >
+  : public std::binary_function<max<float>, float, float>
+{
+    float operator() (const max<float>&, const float& ) const
+    {
+	using std::numeric_limits;
+	return -numeric_limits<float>::max();
+    }
+};
+
+template <>
+struct identity_t< max<double>, double >
+  : public std::binary_function<max<double>, double, double>
+{
+    double operator() (const max<double>&, const double& ) const
+    {
+	using std::numeric_limits;
+	return -numeric_limits<double>::max();
+    }
+};
+
+
+#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
+
+template <>
+struct identity_t< max<long double>, long double >
+  : public std::binary_function<max<long double>, long double, long double>
+{
+    long double operator() (const max<long double>&, const long double& ) const
+    {
+	using std::numeric_limits;
+	return -numeric_limits<long double>::max();
+    }
+};
+
+#endif
+
+
+
+// Identity of min is maximal representable value, for standard types defined in numeric_limits
+template <typename Element>
+struct identity_t< min<Element>, Element >
+  : public std::binary_function<min<Element>, Element, Element>
+{
+    Element operator() (const min<Element>&, const Element& ) const
+    {
+	using std::numeric_limits;
+	return numeric_limits<Element>::max();
+    }
+};
+
+// Identity of bit-wise and
+template <typename Element>
+struct identity_t< bitwise_and<Element>, Element >
+  : public std::binary_function<bitwise_and<Element>, Element, Element>
+{
+    Element operator() (const bitwise_and<Element>&, const Element&) const
+    {
+	return 0;
+    }
+};
+
+// Identity of bit-wise or
+template <typename Element>
+struct identity_t< bitwise_or<Element>, Element >
+  : public std::binary_function<bitwise_or<Element>, Element, Element>
+{
+    Element operator() (const bitwise_or<Element>&, const Element&) const
+    {
+	return 0 - 1;
+    }
+};
+
+#if 0 // ambiguous specialization
+template <template <typename> class Operation, typename First, typename Second>
+struct identity_t< Operation<std::pair<First, Second> >, std::pair<First, Second> >
+{
+    typedef std::pair<First, Second> pt;
+
+    pt operator()(const Operation<pt>&, const pt& ref) const
+    {
+	return std::make_pair(identity(Operation<First>(), ref.first), identity(Operation<Second>(), ref.second));
+    }
+};
+#endif
+
+// Function is shorter than typetrait-like functor
+template <typename Operation, typename Element>
+inline Element identity(const Operation& op, const Element& v)
+{
+    return identity_t<Operation, Element>() (op, v);
+}
+
+#if 1
+// I shouldn't do this (but as functor I'd need too many specializations)
+template <template <typename> class Operation, typename First, typename Second>
+inline std::pair<First, Second> identity(const Operation<std::pair<First, Second> >&, const std::pair<First, Second>& v)
+{
+    return std::pair<First, Second>(math::identity(Operation<First>(), v.first), math::identity(Operation<Second>(), v.second));
+}
+#endif
+
+// Short-cut for additive identity
+template <typename Element>
+inline Element zero(const Element& v)
+{
+    return identity_t<math::add<Element>, Element>() (math::add<Element>(), v);
+}
+
+
+// Short-cut for multiplicative identity
+template <typename Element>
+inline Element one(const Element& v)
+{
+    return identity_t<math::mult<Element>, Element>() (math::mult<Element>(), v);
+}
+
+
+} // namespace math
+
+#endif // MATH_IDENTITY_INCLUDE

AMDiS/lib/mtl4/boost/numeric/linear_algebra/intrinsic_concept_maps.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_INTRINSIC_CONCEPT_MAPS_INCLUDE
+#define MATH_INTRINSIC_CONCEPT_MAPS_INCLUDE
+
+#include <concepts>
+
+namespace math {
+
+
+    // The following concepts are used to classify intrinsic arithmetic types.
+    // The standard concepts already define the syntactic requirements,
+    // i.e. the interface.
+    // However they sey nothing about the semantics.
+    // Therefore, user-defined types can model the syntactic/interface
+    // requirements while still having a different mathematical behavior.
+    // For that reason, we introduce concepts that are only used for intrinsic types.
+    // For them we can define concept_maps regarding semantic behavior as monoids.
+
+    concept IntrinsicSignedIntegral<typename T>
+      : std::SignedIntegralLike<T>
+    {}
+
+    concept IntrinsicUnsignedIntegral<typename T>
+      : std::UnsignedIntegralLike<T>
+    {}
+
+    concept IntrinsicFloatingPoint<typename T>
+      : std::FloatingPointLike<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> { }
+
+    concept Integral<typename T> : std::CopyAssignable<T>
+    {
+        requires std::HasPlus<T> && std::HasMinus<T> && std::HasMultiply<T> && std::HasDivide<T>
+              && std::HasUnaryPlus<T> && std::HasNegate<T>;
+
+
+        T& operator++(T&);
+        T operator++(T& t, int) { T tmp(t); ++t; return tmp; }
+        T& operator--(T&);
+        T operator--(T& t, int) { T tmp(t); --t; return tmp; }
+
+        requires std::Convertible<std::HasUnaryPlus<T>::result_type, T>
+              && std::Convertible<std::HasNegate<T>::result_type, T>
+              && std::Convertible<std::HasPlus<T>::result_type, T>
+              && std::Convertible<std::HasMinus<T>::result_type, T>
+              && std::Convertible<std::HasMultiply<T>::result_type, T>
+              && std::Convertible<std::HasDivide<T>::result_type, T>;
+
+        T& operator*=(T&, T);
+        T& operator/=(T&, T);
+        T& operator+=(T&, T);
+        T& operator-=(T&, T);
+
+        requires std::HasComplement<T> && std::HasModulus<T> && std::HasBitAnd<T>
+              && std::HasBitOr<T> && std::HasBitXor<T> && std::HasLeftShift<T>
+              && std::HasRightShift<T>;
+
+        requires std::Convertible<std::HasComplement<T>::result_type, T>
+              && std::Convertible<std::HasModulus<T>::result_type, T>
+              && std::Convertible<std::HasBitAnd<T>::result_type, T>
+              && std::Convertible<std::HasBitOr<T>::result_type, T>
+              && std::Convertible<std::HasBitXor<T>::result_type, T>
+              && std::Convertible<std::HasLeftShift<T>::result_type, T>
+              && std::Convertible<std::HasRightShift<T>::result_type, T>;
+
+        requires std::LessThanComparable<T> && std::EqualityComparable<T>;
+
+        T& operator%=(T&, T);
+        T& operator&=(T&, T);
+        T& operator|=(T&, T);
+        T& operator^=(T&, T);
+        T& operator<<=(T&, T);
+        T& operator>>=(T&, T);
+    }
+
+
+    template <typename T>
+    requires IntrinsicUnsignedIntegral<T>
+    concept_map Integral<T> {typedef T result_type;}
+
+    template <IntrinsicSignedIntegral T>
+    concept_map Integral<T> {typedef T result_type;}
+
+
+} // namespace math
+
+#endif // MATH_INTRINSIC_CONCEPT_MAPS_INCLUDE

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

+// Copyright 2006. Peter Gottschling, Matthias Troyer, Rolf Bonderer
+// 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_INVERSE_INCLUDE
+#define MATH_INVERSE_INCLUDE
+
+#include <boost/numeric/linear_algebra/operators.hpp>
+#include <boost/numeric/linear_algebra/identity.hpp>
+
+namespace math {
+
+template <typename Operation, typename Element>
+struct inverse_t {} ;
+
+
+template <typename Element>
+struct inverse_t< add<Element>, Element >
+  : public std::binary_function<add<Element>, Element, Element>
+{
+    Element operator()(const add<Element>&, const Element& v) const
+    {
+	return -v;
+    }
+};
+
+
+template <typename Element>
+struct inverse_t< mult<Element>, Element >
+  : public std::binary_function<mult<Element>, Element, Element>
+{
+    Element operator()(const mult<Element>&, const Element& v) const
+    {
+	return one(v) / v ;
+    }
+};
+
+
+// Function is shorter than typetrait-like functor
+template <typename Operation, typename Element>
+inline Element inverse(const Operation& op, const Element& v)
+{
+    return inverse_t<Operation, Element>() (op, v);
+}
+
+
+// Short-cut for multiplicative inverse
+template <typename Element>
+inline Element reciprocal(const Element& v)
+{
+    return inverse(math::mult<Element>(), v);
+}
+
+} // namespace math
+
+#endif // MATH_INVERSE_INCLUDE

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

+// Copyright 2006. Peter Gottschling, Matthias Troyer, Rolf Bonderer
+// 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_IS_INVERTIBLE_INCLUDE
+#define MATH_IS_INVERTIBLE_INCLUDE
+
+#include <boost/numeric/linear_algebra/operators.hpp>
+#include <boost/numeric/linear_algebra/identity.hpp>
+
+namespace math {
+
+template <typename Operation, typename Element>
+struct is_invertible_t
+{
+    // bool operator()(const Operation&, const Element&) const;
+};
+
+
+// By default all elements are invertible w.r.t. addition
+// If only part of the elements are invertible it shall be handled by specialization of the type
+// Whether invertibility is relevant at all shall be concrolled by the user with concept maps
+template <typename Element>
+struct is_invertible_t< add<Element>, Element >
+  : public std::binary_function<add<Element>, Element, Element>
+{
+    bool operator() (const add<Element>&, const Element&) const
+    {
+	return true;
+    }
+};
+
+
+// By default all non-zero elements are invertible w.r.t. multiplication
+// If another part of the elements or all elements are invertible it shall be handled by specialization of the type
+// Whether invertibility is relevant at all shall be concrolled by the user with concept maps
+template <typename Element>
+struct is_invertible_t< mult<Element>, Element >
+  : public std::binary_function<mult<Element>, Element, Element>
+{
+    bool operator() (const mult<Element>&, const Element& v) const
+    {
+	return v != zero(v);
+    }
+};
+
+
+// Function is shorter than typetrait-like functor
+template <typename Operation, typename Element>
+inline bool is_invertible(const Operation& op, const Element& v)
+{
+    return is_invertible_t<Operation, Element>() (op, v);
+}
+
+
+namespace detail {
+
+    // Helper type whose operator returns true if v is not 0
+    // 0 must be convertible into Element and Element must be EqualityComparable
+    template <typename Operation, typename Element>
+    struct non_zero_is_invertible_t
+    {
+	bool operator() (const Operation&, const Element& v)
+	{
+	    return v != Element(0);
+	}
+    };
+
+} // namespace detail
+
+} // namespace math
+
+#endif // MATH_IS_INVERTIBLE_INCLUDE

AMDiS/lib/mtl4/boost/numeric/linear_algebra/linear_operator.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_LINEAR_OPERATOR_INCLUDE
+#define MATH_LINEAR_OPERATOR_INCLUDE
+
+#ifdef __GXX_CONCEPTS__
+#  include <concepts>
+#else
+#  include <boost/numeric/linear_algebra/pseudo_concept.hpp>
+#endif
+
+namespace math {
+
+/** @addtogroup Concepts
+ *  @{
+ */
+
+#ifdef __GXX_CONCEPTS__
+    concept LinearOperator<typename Operator, typename VectorDomain, typename VectorImage>
+    {
+	requires VectorSpace<VectorDomain>;
+	requires VectorSpace<VectorImage>;
+
+	typename result_type;
+	result_type operator* (Operator, VectorDomain);
+
+	typename assign_type;
+	assign_type operator= (VectorImage, result_type);
+
+	// The following two requirements are subject to discussion
+	typename plus_assign_type;
+	plus_assign_type operator+= (VectorImage, result_type);
+
+	typename minus_assign_type;
+	minus_assign_type operator-= (VectorImage, result_type);
+
+	axiom Addability(Operator A, VectorDomain x, VectorDomain y)
+	{
+	    A * (x + y) == A*x  + A*y;
+	}
+
+	// The two vector spaces must be scalable with the same scalar types
+	axiom Scalability(Operator A, VectorSpace<VectorDomain>::scalar_type alpha, VectorDomain x)
+	{
+	    A * (alpha * x) == alpha * (A * x);
+	}
+    };
+#else
+    //! Concept LinearOperator
+    /*!
+        Linear operator from one vector space into another one.
+
+        \param Operator The type of the operator, e.g., some matrix type
+	\param VectorDomain The the type of a vector in the domain vector space
+	\param VectorImage The the type of a vector in the image vector space
+
+	\par Associated Types:
+	- result_type
+	- assign_type
+	- plus_assign_type
+	- minus_assign_type
+
+	\par Requires:
+	- VectorSpace < VectorDomain >
+	- VectorSpace < VectorImage >
+
+        \par Notation:
+        <table summary="notation">
+          <tr>
+            <td>A</td>
+	    <td>Object of type Operation</td>
+	  </tr>
+          <tr>
+            <td>x, y</td>
+	    <td>Objects of type VectorDomain</td>
+          </tr>
+          <tr>
+            <td>u</td>
+	    <td>Object of type VectorImage</td>
+          </tr>
+        </table>
+
+        \par Valid Expressions:
+        <table>
+          <tr>
+            <td>Assign product:</td>
+	    <td>u= A * x</td>
+	  </tr>
+          <tr>
+            <td>Add product:</td>
+	    <td>u+= A * x</td>
+	  </tr>
+          <tr>
+            <td>Subtract product:</td>
+	    <td>u-= A * x</td>
+	  </tr>
+        </table>
+
+        \invariant
+        <table summary="invariants">
+          <tr>
+            <td>Addability</td>
+	    <td>A * (x + y) == A*x + A*y</td>
+          </tr>
+          <tr>
+            <td>Scalability</td>
+	    <td>alpha * (A * x) == A * (alpha * x)</td>
+          </tr>
+        </table>
+
+	\note
+	-# Using matrix vector products in arbitrary expressions requires
+	   storing it in temporary objects to avoid redundant computation.
+	   On the other hand, it is not always obvious to choose an appropriate
+	   type for such temporary depending on arbitrary operator and vector types.
+	   Using the products directly in assignments allows implementation without
+	   temporaries, e.g., by calling a function mult(A, x, u) internally.
+     */
+    template <typename Operator, typename VectorDomain, typename VectorImage>
+    struct LinearOperator
+    {
+	/// Associated type: result of multiplication; automatically deducted
+	typedef associated_type result_type;
+	/// Multiplication of linear operator with vector
+	result_type operator* (Operator, VectorDomain);
+
+	/// Associated type: return type of assigning product to vector.
+	/** Automatically deducted. Using expression templates it can be different from VectorImage& */
+	typedef associated_type assign_type;
+	/// Product must be assignable
+	assign_type operator= (VectorImage, result_type);
+
+	// The following two requirements are subject to discussion
+	/// Associated type: return type of incrementally assigning product to vector.
+	/** Automatically deducted. Using expression templates it can be different from VectorImage& */
+	typedef associated_type plus_assign_type;
+	/// Product must be assignable with increment
+	plus_assign_type operator+= (VectorImage, result_type);
+
+	// The following two requirements are subject to discussion
+	/// Associated type: return type of decrementally assigning product to vector.
+	/** Automatically deducted. Using expression templates it can be different from VectorImage& */
+	typedef associated_type minus_assign_type;
+	/// Product must be assignable with decrement
+	minus_assign_type operator+= (VectorImage, result_type);
+
+	/// Invariant: the linear projection of a sum is the sum of the linear projections
+	axiom Addability(Operator A, VectorDomain x, VectorDomain y)
+	{
+	    A * (x + y) == A*x  + A*y;
+	}
+
+	/// Invariant: the linear projection of a scaled vector is the scaling of the vector's linear projections
+	axiom Scalability(Operator A, VectorSpace<VectorDomain>::scalar_type alpha, VectorDomain x)
+	{
+	    A * (alpha * x) == alpha * (A * x);
+	}
+    };
+#endif
+
+
+#ifdef __GXX_CONCEPTS__
+    concept SelfAdjointOperator<typename Operator, typename VectorDomain, typename VectorImage>
+      : LinearOperator<Operator, VectorDomain, VectorImage>
+    {};
+#else
+    //! Concept SelfAdjointOperator
+    /*!
+
+
+        \param Operator The type of the operator, e.g., some matrix type
+	\param VectorDomain The the type of a vector in the domain vector space
+	\param VectorImage The the type of a vector in the image vector space
+
+        \par Refinement of:
+	- LinearOperator <Operator, VectorDomain, VectorImage>
+    */
+    template <typename Operator, typename VectorDomain, typename VectorImage>
+    struct SelfAdjointOperator
+      : LinearOperator<Operator, VectorDomain, VectorImage>
+    {};
+#endif
+
+
+#ifdef __GXX_CONCEPTS__
+    concept RealOperator<typename Operator, typename VectorDomain, typename VectorImage>
+      : LinearOperator<Operator, VectorDomain, VectorImage>
+    {};
+#else
+    //! Concept RealOperator
+    /*!
+
+
+        \param Operator The type of the operator, e.g., some matrix type
+	\param VectorDomain The the type of a vector in the domain vector space
+	\param VectorImage The the type of a vector in the image vector space
+
+        \par Refinement of:
+	- LinearOperator <Operator, VectorDomain, VectorImage>
+    */
+    template <typename Operator, typename VectorDomain, typename VectorImage>
+    struct RealOperator
+      : LinearOperator<Operator, VectorDomain, VectorImage>
+    {};
+#endif
+
+
+#ifdef __GXX_CONCEPTS__
+    concept SymmetricOperator<typename Operator, typename VectorDomain, typename VectorImage>
+      : SelfAdjointOperator<Operator, VectorDomain, VectorImage>,
+        RealOperator<Operator, VectorDomain, VectorImage>
+    {};
+#else
+    //! Concept SymmetricOperator
+    /*!
+
+
+        \param Operator The type of the operator, e.g., some matrix type
+	\param VectorDomain The the type of a vector in the domain vector space
+	\param VectorImage The the type of a vector in the image vector space
+
+        \par Refinement of:
+	- SelfAdjointOperator <Operator, VectorDomain, VectorImage>
+	- RealOperator <Operator, VectorDomain, VectorImage>
+    */
+    template <typename Operator, typename VectorDomain, typename VectorImage>
+    struct SymmetricOperator
+      : SelfAdjointOperator<Operator, VectorDomain, VectorImage>,
+        RealOperator<Operator, VectorDomain, VectorImage>
+    {};
+#endif
+
+/*@}*/ // end of group Concepts
+
+} // namespace math
+
+#endif // MATH_LINEAR_OPERATOR_INCLUDE

AMDiS/lib/mtl4/boost/numeric/linear_algebra/new_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<Operation, Element, Element>
+{
+    axiom Commutativity(Operation op, Element x, Element y)
+    {
+	op(x, y) == op(y, x);
+    }
+};
+
+
+concept SemiGroup<typename Operation, typename Element>
+  : std::Callable2<Operation, 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 >
+{}
+
+
+#ifdef COMPILER_WITHOUT_OVERLOAD_ERROR  // Uncompilable due to error in compiler
+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);
+    }
+};
+
+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 >
+{}
+
+
+#ifdef COMPILER_WITHOUT_OVERLOAD_ERROR  // Uncompilable due to error in compiler
+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);
+    }
+};
+
+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>
+{}
+
+
+
+
+#if 0
+
+// ====================
+// 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/old_concepts.hpp

+// Copyright 2006. Peter Gottschling, Matthias Troyer, Rolf Bonderer
+// 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_CONCEPTS_INCLUDE
+#define LA_CONCEPTS_INCLUDE
+
+#include <boost/config/concept_macros.hpp>
+
+#ifdef __GXX_CONCEPTS__
+#  include <concepts>
+#else
+#  ifdef LA_SHOW_WARNINGS
+#    warning "Concepts are not used"
+#  endif
+#endif
+
+
+#include <boost/numeric/linear_algebra/identity.hpp>
+#include <boost/numeric/linear_algebra/is_invertible.hpp>
+#include <boost/numeric/linear_algebra/inverse.hpp>
+#include <boost/numeric/linear_algebra/operators.hpp>
+#include <boost/numeric/linear_algebra/algebraic_concepts.hpp>
+#include <complex>
+
+// If desired one can disable the default concept maps with LA_NO_CONCEPT_MAPS
+
+// We consider to change the namespace from math to numeric
+// More precisely, the concepts may be moved into namespace numeric and the standard functions stay in math
+
+/// Namespace for mathematical concepts
+/** In contrast to the ones in algebra the concepts can require
+    basic implementation concepts like std::CopyAssignable */
+namespace math {
+
+#ifdef __GXX_CONCEPTS__
+
+// ==================================
+// Classification of Arithmetic Types
+// ==================================
+
+// We omit this now because it became part of the standard concepts
+
+#if 0
+
+// In addtion to std::Integral
+concept Float<typename T>
+  : std::DefaultConstructible<T>, std::CopyConstructible<T>,
+    std::LessThanComparable<T>, std::EqualityComparable<T>
+{
+  T operator+(T);
+  T operator+(T, T);
+  T& operator+=(T&, T);
+  T operator-(T, T);
+  T operator-(T);
+  T& operator-=(T&, T);
+  T operator*(T, T);
+  T& operator*=(T&, T);
+  T operator/(T, T);
+  T& operator/=(T&, T);
+
+  requires std::CopyAssignable<T>
+        && std::SameType<std::CopyAssignable<T>::result_type, T&>;
+}
+
+concept_map Float<float> {}
+concept_map Float<double> {}
+concept_map Float<long double> {}
+
+// The difference to Float is the lack of LessThanComparable
+concept Complex<typename T>
+  : std::DefaultConstructible<T>, std::CopyConstructible<T>,
+    std::EqualityComparable<T>
+{
+  T operator+(T);
+  T operator+(T, T);
+  T& operator+=(T&, T);
+  T operator-(T, T);
+  T operator-(T);
+  T& operator-=(T&, T);
+  T operator*(T, T);
+  T& operator*=(T&, T);
+  T operator/(T, T);
+  T& operator/=(T&, T);
+
+  requires std::CopyAssignable<T>
+        && std::SameType<std::CopyAssignable<T>::result_type, T&>;
+}
+
+template <typename T>
+  requires Float<T>
+concept_map Complex<std::complex<T> > {}
+
+
+// TBD: Concept Arithmetic is useless like this, it should have operations and then be the base for
+// Integral, Float and Complex
+concept Arithmetic<typename T> {}
+
+template <typename T>
+  requires std::Integral<T>
+concept_map Arithmetic<T> {}
+
+template <typename T>
+  requires Float<T>
+concept_map Arithmetic<T> {}
+
+template <typename T>
+  requires Arithmetic<T>
+concept_map Arithmetic< std::complex<T> > {}
+
+#endif
+
+// The following concepts are used to classify intrinsic arithmetic types.
+// The standard concepts already define the syntactic requirements,
+// i.e. the interface.
+// However they sey nothing about the semantics.
+// Therefore, user-defined types can model the syntactic/interface
+// requirements while still having a different mathematical behavior.
+// For that reason, we introduce concepts that are only used for intrinsic types.
+// For them we can define concept_maps regarding semantic behavior as monoids.
+
+concept IntrinsicSignedIntegral<typename T>
+  : std::SignedIntegralLike<T>
+{}
+
+concept IntrinsicUnsignedIntegral<typename T>
+  : std::UnsignedIntegralLike<T>
+{}
+
+concept IntrinsicFloatingPoint<typename T>
+  : std::FloatingPointLike<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> {}
+
+
+
+// ================
+// Utility Concepts
+// ================
+
+// Concepts for functions mapping to same type or convertible
+auto concept UnaryIsoFunction<typename Operation, typename Element>
+{
+    requires std::Callable1<Operation, Element>;
+    requires std::Convertible<std::Callable1<Operation, Element>::result_type, Element>;
+
+    typename result_type = std::Callable1<Operation, Element>::result_type;
+};
+
+
+auto concept BinaryIsoFunction<typename Operation, typename Element>
+{
+    requires std::Callable2<Operation, Element, Element>;
+    requires std::Convertible<std::Callable2<Operation, Element, Element>::result_type, Element>;
+
+    typename result_type = std::Callable2<Operation, Element, Element>::result_type;
+};
+
+#if 0
+auto concept CompatibleBinaryFunction<typename A1, typename A2, typename Result>
+{
+    typename result_type;
+    result_type F(A1, A2);
+    requires std::Convertible<result_type, Result>;
+}
+#endif
+
+// ==================
+// Algebraic Concepts
+// ==================
+
+
+auto concept Magma<typename Operation, typename Element>
+    : BinaryIsoFunction<Operation, Element>
+{
+    requires std::CopyAssignable<Element>
+	  && std::CopyAssignable<Element, BinaryIsoFunction<Operation, Element>::result_type>;
+};
+
+
+// For algebraic structures that are commutative but not associative
+// As an example floating point numbers are commutative but not associative
+//   w.r.t. addition and multiplication
+auto concept CommutativeMagma<typename Operation, typename Element>
+  : Magma<Operation, Element>,
+    algebra::Commutative<Operation, Element>
+{};
+
+
+// SemiGroup is a refinement which must be nominal
+auto concept SemiGroup<typename Operation, typename Element>
+  : Magma<Operation, Element>,
+    algebra::SemiGroup<Operation, Element>
+{};
+
+
+auto concept CommutativeSemiGroup<typename Operation, typename Element>
+  : SemiGroup<Operation, Element>,
+    CommutativeMagma<Operation, Element>
+{};
+
+
+// Adding identity
+// auto
+concept Monoid<typename Operation, typename Element>
+  : SemiGroup<Operation, Element>,
+    algebra::Monoid<Operation, Element>
+{
+    requires std::Convertible<identity_result_type, Element>;
+};
+
+
+auto concept CommutativeMonoid<typename Operation, typename Element>
+  : CommutativeSemiGroup<Operation, Element>,
+    Monoid<Operation, Element>
+{};
+
+
+concept PartiallyInvertibleMonoid<typename Operation, typename Element>
+  : Monoid<Operation, Element>,
+    algebra::Inversion<Operation, Element>
+{
+    typename is_invertible_result_type;
+    is_invertible_result_type is_invertible(Operation, Element);
+    requires std::Convertible<is_invertible_result_type, bool>;
+
+    requires std::Convertible<inverse_result_type, Element>;
+
+    // Does it overwrites the axiom from algebra::Inversion
+    axiom Inversivity(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);
+    }
+};
+
+
+auto concept PartiallyInvertibleCommutativeMonoid<typename Operation, typename Element>
+  : PartiallyInvertibleMonoid<Operation, Element>,
+    CommutativeMonoid<Operation, Element>
+{};
+
+
+concept Group<typename Operation, typename Element>
+  : PartiallyInvertibleMonoid<Operation, Element>,
+    algebra::Group<Operation, Element>
+{
+    axiom AlwaysInvertible(Operation op, Element x)
+    {
+	is_invertible(op, x);
+    }
+
+    axiom GlobalInversivity(Operation op, Element x)
+    {
+	// In fact this is implied by AlwaysInvertible and inherited Inversion axiom
+	// However, we don't rely on the compiler to deduce this
+	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>,
+    PartiallyInvertibleCommutativeMonoid<Operation, Element>,
+    algebra::AbelianGroup<Operation, Element>
+{};
+
+
+// ========================
+// Additive scalar concepts
+// ========================
+
+
+concept AdditiveMagma<typename Element>
+  : Magma< math::add<Element>, Element >
+{
+    typename plus_assign_result_type;
+    plus_assign_result_type operator+=(Element& x, Element y);
+    // requires std::Convertible<plus_assign_result_type, Element>;
+
+    // Operator + is by default defined with +=
+    typename addition_result_type;
+    addition_result_type operator+(Element x, Element y);
+#if 0
+    {
+	Element tmp(x);
+	return tmp += y;                      defaults NYS
+    }
+#endif
+    requires std::Convertible<addition_result_type, Element>;
+
+    // Type consistency with Magma
+    requires std::Convertible< addition_result_type,
+                            Magma< math::add<Element>, Element >::result_type>;
+
+    // SameType requires more rigorous specializations on pure algebraic functors
+    // requires std::SameType< addition_result_type,
+    // 	                    Magma< math::add<Element>, Element >::result_type>;
+
+    axiom Consistency(math::add<Element> op, Element x, Element y)
+    {
+	op(x, y) == x + y;
+
+	// Consistency definition between + and += might change later
+        x + y == x += y;
+	// Element tmp = x; tmp+= y; tmp == x + y; not proposal-compliant
+    }
+}
+
+
+auto concept AdditiveCommutativeMagma<typename Element>
+  : AdditiveMagma<Element>,
+    CommutativeMagma< math::add<Element>, Element >
+{};
+
+
+auto concept AdditiveSemiGroup<typename Element>
+  : AdditiveMagma<Element>,
+    SemiGroup< math::add<Element>, Element >
+{};
+
+
+// We really need only one of the additive concepts for the requirements,
+// the requirements of the other would be implied.
+// Vice versa, to derive concept maps of nested concepts from
+// concept maps of refined concepts, they are needed all.
+auto concept AdditiveCommutativeSemiGroup<typename Element>
+  : AdditiveSemiGroup<Element>,
+    AdditiveCommutativeMagma<Element>,
+    CommutativeSemiGroup< math::add<Element>, Element >
+{};
+
+
+concept AdditiveMonoid<typename Element>
+  : AdditiveSemiGroup<Element>,
+    Monoid< math::add<Element>, Element >
+{
+    Element zero(Element v);
+
+    axiom Consistency (math::add<Element> op, Element x)
+    {
+	zero(x) == identity(op, x);
+    }
+};
+
+
+// We really need only one of the additive concepts for the requirements,
+// the requirements of the other would be implied.
+// Vice versa, to derive concept maps of nested concepts from
+// concept maps of refined concepts, they are needed all.
+auto concept AdditiveCommutativeMonoid<typename Element>
+  : AdditiveMonoid<Element>,
+    AdditiveCommutativeSemiGroup<Element>,
+    CommutativeMonoid< math::add<Element>, Element >
+{};
+
+
+concept AdditivePartiallyInvertibleMonoid<typename Element>
+  : AdditiveMonoid<Element>,
+    PartiallyInvertibleMonoid< math::add<Element>, Element >
+{
+    typename minus_assign_result_type;
+    minus_assign_result_type operator-=(Element& x, Element y);
+    // requires std::Convertible<minus_assign_result_type, Element>;
+
+    // Operator - by default defined with -=
+    typename subtraction_result_type;
+    subtraction_result_type operator-(Element& x, Element y);
+#if 0
+    {
+	Element tmp(x);
+	return tmp -= y;                      defaults NYS
+    }
+#endif
+    requires std::Convertible<subtraction_result_type, Element>;
+
+
+    typename unary_result_type;
+    unary_result_type operator-(Element x);
+#if 0
+    {
+	return zero(x) - x;      defaults NYS
+    }
+#endif
+    requires std::Convertible<unary_result_type, Element>;
+
+    axiom Consistency(math::add<Element> op, Element x, Element y)
+    {
+	// consistency between additive and pure algebraic concept
+	if ( is_invertible(op, y) )
+	    op(x, inverse(op, y)) == x - y;
+	if ( is_invertible(op, y) )
+	    inverse(op, y) == -y;
+
+	// consistency between unary and binary -
+	if ( is_invertible(op, x) )
+	    identity(op, x) - x == -x;
+
+	// Might change later
+	if ( is_invertible(op, y) )
+	    x - y == x -= y;
+	// Element tmp = x; tmp-= y; tmp == x - y; not proposal-compliant
+    }
+
+};
+
+
+auto concept AdditivePartiallyInvertibleCommutativeMonoid<typename Element>
+  : AdditivePartiallyInvertibleMonoid<Element>,
+    AdditiveCommutativeMonoid<Element>,
+    PartiallyInvertibleCommutativeMonoid< math::add<Element>, Element >
+{};
+
+
+
+auto concept AdditiveGroup<typename Element>
+  : AdditivePartiallyInvertibleMonoid<Element>,
+    Group< math::add<Element>, Element >
+{};
+
+
+auto concept AdditiveAbelianGroup<typename Element>
+  : AdditiveGroup<Element>,
+    AdditiveCommutativeMonoid<Element>,
+    AbelianGroup< math::add<Element>, Element >
+{};
+
+
+// ============================
+// Multiplitive scalar concepts
+// ============================
+
+
+concept MultiplicativeMagma<typename Element>
+  : Magma< math::mult<Element>, Element >
+{
+    typename mult_assign_result_type;
+    mult_assign_result_type operator*=(Element& x, Element y);
+    // requires std::Convertible<mult_assign_result_type, Element>;
+
+    // Operator * is by default defined with *=
+    typename mult_result_type;
+    mult_result_type operator*(Element x, Element y);
+#if 0
+    {
+	Element tmp(x);
+	return tmp *= y;                      defaults NYS
+    }
+#endif
+    requires std::Convertible<mult_result_type, Element>;
+
+    // Type consistency with Magma
+    requires std::Convertible< mult_result_type,
+                            Magma< math::mult<Element>, Element >::result_type>;
+
+    // SameType requires more rigorous specializations on pure algebraic functors
+    // requires std::SameType< mult_result_type,
+    // 	                    Magma< math::mult<Element>, Element >::result_type>;
+
+
+    axiom Consistency(math::mult<Element> op, Element x, Element y)
+    {
+	op(x, y) == x * y;
+
+	// Consistency definition between * and *= might change later
+        x * y == x *= y;
+	// Element tmp = x; tmp*= y; tmp == x * y; not proposal-compliant
+    }
+
+}
+
+
+auto concept MultiplicativeSemiGroup<typename Element>
+  : MultiplicativeMagma<Element>,
+    SemiGroup< math::mult<Element>, Element >
+{};
+
+
+auto concept MultiplicativeCommutativeSemiGroup<typename Element>
+  : MultiplicativeSemiGroup<Element>,
+    CommutativeSemiGroup< math::mult<Element>, Element >
+{};
+
+
+concept MultiplicativeMonoid<typename Element>
+  : MultiplicativeSemiGroup<Element>,
+    Monoid< math::mult<Element>, Element >
+{
+    Element one(Element v);
+
+    axiom Consistency (math::mult<Element> op, Element x)
+    {
+	one(x) == identity(op, x);
+    }
+};
+
+
+auto concept MultiplicativeCommutativeMonoid<typename Element>
+  : MultiplicativeMonoid<Element>,
+    MultiplicativeCommutativeSemiGroup<Element>,
+    CommutativeMonoid< math::mult<Element>, Element >
+{};
+
+
+concept MultiplicativePartiallyInvertibleMonoid<typename Element>
+  : MultiplicativeMonoid<Element>,
+    PartiallyInvertibleMonoid< math::mult<Element>, Element >
+{
+    typename divide_assign_result_type;
+    divide_assign_result_type operator/=(Element& x, Element y);
+    // requires std::Convertible<divide_assign_result_type, Element>;
+
+    // Operator / by default defined with /=
+    typename division_result_type = Element;
+    division_result_type operator/(Element x, Element y);
+#if 0
+    {
+	Element tmp(x);
+	return tmp /= y;                      defaults NYS
+    }
+#endif
+    requires std::Convertible<division_result_type, Element>;
+
+    axiom Consistency(math::mult<Element> op, Element x, Element y)
+    {
+	// consistency between multiplicative and pure algebraic concept
+	if ( is_invertible(op, y) )
+	    op(x, inverse(op, y)) == x / y;
+
+	// Consistency between / and /=, might change later
+	if ( is_invertible(op, y) )
+	    x / y == x /= y;
+	// Element tmp = x; tmp/= y; tmp == x / y; not proposal-compliant
+    }
+};
+
+
+auto concept MultiplicativePartiallyInvertibleCommutativeMonoid<typename Element>
+  : MultiplicativePartiallyInvertibleMonoid<Element>,
+    MultiplicativeCommutativeMonoid<Element>,
+    PartiallyInvertibleCommutativeMonoid< math::mult<Element>, Element >
+{};
+
+
+auto concept MultiplicativeGroup<typename Element>
+  : MultiplicativeMonoid<Element>,
+    Group< math::mult<Element>, Element >
+{};
+
+
+auto concept MultiplicativeAbelianGroup<typename Element>
+  : MultiplicativeGroup<Element>,
+    MultiplicativeCommutativeMonoid<Element>,
+    AbelianGroup< math::mult<Element>, Element >
+{};
+
+
+// ======================================
+// Algebraic concepts with two connectors
+// ======================================
+
+// -----------------
+// Based on functors
+// -----------------
+
+// More generic, less handy to use
+
+auto concept GenericRing<typename AddOp, typename MultOp, typename Element>
+  : AbelianGroup<AddOp, Element>,
+    SemiGroup<MultOp, Element>,
+    algebra::Ring<AddOp, MultOp, Element>
+{};
+
+
+auto concept GenericCommutativeRing<typename AddOp, typename MultOp, typename Element>
+  : GenericRing<AddOp, MultOp, Element>,
+    CommutativeSemiGroup<MultOp, Element>
+{};
+
+
+auto concept GenericRingWithIdentity<typename AddOp, typename MultOp, typename Element>
+  : GenericRing<AddOp, MultOp, Element>,
+    Monoid<MultOp, Element>,
+    algebra::RingWithIdentity<AddOp, MultOp, Element>
+{};
+
+
+auto concept GenericCommutativeRingWithIdentity<typename AddOp, typename MultOp, typename Element>
+  : GenericRingWithIdentity<AddOp, MultOp, Element>,
+    GenericCommutativeRing<AddOp, MultOp, Element>,
+    CommutativeMonoid<MultOp, Element>
+{};
+
+
+// auto
+concept GenericDivisionRing<typename AddOp, typename MultOp, typename Element>
+  : GenericRingWithIdentity<AddOp, MultOp, Element>,
+    algebra::DivisionRing<AddOp, MultOp, Element>
+{
+    requires std::Convertible<inverse_result_type, Element>;
+};
+
+
+auto concept GenericField<typename AddOp, typename MultOp, typename Element>
+  : GenericDivisionRing<AddOp, MultOp, Element>,
+    GenericCommutativeRingWithIdentity<AddOp, MultOp, Element>,
+    algebra::Field<AddOp, MultOp, Element>
+{};
+
+
+// ------------------
+// Based on operators
+// ------------------
+
+// Handier, less generic
+
+// Alternative definitions use MultiplicativeMonoid<Element> for Ring
+// and call such concepts Pseudo-Ring
+
+
+auto concept Ring<typename Element>
+  : AdditiveAbelianGroup<Element>,
+    MultiplicativeSemiGroup<Element>,
+    GenericRing<math::add<Element>, math::mult<Element>, Element>
+{};
+
+
+auto concept CommutativeRing<typename Element>
+  : Ring<Element>,
+    MultiplicativeCommutativeSemiGroup<Element>,
+    GenericCommutativeRing<math::add<Element>, math::mult<Element>, Element>
+{};
+
+
+auto concept RingWithIdentity<typename Element>
+  : Ring<Element>,
+    MultiplicativeMonoid<Element>,
+    GenericRingWithIdentity<math::add<Element>, math::mult<Element>, Element>
+{};
+
+
+auto concept CommutativeRingWithIdentity<typename Element>
+  : RingWithIdentity<Element>,
+    CommutativeRing<Element>,
+    MultiplicativeCommutativeMonoid<Element>,
+    GenericCommutativeRingWithIdentity<math::add<Element>, math::mult<Element>, Element>
+{};
+
+
+concept DivisionRing<typename Element>
+  : RingWithIdentity<Element>,
+    MultiplicativePartiallyInvertibleMonoid<Element>,
+    GenericDivisionRing<math::add<Element>, math::mult<Element>, Element>
+{
+    axiom NonZeroDivisibility(Element x)
+    {
+	if (x != zero(x))
+	    x / x == one(x);
+    }
+};
+
+
+auto concept Field<typename Element>
+  : DivisionRing<Element>,
+    CommutativeRingWithIdentity<Element>,
+    GenericField<math::add<Element>, math::mult<Element>, Element>
+{};
+
+
+// ======================
+// Miscellaneous concepts
+// ======================
+
+// that shall find a better place later
+
+
+// EqualityComparable will have the != when defaults are supported
+// At this point the following won't needed anymore
+auto concept FullEqualityComparable<typename T, typename U = T>
+{
+  //requires std::EqualityComparable<T, U>;
+
+    bool operator==(const T&, const U&);
+    bool operator!=(const T&, const U&);
+};
+
+// Closure of EqualityComparable under a binary operation:
+// That is, the result of this binary operation is also EqualityComparable
+// with itself and with the operand type.
+auto concept Closed2EqualityComparable<typename Operation, typename Element>
+  : BinaryIsoFunction<Operation, Element>
+{
+    requires FullEqualityComparable<Element>;
+    requires FullEqualityComparable< BinaryIsoFunction<Operation, Element>::result_type >;
+    requires FullEqualityComparable< Element, BinaryIsoFunction<Operation, Element>::result_type >;
+    requires FullEqualityComparable< BinaryIsoFunction<Operation, Element>::result_type, Element >;
+};
+
+
+// LessThanComparable will have the other operators when defaults are supported
+// At this point the following won't needed anymore
+auto concept FullLessThanComparable<typename T, typename U = T>
+{
+    bool operator<(const T&, const U&);
+    bool operator<=(const T&, const U&);
+    bool operator>(const T&, const U&);
+    bool operator>=(const T&, const U&);
+};
+
+
+// Same for LessThanComparable
+auto concept Closed2LessThanComparable<typename Operation, typename Element>
+  : BinaryIsoFunction<Operation, Element>
+{
+    requires FullLessThanComparable<Element>;
+    requires FullLessThanComparable< BinaryIsoFunction<Operation, Element>::result_type >;
+    requires FullLessThanComparable< Element, BinaryIsoFunction<Operation, Element>::result_type >;
+    requires FullLessThanComparable< BinaryIsoFunction<Operation, Element>::result_type, Element >;
+};
+
+#if 0
+auto concept NumericOperatorResultConvertible<typename T>
+  : AddableWithAssign<T>,
+    SubtractableWithAssign<T>,
+    MultiplicableWithAssign<T>,
+    DivisibleWithAssign<T>
+{
+    requires std::Convertible< AddableWithAssign<T>::result_type, T>;
+    requires std::Convertible< SubtractableWithAssign<T>::result_type, T>;
+    requires std::Convertible< MultiplicableWithAssign<T>::result_type, T>;
+    requires std::Convertible< DivisibleWithAssign<T>::result_type, T>;
+}
+#endif
+
+auto concept AdditionResultConvertible<typename T>
+{
+    typename result_type;
+    result_type operator+(T t, T u);
+    requires std::Convertible<result_type, T>;
+
+    typename result_type;
+    result_type operator+=(T& t, T u);
+    requires std::Convertible<result_type, T>;
+};
+
+
+auto concept SubtractionResultConvertible<typename T>
+{
+    typename result_type;
+    result_type operator-(T t, T u);
+    requires std::Convertible<result_type, T>;
+
+    typename result_type;
+    result_type operator-=(T& t, T u);
+    requires std::Convertible<result_type, T>;
+};
+
+auto concept NumericOperatorResultConvertible<typename T>
+  : AdditionResultConvertible<T>,
+    SubtractionResultConvertible<T>
+{};
+
+// ====================
+// Default Concept Maps
+// ====================
+
+#ifndef LA_NO_CONCEPT_MAPS
+
+// ==============
+// Integral Types
+// ==============
+
+template <typename T>
+  requires IntrinsicSignedIntegral<T>
+concept_map CommutativeRingWithIdentity<T> {}
+
+
+template <typename T>
+  requires IntrinsicUnsignedIntegral<T>
+concept_map AdditiveCommutativeMonoid<T> {}
+
+template <typename T>
+  requires IntrinsicUnsignedIntegral<T>
+concept_map MultiplicativeCommutativeMonoid<T> {}
+
+
+// ====================
+// Floating Point Types
+// ====================
+
+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
+// ===========
+
+// Draft version: defined generously unless there will be problems with some types
+
+template <typename Element>
+concept_map CommutativeMonoid< max<Element>, Element >
+{
+    // Why do we need this?
+    typedef Element identity_result_type;
+}
+
+template <typename Element>
+concept_map CommutativeMonoid< min<Element>, Element >
+{
+    // Why do we need this?
+    typedef Element identity_result_type;
+}
+
+#endif // LA_NO_CONCEPT_MAPS
+
+// Here come some mathematical concepts to be defined later
+
+/// Specify the semantic Behavior of natural numbers (TBD)
+/** Mathematic properties are in most cases only approximated
+    but not held exactly. Rigorous definition would impede
+    usage of real processors. **/
+concept NaturalNumber<typename T> {}
+
+/// Specify the semantic Behavior of integral numbers (TBD)
+/** Mathematic properties are in most cases only approximated
+    but not held exactly. Rigorous definition would impede
+    usage of real processors.
+    It is arguable if this is really a refinement **/
+concept IntegralNumber<typename T> : NaturalNumber<T> {}
+
+/// Specify the semantic Behavior of complex numbers (TBD)
+/** Mathematic properties are in most cases only approximated
+    but not held exactly. Rigorous definition would impede
+    usage of real processors. **/
+concept ComplexNumber<typename T> {}
+
+/// Specify the semantic Behavior of real numbers (TBD)
+/** Mathematic properties are in most cases only approximated
+    but not held exactly. Rigorous definition would impede
+    usage of real processors. **/
+concept RealNumber<typename T> : ComplexNumber<T> {}
+
+#endif // __GXX_CONCEPTS__
+
+
+
+
+// =================================================
+// Concept to specify return type of abs (and norms)
+// =================================================
+
+
+#ifdef __GXX_CONCEPTS__
+
+// Concept to specify to specify projection of scalar value to comparable type
+// For instance as return type of abs
+// Minimalist definition for maximal applicability
+auto concept Magnitude<typename T>
+{
+    typename type = T;
+};
+
+template <typename T>
+concept_map Magnitude<std::complex<T> >
+{
+    typedef T type;
+}
+
+
+// Concept for norms etc., which are real values in mathematical definitions
+auto concept RealMagnitude<typename T>
+  : Magnitude<T>
+{
+    requires FullEqualityComparable<type>;
+    requires FullLessThanComparable<type>;
+
+    requires Field<type>;
+
+    type sqrt(type);
+    // typename sqrt_result;
+    // sqrt_result sqrt(type);
+    // requires std::Convertible<sqrt_result, type>;
+
+    // using std::abs;
+    type abs(T);
+}
+
+#else  // now without concepts
+
+template <typename T>
+struct Magnitude
+{
+    typename type = T;
+};
+
+template <typename T>
+struct Magnitude<std::complex<T> >
+{
+    typedef T type;
+}
+
+template <typename T> struct RealMagnitude
+  : public Magnitude<T>
+{}
+
+#endif  // __GXX_CONCEPTS__
+
+// Type trait version both available with and w/o concepts (TBD: Macro finally :-( )
+// For the moment everything is its own magnitude type, unless stated otherwise
+template <typename T>
+struct magnitude_type_trait
+{
+    typedef T type;
+};
+
+template <typename T>
+struct magnitude_type_trait< std::complex<T> >
+{
+    typedef T type;
+};
+
+
+// =========================================
+// Concepts for convenience (many from Rolf)
+// =========================================
+
+
+#ifdef __GXX_CONCEPTS__
+
+//The following concepts Addable, Subtractable etc. differ from std::Addable, std::Subtractable
+//etc. in so far that no default for result_type is provided, thus allowing automated return type deduction
+
+auto concept Addable<typename T, typename U = T>
+{
+    typename result_type;
+    result_type operator+(const T& t, const U& u);
+};
+
+
+// Usually + and += are both defined
+// + can be efficiently derived from += but not vice versa
+auto concept AddableWithAssign<typename T, typename U = T>
+{
+    typename assign_result_type;
+    assign_result_type operator+=(T& x, U y);
+
+    // Operator + is by default defined with +=
+    typename result_type;
+    result_type operator+(T x, U y);
+#if 0
+    {
+	// Default requires std::CopyConstructible, without default not needed
+	Element tmp(x);
+	return tmp += y;                      defaults NYS
+    }
+#endif
+};
+
+
+auto concept Subtractable<typename T, typename U = T>
+{
+    typename result_type;
+    result_type operator-(const T& t, const U& u);
+};
+
+
+// Usually - and -= are both defined
+// - can be efficiently derived from -= but not vice versa
+auto concept SubtractableWithAssign<typename T, typename U = T>
+{
+    typename assign_result_type;
+    assign_result_type operator-=(T& x, U y);
+
+    // Operator - is by default defined with -=
+    typename result_type;
+    result_type operator-(T x, U y);
+#if 0
+    {
+	// Default requires std::CopyConstructible, without default not needed
+	Element tmp(x);
+	return tmp -= y;                      defaults NYS
+    }
+#endif
+};
+
+
+auto concept Multiplicable<typename T, typename U = T>
+{
+    typename result_type;
+    result_type operator*(const T& t, const U& u);
+};
+
+
+// Usually * and *= are both defined
+// * can be efficiently derived from *= but not vice versa
+auto concept MultiplicableWithAssign<typename T, typename U = T>
+{
+    typename assign_result_type;
+    assign_result_type operator*=(T& x, U y);
+
+    // Operator * is by default defined with *=
+    typename result_type;
+    result_type operator*(T x, U y);
+#if 0
+    {
+	// Default requires std::CopyConstructible, without default not needed
+	Element tmp(x);
+	return tmp *= y;                      defaults NYS
+    }
+#endif
+};
+
+
+auto concept Divisible<typename T, typename U = T>
+{
+    typename result_type;
+    result_type operator / (const T&, const U&);
+};
+
+
+// Usually * and *= are both defined
+// * can be efficiently derived from *= but not vice versa
+auto concept DivisibleWithAssign<typename T, typename U = T>
+{
+    typename assign_result_type;
+    assign_result_type operator*=(T& x, U y);
+
+    // Operator * is by default defined with *=
+    typename result_type;
+    result_type operator*(T x, U y);
+#if 0
+    {
+	// Default requires std::CopyConstructible, without default not needed
+	Element tmp(x);
+	return tmp *= y;                      defaults NYS
+    }
+#endif
+};
+
+
+auto concept Transposable<typename T>
+{
+    typename result_type;
+    result_type trans(T&);
+};
+
+
+// Unary Negation -> Any suggestions for better names?! Is there a word as "negatable"?!
+auto concept Negatable<typename S>
+{
+    typename result_type = S;
+    result_type operator-(const S&);
+};
+
+// Or HasAbs?
+using std::abs;
+auto concept AbsApplicable<typename S>
+{
+    // There are better ways to define abs than the way it is done in std
+    // Likely we replace the using one day
+    typename result_type;
+    result_type abs(const S&);
+};
+
+
+using std::conj;
+auto concept HasConjugate<typename S>
+{
+    typename result_type;
+    result_type conj(const S&);
+};
+
+
+// We need the following; might be placed somewhere else later
+template <Float T>
+concept_map HasConjugate<T>
+{
+    typedef T result_type;
+    result_type conj(const T& s) {return s;}
+}
+
+
+
+// Dot product to be defined:
+auto concept Dottable<typename T, typename U = T>
+{
+    typename result_type = T;
+    result_type dot(const T&t, const U& u);
+};
+
+
+auto concept OneNormApplicable<typename V>
+{
+    typename result_type;
+    result_type one_norm(const V&);
+};
+
+
+auto concept TwoNormApplicable<typename V>
+{
+    typename result_type;
+    result_type two_norm(const V&);
+};
+
+
+auto concept InfinityNormApplicable<typename V>
+{
+    typename result_type;
+    result_type inf_norm(const V&);
+};
+
+
+
+
+#endif  // __GXX_CONCEPTS__
+
+
+} // namespace math
+
+
+
+#endif // LA_CONCEPTS_INCLUDE

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

+// Copyright 2006. Peter Gottschling, Matthias Troyer, Rolf Bonderer
+// 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_OPERATORS_INCLUDE
+#define MATH_OPERATORS_INCLUDE
+
+#include <functional>
+
+#ifndef MATH_DEFAULT_FUNCTORS_WITH_CONCEPTS
+namespace math {
+
+template <typename Element>
+struct add : public std::binary_function<Element, Element, Element>
+{
+    Element operator() (const Element& x, const Element& y)
+    {
+	return x + y;
+    }
+};
+
+
+// Heterogeneous addition, i.e. addends and result type may be different
+template <typename A1, typename A2, typename R>
+struct heterogeneous_add
+  : public std::binary_function<A1, A2, R>
+{
+    R operator() (const A1& x, const A2& y)
+    {
+	return x + y;
+    }
+};
+
+
+// The results of char and short additions are int, dito unsigned
+template <> struct add<char> : heterogeneous_add<char, char, int> {};
+template <> struct add<short> : heterogeneous_add<short, short, int> {};
+template <> struct add<unsigned char> : heterogeneous_add<unsigned char, unsigned char, unsigned int> {};
+template <> struct add<unsigned short> : heterogeneous_add<unsigned short, unsigned short, unsigned int> {};
+
+
+template <typename Element>
+struct mult : public std::binary_function<Element, Element, Element>
+{
+    Element operator() (const Element& x, const Element& y)
+    {
+	return x * y;
+    }
+};
+
+
+template <typename A1, typename A2, typename R>
+struct heterogeneous_mult
+  : public std::binary_function<A1, A2, R>
+{
+    R operator() (const A1& x, const A2& y)
+    {
+	return x * y;
+    }
+};
+
+
+// The results of char and short multiplications are int, dito unsigned
+template <> struct mult<char> : heterogeneous_mult<char, char, int> {};
+template <> struct mult<short> : heterogeneous_mult<short, short, int> {};
+template <> struct mult<unsigned char> : heterogeneous_mult<unsigned char, unsigned char, unsigned int> {};
+template <> struct mult<unsigned short> : heterogeneous_mult<unsigned short, unsigned short, unsigned int> {};
+
+#else
+
+// Now the same with concepts
+
+template <typename Element>
+  requires std::HasPlus<Element>
+struct add : public std::binary_function<Element, Element, result_type>
+{
+    result_type operator() (const Element& x, const Element& y)
+    {
+	return x + y;
+    }
+};
+
+template <typename Element>
+  requires std::HasMultiply<Element>
+struct mult : public std::binary_function<Element, Element, result_type>
+{
+    result_type operator() (const Element& x, const Element& y)
+    {
+	return x * y;
+    }
+};
+
+
+
+
+#endif // MATH_DEFAULT_FUNCTORS_WITH_CONCEPTS
+
+template <typename Element>
+struct min : public std::binary_function<Element, Element, Element>
+{
+    Element operator() (const Element& x, const Element& y)
+    {
+	return x <= y ? x : y;
+    }
+};
+
+
+template <typename Element>
+struct max : public std::binary_function<Element, Element, Element>
+{
+    Element operator() (const Element& x, const Element& y)
+    {
+	return x >= y ? x : y;
+    }
+};
+
+template <typename Element>
+struct bitwise_and : public std::binary_function<Element, Element, Element>
+{
+    Element operator() (const Element& x, const Element& y)
+    {
+	return x & y;
+    }
+};
+
+template <typename Element>
+struct bitwise_or : public std::binary_function<Element, Element, Element>
+{
+    Element operator() (const Element& x, const Element& y)
+    {
+	return x | y;
+    }
+};
+
+template <typename Element>
+struct bitwise_xor : public std::binary_function<Element, Element, Element>
+{
+    Element operator() (const Element& x, const Element& y)
+    {
+	return x ^ y;
+    }
+};
+
+
+} // namespace math
+
+#endif // MATH_OPERATORS_INCLUDE

AMDiS/lib/mtl4/boost/numeric/linear_algebra/power.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_POWER_INCLUDE
+#define MATH_POWER_INCLUDE
+
+#include <concepts>
+#include <boost/numeric/linear_algebra/concepts.hpp>
+#include <boost/numeric/linear_algebra/identity.hpp>
+#include <stdexcept>
+
+
+namespace math {
+
+    template <typename Op, std::Semiregular Element, Integral Exponent>
+        requires std::Callable2<Op, Element, Element>
+              && std::Convertible<std::Callable2<Op, Element, Element>::result_type, Element>
+    inline Element power(const Element& a, Exponent n, Op op)
+    {
+	std::cout << "[Magma] ";
+	if (n < 1) throw std::range_error("power [magma]: n must be > 0");
+
+	Element value= a;
+	for (; n > 1; --n)
+	    value= op(value, a);
+	return value;
+    }
+
+#if 0
+    template <typename Op, std::Semiregular Element, Integral Exponent>
+        requires SemiGroup<Op, Element>
+              && std::Callable2<Op, Element, Element>
+              && std::Convertible<std::Callable2<Op, Element, Element>::result_type, Element>
+    inline Element multiply_and_square_horner(const Element& a, Exponent n, Op op)
+    {
+	if (n < 1) throw std::range_error("mult&square Horner: n must be > 0");
+
+        // Set mask to highest bit
+        Exponent mask= 1 << (8 * sizeof(mask) - 1);
+
+        // If this is a negative number right shift can insert 1s instead of 0s -> infinite loop
+        // Therefore we take the 2nd-highest bit
+        if (mask < 0)
+	    mask= 1 << (8 * sizeof(mask) - 2);
+
+        // Find highest 1 bit
+        while(!bool(mask & n)) mask>>= 1;
+
+        Element value= a;
+        for (mask>>= 1; mask; mask>>= 1) {
+	    value= op(value, value);
+	    if (n & mask)
+		value= op(value, a);
+        }
+        return value;
+    }
+
+    template <typename Op, std::Semiregular Element, Integral Exponent>
+        requires SemiGroup<Op, Element>
+              && std::Callable2<Op, Element, Element>
+              && std::Convertible<std::Callable2<Op, Element, Element>::result_type, Element>
+    inline Element power(const Element& a, Exponent n, Op op)
+    {
+	return multiply_and_square_horner(a, n, op);
+    }
+#endif
+
+
+#if 1
+    // With Horner scheme we can avoid recursion
+    // This one is more intuitive (I believe)
+    template <typename Op, std::Semiregular Element, Integral Exponent>
+        requires SemiGroup<Op, Element>
+              && std::Callable2<Op, Element, Element>
+              && std::Convertible<std::Callable2<Op, Element, Element>::result_type, Element>
+    inline Element power(const Element& a, Exponent n, Op op)
+    {
+	std::cout << "[SemiGroup] ";
+	if (n < 1) throw std::range_error("power [SemiGroup]: n must be > 0");
+
+	Exponent half(n / 2);
+        // If half is 0 then n must be 1 and the result is a
+        if (half == 0)
+	    return a;
+
+        // Compute power of downward rounded exponent and "square" the result
+        Element value= power(a, half, op);
+        value= op(value, value);
+
+        // If n is odd another operation with a is needed
+        if (n & 1)
+	    value= op(value, a);
+        return value;
+    }
+#endif
+
+
+
+    template <typename Op, std::Semiregular Element, Integral Exponent>
+        requires Monoid<Op, Element>
+              && std::Callable2<Op, Element, Element>
+              && std::Convertible<std::Callable2<Op, Element, Element>::result_type, Element>
+    inline Element multiply_and_square(const Element& a, Exponent n, Op op)
+    {
+	// Same as the simpler form except that the first multiplication is made before
+	// the loop and one squaring is saved this way
+	if (n < 0) throw std::range_error("mult&square: n must be >= 0");
+
+	using math::identity;
+	Element value= bool(n & 1) ? Element(a) : Element(identity(op, a)), square= a;
+
+	for (n>>= 1; n > 0; n>>= 1) {
+	    square= op(square, square);
+	    if (n & 1)
+		value= op(value, square);
+	}
+	return value;
+    }
+
+
+    template <typename Op, std::Semiregular Element, Integral Exponent>
+        requires Monoid<Op, Element>
+              && std::Callable2<Op, Element, Element>
+              && std::Convertible<std::Callable2<Op, Element, Element>::result_type, Element>
+    inline Element power(const Element& a, Exponent n, Op op)
+    {
+	std::cout << "[Monoid] ";
+	return multiply_and_square(a, n, op);
+    }
+
+
+
+
+    template <typename Op, std::Semiregular Element, Integral Exponent>
+        requires PIMonoid<Op, Element>
+              && std::Callable2<Op, Element, Element>
+              && std::Convertible<std::Callable2<Op, Element, Element>::result_type, Element>
+    inline Element power(const Element& a, Exponent n, Op op)
+    {
+	std::cout << "[PIMonoid] ";
+	if (n < 0 && !is_invertible(op, a))
+	    throw std::range_error("power [PIMonoid]: a must be invertible with n < 0");
+
+	return n < 0 ? multiply_and_square(Element(inverse(op, a)), Exponent(-n), op)
+	             : multiply_and_square(a, n, op);
+    }
+
+#if 1
+    template <typename Op, std::Semiregular Element, Integral Exponent>
+        requires Group<Op, Element>
+              && std::Callable2<Op, Element, Element>
+              && std::Convertible<std::Callable2<Op, Element, Element>::result_type, Element>
+    inline Element power(const Element& a, Exponent n, Op op)
+    {
+	std::cout << "[Group] ";
+	// For groups we don't need any range test
+
+	return n < 0 ? multiply_and_square(Element(inverse(op, a)), Exponent(-n), op)
+	             : multiply_and_square(a, n, op);
+    }
+#endif
+
+
+#if 0
+    template <typename Op, typename Element, typename Exponent>
+        requires Group<Op, Element>
+              && Integral<Exponent>
+              && std::Semiregular<Element>
+              && std::Callable2<Op, Element, Element>
+              && std::Convertible<std::Callable2<Op, Element, Element>::result_type, Element>
+              && std::Semiregular<math::Inversion<Op, Element>::result_type>
+              && std::HasNegate<Exponent>
+              && math::Monoid<Op, math::Inversion<Op, Element>::result_type>
+              && Integral< std::HasNegate<Exponent>::result_type>
+              && std::Callable2<Op, math::Inversion<Op, Element>::result_type,
+				math::Inversion<Op, Element>::result_type>
+              && std::Convertible<std::Callable2<Op, math::Inversion<Op, Element>::result_type,
+						 math::Inversion<Op, Element>::result_type>::result_type,
+				  math::Inversion<Op, Element>::result_type>
+    inline Element power(const Element& a, Exponent n, Op op)
+    {
+	return n < 0 ? multiply_and_square(inverse(op, a), -n, op)
+	             : multiply_and_square(a, n, op);
+    }
+#endif
+
+} // namespace math
+
+#endif // MATH_POWER_INCLUDE

AMDiS/lib/mtl4/boost/numeric/linear_algebra/pseudo_concept.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 PSEUDO_CONCEPT_INCLUDE
+#define PSEUDO_CONCEPT_INCLUDE
+
+/** @addtogroup Concepts
+ *  @{
+ */
+
+#ifndef __GXX_CONCEPTS__
+
+//! Pseudo type for invariants in concepts
+/// Pseudo type used to document invariants in concepts
+struct axiom {};
+
+//! Pseudo type used to document associated types in concepts
+/// Pseudo type for associated types in concepts
+struct associated_type {};
+
+#endif // __GXX_CONCEPTS__
+
+/*@}*/ // end of group Concepts
+
+#endif // PSEUDO_CONCEPT_INCLUDE

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

+// Copyright 2006. Peter Gottschling, Matthias Troyer, Rolf Bonderer
+// 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_VECTOR_CONCEPTS_INCLUDE
+#define LA_VECTOR_CONCEPTS_INCLUDE
+
+
+#include <boost/numeric/linear_algebra/concepts.hpp>
+#include <boost/numeric/linear_algebra/ets_concepts.hpp>
+
+#ifdef __GXX_CONCEPTS__
+#  include <concepts>
+#else
+#  include <boost/numeric/linear_algebra/pseudo_concept.hpp>
+#endif
+
+
+namespace math {
+
+/** @addtogroup Concepts
+ *  @{
+ */
+
+#ifdef __GXX_CONCEPTS__
+concept VectorSpace<typename Vector, typename Scalar = typename Vector::value_type>
+: AdditiveAbelianGroup<Vector>
+{
+    requires Field<Scalar>;
+    requires Multiplicable<Scalar, Vector>;
+    requires MultiplicableWithAssign<Vector, Scalar>;
+    requires DivisibleWithAssign<Vector, Scalar>;
+
+    requires std::Assignable<Vector, Multiplicable<Scalar, Vector>::result_type>;
+    requires std::Assignable<Vector, Multiplicable<Vector, Scalar>::result_type>;
+    requires std::Assignable<Vector, Divisible<Vector, Scalar>::result_type>;
+
+    // Associated types of Field<Scalar> and AdditiveAbelianGroup<Vector> collide
+    // typename result_type = AdditiveAbelianGroup<Vector>::result_type;
+    // typename assign_result_type = AdditiveAbelianGroup<Vector>::assign_result_type;
+
+    axiom Distributivity(Vector v, Vector w, Scalar a, Scalar b)
+    {
+	a * (v + w) == a * v + a * w;
+	(a + b) * v == a * v + b * v;
+	// The following properties are implied by the above, Field and Abelian group
+	// Can we be sure that compilers can deduce/interfere it?
+	(v + w) * a == v * a + w * a;
+	v * (a + b) == v * a + v * b;
+    }
+}
+#else
+    //! Concept VectorSpace
+    /*!
+	\param Vector   The the type of a vector or a collection
+        \param Scalar   The scalar over which the vector field is defined
+
+	\par Requires:
+	- Field < Scalar >;
+	- Multiplicable <Scalar, Vector>;
+	- MultiplicableWithAssign <Vector, Scalar>;
+	- DivisibleWithAssign <Vector, Scalar>;
+	- std::Assignable <Vector, Multiplicable<Scalar, Vector>::result_type>;
+	- std::Assignable <Vector, Multiplicable<Vector, Scalar>::result_type>;
+	- std::Assignable <Vector, Divisible<Vector, Scalar>::result_type>;
+
+    */
+    template <typename Vector, typename Scalar = typename Vector::value_type>
+    struct VectorSpace
+      : AdditiveAbelianGroup<Vector>
+    {
+	/// Invariant: Distributivity of scalars and vectors from left and from right
+	axiom Distributivity(Vector v, Vector w, Scalar a, Scalar b)
+	{
+	    /// a * (v + w) == a * v + a * w;       // Scalar from left
+
+	    /// Vector from right: (a + b) * v == a * v + b * v;
+
+	    /// Scalar from right: (v + w) * a == v * a + w * a;
+
+	    /// Vector from left:  v * (a + b) == v * a + v * b;
+	}
+    };
+#endif
+
+#ifdef __GXX_CONCEPTS__
+concept Norm<typename N, typename Vector,
+	     typename Scalar = typename Vector::value_type>
+  : std::Callable1<N, Vector>
+{
+    requires VectorSpace<Vector, Scalar>;
+    requires RealMagnitude<Scalar>;
+    typename magnitude_type = MagnitudeType<Scalar>::type;
+    requires std::Convertible<magnitude_type, Scalar>;
+
+    typename result_type_norm = std::Callable1<N, Vector>::result_type;
+    requires std::Convertible<result_type_norm, RealMagnitude<Scalar>::magnitude_type>;
+    requires std::Convertible<result_type_norm, Scalar>;
+
+    // Version with function instead functor, as used by Rolf and Matthias
+    // Axioms there defined without norm functor and concept has only 2 types
+#if 0
+    typename result_type_norm;
+    result_type_norm norm(const Vector&);
+    requires std::Convertible<result_type_norm, magnitude_type>;
+    requires std::Convertible<result_type_norm, Scalar>;
+#endif
+
+    axiom Positivity(N norm, Vector v, magnitude_type ref)
+    {
+	norm(v) >= zero(ref);
+    }
+
+    // The following is covered by RealMagnitude
+    // requires AbsApplicable<Scalar>;
+    // requires std::Convertible<AbsApplicable<Scalar>::result_type, magnitude_type>;
+    // requires Multiplicable<magnitude_type>;
+
+    axiom PositiveHomogeneity(N norm, Vector v, Scalar a)
+    {
+	norm(a * v) == abs(a) * norm(v);
+    }
+
+    axiom TriangleInequality(N norm, Vector u, Vector v)
+    {
+	norm(u + v) <= norm(u) + norm(v);
+    }
+}
+#else
+    //! Concept Norm
+    /*!
+        Semantic requirements of a norm
+
+	\param N        Norm functor
+	\param Vector   The the type of a vector or a collection
+        \param Scalar   The scalar over which the vector field is defined
+
+	\par Refinement of:
+	- std::Callable1 <N, Vector>
+
+	\par Associated types:
+	- magnitude_type
+	- result_type_norm
+
+	\par Requires:
+	- VectorSpace <Vector, Scalar>;
+	- RealMagnitude < Scalar >;
+	- std::Convertible <magnitude_type, Scalar>;
+	- std::Convertible <result_type_norm, RealMagnitude<Scalar>::magnitude_type>;
+	- std::Convertible <result_type_norm, Scalar>;
+
+    */
+template <typename N, typename Vector,
+	  typename Scalar = typename Vector::value_type>
+struct Norm
+  : std::Callable1<N, Vector>
+{
+    /// Associated type to represent real values in teh Field of scalar (with default)
+    /** By default MagnitudeType<Scalar>::type */
+    typedef associated_type magnitude_type;
+
+    /// Associated type for result of norm functor
+    /** Automatically detected */
+    typedef associated_type result_type_norm;
+
+    /// Invariant: norm of vector is larger than zero
+    axiom Positivity(N norm, Vector v, magnitude_type ref)
+    {
+	/// norm(v) >= zero(ref);
+    }
+
+    /// Invariant: positive homogeneity with scalar
+    axiom PositiveHomogeneity(N norm, Vector v, Scalar a)
+    {
+	/// norm(a * v) == abs(a) * norm(v);
+    }
+
+    /// Invariant: triangle inequality
+    axiom TriangleInequality(N norm, Vector u, Vector v)
+    {
+	/// norm(u + v) <= norm(u) + norm(v);
+    }
+};
+#endif
+
+
+#ifdef __GXX_CONCEPTS__
+concept SemiNorm<typename N, typename Vector,
+		 typename Scalar = typename Vector::value_type>
+  : Norm<N, Vector, Scalar>
+{
+    axiom PositiveDefiniteness(N norm, Vector v, magnitude_type ref)
+    {
+	if (norm(v) == zero(ref))
+	    v == zero(v);
+	if (v == zero(v))
+	    norm(v) == zero(ref);
+    }
+}
+#else
+    //! Concept SemiNorm
+    /*!
+        Semantic requirements of a semi-norm
+
+	\param N        Norm functor
+	\param Vector   The the type of a vector or a collection
+        \param Scalar   The scalar over which the vector field is defined
+
+	\par Refinement of:
+	- Norm <N, Vector, Scalar>
+    */
+template <typename N, typename Vector,
+	  typename Scalar = typename Vector::value_type>
+struct SemiNorm
+  : Norm<N, Vector, Scalar>
+{
+    /// The norm of a vector is zero if and only if the vector is the zero vector
+    axiom PositiveDefiniteness(N norm, Vector v, magnitude_type ref)
+    {
+	/// if (norm(v) == zero(ref)) v == zero(v);
+
+	/// if (v == zero(v)) norm(v) == zero(ref);
+    }
+};
+#endif
+
+#ifdef __GXX_CONCEPTS__
+concept BanachSpace<typename N, typename Vector,
+		    typename Scalar = typename Vector::value_type>
+  : Norm<N, Vector, Scalar>,
+    VectorSpace<Vector, Scalar>
+{};
+#else
+    //! Concept BanachSpace
+    /*!
+        A Banach space is a vector space with a norm
+
+	\param N        Norm functor
+	\param Vector   The the type of a vector or a collection
+        \param Scalar   The scalar over which the vector field is defined
+
+	\par Refinement of:
+	- Norm <N, Vector, Scalar>
+	- VectorSpace <Vector, Scalar>
+
+	\note
+	- The (expressible) requirements of Banach Space are already given in Norm.
+	- The difference between the requirements is the completeness of the
+	  Banach space, i.e. that every Cauchy sequence w.r.t. norm(v-w) has a limit
+	  in the space. Unfortunately, completeness is never satisfied for
+	  finite precision arithmetic types.
+	- Another subtle difference is that Norm is not a refinement of Vectorspace
+    */
+template <typename N, typename Vector,
+	  typename Scalar = typename Vector::value_type>
+struct BanachSpace
+  : Norm<N, Vector, Scalar>,
+    VectorSpace<Vector, Scalar>
+{};
+#endif
+
+
+#ifdef __GXX_CONCEPTS__
+concept InnerProduct<typename I, typename Vector,
+		     typename Scalar = typename Vector::value_type>
+  : std::Callable2<I, Vector, Vector>
+{
+    // Result of the inner product must be convertible to Scalar
+    requires std::Convertible<std::Callable2<I, Vector, Vector>::result_type, Scalar>;
+
+    // Let's try without this
+    // requires ets::InnerProduct<I, Vector, Scalar>;
+
+    requires HasConjugate<Scalar>;
+
+    axiom ConjugateSymmetry(I inner, Vector v, Vector w)
+    {
+	inner(v, w) == conj(inner(w, v));
+    }
+
+    axiom SequiLinearity(I inner, Scalar a, Scalar b, Vector u, Vector v, Vector w)
+    {
+	inner(v, b * w) == b * inner(v, w);
+	inner(u, v + w) == inner(u, v) + inner(u, w);
+	// This implies the following (will compilers infere/deduce?)
+	inner(a * v, w) == conj(a) * inner(v, w);
+	inner(u + v, w) == inner(u, w) + inner(v, w);
+    }
+
+    requires RealMagnitude<Scalar>;
+    typename magnitude_type = RealMagnitude<Scalar>::type;
+    // requires FullLessThanComparable<magnitude_type>;
+
+    axiom NonNegativity(I inner, Vector v, MagnitudeType<Scalar>::type magnitude)
+    {
+	// inner(v, v) == conj(inner(v, v)) implies inner(v, v) is real
+	// ergo representable as magnitude type
+	magnitude_type(inner(v, v)) >= zero(magnitude)
+    }
+
+    axiom NonDegeneracy(I inner, Vector v, Vector w, Scalar s)
+    {
+	if (v == zero(v))
+	    inner(v, w) == zero(s);
+	if (inner(v, w) == zero(s))
+	    v == zero(v);
+    }
+};
+#else
+    //! Concept InnerProduct
+    /*!
+        Semantic requirements of a inner product
+
+	\param I        The inner product functor
+	\param Vector   The the type of a vector or a collection
+        \param Scalar   The scalar over which the vector field is defined
+
+	\par Refinement of:
+	- std::Callable2 <I, Vector, Vector>
+
+	\par Associated types:
+	- magnitude_type
+
+	\par Requires:
+	- std::Convertible<std::Callable2 <I, Vector, Vector>::result_type, Scalar> ;
+	  result of inner product convertible to scalar to be used in expressions
+	- HasConjugate < Scalar >
+	- RealMagnitude < Scalar > ; the scalar value needs a real magnitude type
+    */
+template <typename I, typename Vector,
+          typename Scalar = typename Vector::value_type>
+struct InnerProduct
+  : std::Callable2<I, Vector, Vector>
+{
+    /// Associated type: the  real magnitude type of the scalar
+    /** By default RealMagnitude<Scalar>::type */
+    typename associated_type magnitude_type;
+    // requires FullLessThanComparable<magnitude_type>;
+
+    /// The arguments can be changed and the result is then the complex conjugate
+    axiom ConjugateSymmetry(I inner, Vector v, Vector w)
+    {
+	/// inner(v, w) == conj(inner(w, v));
+    }
+
+    /// The inner product is linear in the second argument and conjugate linear in the first one
+    /** The equalities are partly redundant with ConjugateSymmetry */
+    axiom SequiLinearity(I inner, Scalar a, Scalar b, Vector u, Vector v, Vector w)
+    {
+	/// inner(v, b * w) == b * inner(v, w);
+
+	/// inner(u, v + w) == inner(u, v) + inner(u, w);
+
+	/// inner(a * v, w) == conj(a) * inner(v, w);
+
+	/// inner(u + v, w) == inner(u, w) + inner(v, w);
+    }
+
+    /// The inner product of a vector with itself is not negative
+    /** inner(v, v) == conj(inner(v, v)) implies inner(v, v) is representable as real */
+    axiom NonNegativity(I inner, Vector v, MagnitudeType<Scalar>::type magnitude)
+    {
+	/// magnitude_type(inner(v, v)) >= zero(magnitude);
+    }
+
+    /// Non-degeneracy not representable with axiom
+    axiom NonDegeneracy(I inner, Vector v, Vector w, Scalar s)
+    {
+	/// \f$\langle v, w\rangle = 0 \forall w \Leftrightarrow v = \vec{0}\f$
+    }
+};
+#endif
+
+
+
+
+#ifdef __GXX_CONCEPTS_
+// A dot product is only a semantically special case of an inner product
+// Questionable if we want such a concept
+concept DotProduct<typename I, typename Vector,
+		   typename Scalar = typename Vector::value_type>
+  : InnerProduct<I, Vector, Scalar>
+{};
+#else
+    //! Concept DotProduct
+    /*!
+        Semantic requirements of dot product. The dot product is a specific inner product.
+
+	\param I        Norm functor
+	\param Vector   The the type of a vector or a collection
+        \param Scalar   The scalar over which the vector field is defined
+
+	\par Refinement of:
+	- InnerProduct <I, Vector, Scalar>
+    */
+template <typename I, typename Vector,
+	  typename Scalar = typename Vector::value_type>
+struct DotProduct
+  : InnerProduct<I, Vector, Scalar>
+{};
+#endif
+
+
+
+
+// Norm induced by inner product
+// Might be moved to another place later
+// Definition as class and function
+// Conversion from scalar to magnitude_type is covered by norm concept
+template <typename I, typename Vector,
+	  typename Scalar = typename Vector::value_type>
+  _GLIBCXX_WHERE(InnerProduct<I, Vector, Scalar>
+		 && RealMagnitude<Scalar>)
+struct induced_norm_t
+{
+    // Return type evtl. with macro to use concept definition
+    typename magnitude_type_trait<Scalar>::type
+    operator() (const I& inner, const Vector& v)
+    {
+	// Check whether inner product is positive real
+	// assert(Scalar(abs(inner(v, v))) == inner(v, v));
+
+	// Similar check while accepting small imaginary values
+	// assert( (abs(inner(v, v)) - inner(v, v)) / abs(inner(v, v)) < 1e-6; )
+
+	// Could also be defined with abs but that might introduce extra ops
+	// typedef RealMagnitude<Scalar>::type magnitude_type;
+
+	typedef typename magnitude_type_trait<Scalar>::type magnitude_type;
+	return sqrt(static_cast<magnitude_type> (inner(v, v)));
+    }
+};
+
+
+#if 0
+template <typename I, typename Vector,
+	  typename Scalar = typename Vector::value_type>
+  LA_WHERE( InnerProduct<I, Vector, Scalar>
+	    && RealMagnitude<Scalar> )
+magnitude_type_trait<Scalar>::type
+induced_norm(const I& inner, const Vector& v)
+{
+    return induced_norm_t<I, Vector, Scalar>() (inner, v);
+}
+#endif
+
+#ifdef __GXX_CONCEPTS__
+
+
+concept HilbertSpace<typename I, typename Vector,
+		     typename Scalar = typename Vector::value_type,
+		     typename N = induced_norm_t<I, Vector, Scalar> >
+  : InnerProduct<I, Vector, Scalar>,
+    BanachSpace<N, Vector, Scalar>
+{
+    axiom Consistency(Vector v)
+    {
+	math::induced_norm_t<I, Vector, Scalar>()(v) == N()(v);
+    }
+};
+#else
+    //! Concept HilbertSpace
+    /*!
+        A Hilbert space is a vector space with an inner product that induces a norm
+
+	\param I        Inner product functor
+	\param Vector   The the type of a vector or a collection
+        \param Scalar   The scalar over which the vector field is defined
+	\param N        Norm functor
+
+	\par Refinement of:
+	- InnerProduct <I, Vector, Scalar>
+	- BanachSpace <N, Vector, Scalar>
+
+	\note
+	- The (expressible) requirements of Banach Space are already given in InnerProduct
+	  (besides consistency of the functors).
+	- A difference is that InnerProduct is not a refinement of Vectorspace
+    */
+template <typename I, typename Vector,
+          typename Scalar = typename Vector::value_type,
+	  typename N = induced_norm_t<I, Vector, Scalar> >
+struct HilbertSpace
+  : InnerProduct<I, Vector, Scalar>,
+    BanachSpace<N, Vector, Scalar>
+{
+    /// Consistency between norm and induced norm
+    axiom Consistency(Vector v)
+    {
+	/// math::induced_norm_t<I, Vector, Scalar>()(v) == N()(v);
+    }
+};
+#endif // __GXX_CONCEPTS__
+
+/*@}*/ // end of group Concepts
+
+} // namespace math
+
+#endif // LA_VECTOR_CONCEPTS_INCLUDE

AMDiS/lib/mtl4/boost/numeric/meta_math/abs.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 META_MATH_ABS_INCLUDE
+#define META_MATH_ABS_INCLUDE
+
+namespace meta_math {
+
+template <long int x>
+struct abs
+{
+  static long int const value = x < 0 ? -x : x;
+};
+
+
+} // namespace meta_math
+
+#endif // META_MATH_ABS_INCLUDE

AMDiS/lib/mtl4/boost/numeric/meta_math/is_power_of_2.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 META_MATH_IS_POWER_OF_2_INCLUDE
+#define META_MATH_IS_POWER_OF_2_INCLUDE
+
+#include <boost/numeric/meta_math/least_significant_one_bit.hpp>
+
+namespace meta_math {
+
+template <unsigned long X>
+struct is_power_of_2
+{
+    static const bool value= X == least_significant_one_bit<X>::value;
+};
+
+} // namespace meta_math
+
+#endif // META_MATH_IS_POWER_OF_2_INCLUDE

AMDiS/lib/mtl4/boost/numeric/meta_math/is_prime.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 META_MATH_IS_PRIME_INCLUDE
+#define META_MATH_IS_PRIME_INCLUDE
+
+#include <boost/mpl/bool.hpp>
+#include <boost/numeric/meta_math/sqrt.hpp>
+// #include <boost/config/concept_macros.hpp>
+#ifdef __GXX_CONCEPTS__
+#   include <concepts>
+#endif
+
+
+namespace meta_math {
+
+namespace mpl = boost::mpl;
+
+namespace impl {
+
+    // Checks if 'x' is divisible by some odd number >= max_odd, checking decreasingly
+    template <long int x, long int max_odd>
+    struct is_prime_to_max_odd
+    {
+	static bool const value = x % max_odd != 0
+                                  && is_prime_to_max_odd<x, max_odd-2>::value;
+    };
+
+    // Once we reach 1, it's prime
+    template <long int x> struct is_prime_to_max_odd<x, 1> : mpl::true_ {};
+
+
+    // Returns the largest number that x is tried to divided by.
+    // This is a odd number slightly larger than the approximated square root.
+    template <long int x>
+    struct max_prime_compare
+    {
+	static long int const tmp = sqrt<x>::value,
+                              value = tmp % 2 == 0 ? tmp + 1 : tmp + 2;
+    };
+
+
+    // Checks if there is an odd number between 3 and sqrt(x)+1 that divides x
+    // if there is no divisor found then x is prime (otherwise not)
+    // must be disabled when x is even
+    template <long int x, bool Disable>
+    struct check_odd
+    {
+	static bool const value = is_prime_to_max_odd<x, max_prime_compare<x>::value>::value;
+    };
+
+    // Intended for even numbers (> 2) which are always prime
+    template <long int x>
+    struct check_odd<x, true>
+    {
+	static bool const value = false;
+    };
+
+}
+
+template <long int x>
+struct is_prime
+{
+  static bool const value = impl::check_odd<x, x % 2 == 0>::value;
+};
+
+template <> struct is_prime<0> : mpl::false_ {};
+template <> struct is_prime<1> : mpl::false_ {};
+template <> struct is_prime<2> : mpl::true_ {};
+template <> struct is_prime<3> : mpl::true_ {};
+template <> struct is_prime<5> : mpl::true_ {};
+
+
+#ifdef __GXX_CONCEPTS__
+    concept Prime<long int N> {}
+
+    template <long int N>
+        where std::True<is_prime<N>::value>
+    concept_map Prime<N> {}
+#endif
+
+
+
+} // namespace meta_math
+
+#endif // META_MATH_IS_PRIME_INCLUDE

AMDiS/lib/mtl4/boost/numeric/meta_math/least_significant_one_bit.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 META_MATH_LEAST_SIGNIFICANT_ONE_BIT_INCLUDE
+#define META_MATH_LEAST_SIGNIFICANT_ONE_BIT_INCLUDE
+
+namespace meta_math {
+
+template <unsigned long X>
+struct least_significant_one_bit
+{
+    static const unsigned long value= ((X ^ (X-1)) + 1) >> 1;
+};
+
+
+} // namespace meta_math
+
+#endif // META_MATH_LEAST_SIGNIFICANT_ONE_BIT_INCLUDE

AMDiS/lib/mtl4/boost/numeric/meta_math/log_2.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 META_MATH_LOG_2_INCLUDE
+#define META_MATH_LOG_2_INCLUDE
+
+#include <boost/numeric/meta_math/is_power_of_2.hpp>
+
+namespace meta_math {
+
+// Computes the logarithm to the basis 2
+// Without testing if power of 2 it rounds values down to next integer
+template <unsigned long X>
+struct log_2
+{
+    // BOOST_STATIC_ASSERT(is_power_of_2_meta<X>::value);
+    static const unsigned long tmp= X >> 1, value= log_2<tmp>::value + 1;
+};
+
+template <> struct log_2<1>
+{
+    static const unsigned long value= 0;
+};
+
+template <> struct log_2<0>
+{
+  // #error "Logarithm of 0 is undefined"
+  BOOST_STATIC_ASSERT(true); // Logarithm of 0 is undefined
+};
+
+
+} // namespace meta_math
+
+#endif // META_MATH_LOG_2_INCLUDE

AMDiS/lib/mtl4/boost/numeric/meta_math/loop.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 META_MATH_LOOP_INCLUDE
+#define META_MATH_LOOP_INCLUDE
+
+#include <boost/numeric/meta_math/loop1.hpp>
+#include <boost/numeric/meta_math/loop2.hpp>
+#include <boost/numeric/meta_math/loop3.hpp>
+
+#endif // META_MATH_LOOP_INCLUDE

AMDiS/lib/mtl4/boost/numeric/meta_math/loop1.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 META_MATH_LOOP1_INCLUDE
+#define META_MATH_LOOP1_INCLUDE
+
+// See loop3.hpp for example
+
+namespace meta_math {
+
+template <std::size_t Index0, std::size_t Max0>
+struct loop1
+{
+    static std::size_t const index0= Index0 - 1, next_index0= Index0 + 1;
+};
+
+
+template <std::size_t Max0>
+struct loop1<Max0, Max0>
+{
+    static std::size_t const index0= Max0 - 1;
+};
+
+
+} // namespace meta_math
+
+#endif // META_MATH_LOOP1_INCLUDE