AMDiS/lib/mtl4/boost/numeric/mtl/operation/entry_similar.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_ENTRY_SIMILAR_INCLUDE
+#define MTL_ENTRY_SIMILAR_INCLUDE
+
+
+
+#include <boost/numeric/mtl/utility/category.hpp>
+#include <boost/numeric/mtl/concept/collection.hpp>
+
+
+namespace mtl {
+
+    namespace matrix {
+
+	template <typename Matrix, typename Value>
+	bool inline entry_similar(const Matrix& A, typename Collection<Matrix>::size_type i,
+				  typename Collection<Matrix>::size_type j, const Value& v, const Value& eps, tag::universe)
+	{
+	    using std::abs;
+	    return abs(A[i][j] - v) <= eps;
+	}
+
+	/// Compares A[i][j] with \p v returns false if larger than eps
+	/** Function works on distributed matrices where only one process tests and the other return true.
+	    The functions enables writing tests that work in parallel (and all other platforms). **/
+	template <typename Matrix, typename Value>
+        bool inline entry_similar(const Matrix& A, typename Collection<Matrix>::size_type i,
+				  typename Collection<Matrix>::size_type j, const Value& v, const Value& eps)
+	{
+	    return entry_similar(A, i, j, v, eps,  typename mtl::traits::category<Matrix>::type());
+	}
+
+    } // namespacs matrix
+
+    namespace vector {
+
+	template <typename Vector, typename Value>
+	bool inline entry_similar(const Vector& x, typename Collection<Vector>::size_type i,
+				  const Value& v, const Value& eps, tag::universe)
+	{
+	    using std::abs;
+	    return abs(x[i] - v) <= eps;
+	}
+
+	/// Compares x[i] with \p v returns false if larger than eps
+	/** Function works on distributed matrices where only one process tests and the other return true.
+	    The functions enables writing tests that work in parallel (and all other platforms). **/
+	template <typename Vector, typename Value>
+        bool inline entry_similar(const Vector& x, typename Collection<Vector>::size_type i,
+				  const Value& v, const Value& eps)
+	{
+	    return entry_similar(x, i, v, eps,  typename mtl::traits::category<Vector>::type());
+	}
+
+    }
+
+} // namespace mtl
+
+#endif // MTL_ENTRY_SIMILAR_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/evaluate_lazy.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_EVALUATE_LAZY_INCLUDE
+#define MTL_EVALUATE_LAZY_INCLUDE
+
+#include <boost/numeric/mtl/operation/lazy_assign.hpp>
+#include <boost/numeric/mtl/vector/lazy_reduction.hpp>
+
+
+namespace mtl {
+
+/// Overloaded function to evaluate lazy expressions
+template <typename T, typename U, typename Assign>
+void inline evaluate_lazy(lazy_assign<T, U, Assign>& lazy)
+{
+    Assign::first_update(lazy.first, lazy.second);
+}
+
+template <typename T, typename Vector, typename Functor, typename Assign>
+void inline evaluate_lazy(lazy_assign<T, vector::lazy_reduction<Vector, Functor>, Assign>& lazy)
+{
+    lazy.first= Functor::post_reduction(vector::reduction<4, Functor, T>::apply(lazy.second.v));
+}
+
+} // namespace mtl
+
+#endif // MTL_EVALUATE_LAZY_INCLUDE

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

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef MTL_OPERATION_EXTENDED_COMPLEX_INCLUDE
+#define MTL_OPERATION_EXTENDED_COMPLEX_INCLUDE
+
+#include <complex>
+#include <boost/utility/enable_if.hpp>
+#include <boost/type_traits/is_unsigned.hpp>
+#include <boost/numeric/mtl/utility/different_non_complex.hpp>
+#include <boost/numeric/mtl/utility/extended_complex.hpp>
+
+
+// Now we dare writing in the standard namespace
+namespace std {
+
+// ========
+// Addition
+// ========
+
+template <typename T, typename U>
+typename mtl::traits::extended_complex<T, U>::type // implicit enable_if
+inline operator+(const T& x, const complex<U>& y)
+{
+    typedef typename mtl::traits::extended_complex<T, U>::type type;
+    return type(x + real(y), imag(y));
+}
+
+template <typename T, typename U>
+typename mtl::traits::extended_complex<T, U>::type // implicit enable_if
+inline operator+(const complex<T>& x, const U& y)
+{
+    typedef typename mtl::traits::extended_complex<T, U>::type type;
+    return type(real(x) + y, imag(x));
+}
+
+template <typename T, typename U>
+typename mtl::traits::extended_complex<T, U>::type // implicit enable_if
+inline operator+(const complex<T>& x, const complex<U>& y)
+{
+    typedef typename mtl::traits::extended_complex<T, U>::type type;
+    return type(real(x) + real(y), imag(x) + imag(y));
+}
+
+// ===========
+// Subtraction
+// ===========
+
+namespace detail {
+
+    // To avoid stupid warnings on unary - for unsigned int specialize it
+    template <typename T>
+    typename boost::disable_if<boost::is_unsigned<T>, T>::type
+    inline negate_helper(const T& x)
+    { return -x; }
+
+    template <typename T>
+    typename boost::enable_if<boost::is_unsigned<T>, T>::type
+    inline negate_helper(const T& x)
+    { return T(0) - x; }
+}
+
+template <typename T, typename U>
+typename mtl::traits::extended_complex<T, U>::type // implicit enable_if
+inline operator-(const T& x, const complex<U>& y)
+{
+    typedef typename mtl::traits::extended_complex<T, U>::type type;
+    return type(x - real(y), detail::negate_helper(imag(y))); // see above for helper
+}
+
+template <typename T, typename U>
+typename mtl::traits::extended_complex<T, U>::type // implicit enable_if
+inline operator-(const complex<T>& x, const U& y)
+{
+    typedef typename mtl::traits::extended_complex<T, U>::type type;
+    return type(real(x) - y, imag(x));
+}
+
+template <typename T, typename U>
+typename mtl::traits::extended_complex<T, U>::type // implicit enable_if
+inline operator-(const complex<T>& x, const complex<U>& y)
+{
+    typedef typename mtl::traits::extended_complex<T, U>::type type;
+    return type(real(x) - real(y), imag(x) - imag(y));
+}
+
+// ==============
+// Multiplication
+// ==============
+
+template <typename T, typename U>
+typename mtl::traits::extended_complex<T, U>::type // implicit enable_if
+inline operator*(const T& x, const complex<U>& y)
+{
+    typedef typename mtl::traits::extended_complex<T, U>::type type;
+    return type(x * real(y), x * imag(y));
+}
+
+template <typename T, typename U>
+typename mtl::traits::extended_complex<T, U>::type // implicit enable_if
+inline operator*(const complex<T>& x, const U& y)
+{
+    typedef typename mtl::traits::extended_complex<T, U>::type type;
+    return type(real(x) * y, imag(x) * y);
+}
+
+template <typename T, typename U>
+typename mtl::traits::extended_complex<T, U>::type // implicit enable_if
+inline operator*(const complex<T>& x, const complex<U>& y)
+{
+    typedef typename mtl::traits::extended_complex<T, U>::type type;
+    return type(real(x) * real(y) - imag(x) * imag(y),
+		real(x) * imag(y) + imag(x) * real(y));
+}
+
+// ========
+// Division
+// ========
+// Not necessarily the most efficient implementations
+// Dealing primarily with types
+
+template <typename T, typename U>
+typename mtl::traits::extended_complex<T, U>::type // implicit enable_if
+inline operator/(const T& x, const complex<U>& y)
+{
+    typedef typename mtl::traits::extended_complex<T, U>::type type;
+    type r(x);
+    return r/= type(real(y), imag(y));
+}
+
+template <typename T, typename U>
+typename mtl::traits::extended_complex<T, U>::type // implicit enable_if
+inline operator/(const complex<T>& x, const U& y)
+{
+    typedef typename mtl::traits::extended_complex<T, U>::type type;
+    return type(real(x) / y, imag(x) / y);
+}
+
+template <typename T, typename U>
+typename mtl::traits::extended_complex<T, U>::type // implicit enable_if
+inline operator/(const complex<T>& x, const complex<U>& y)
+{
+    typedef typename mtl::traits::extended_complex<T, U>::type type;
+    type r(real(x), imag(x));
+    return r/= type(real(y), imag(y));
+}
+
+
+// ==========
+// Comparison
+// ==========
+
+template <typename T, typename U>
+typename boost::enable_if<mtl::traits::different_non_complex<T, U>, bool>::type
+inline operator==(const complex<T>& x, const U& y)
+{
+    return real(x) == y && imag(x) == T(0);
+}
+
+template <typename T, typename U>
+typename boost::enable_if<mtl::traits::different_non_complex<T, U>, bool>::type
+inline operator==(const T& x, const complex<U>& y)
+{
+    return x == real(y) && imag(y) == T(0);
+}
+
+template <typename T, typename U>
+typename boost::enable_if<mtl::traits::different_non_complex<T, U>, bool>::type
+inline operator==(const complex<T>& x, const complex<U>& y)
+{
+    return real(x) == real(y) && imag(x) == imag(y);
+}
+
+template <typename T, typename U>
+typename boost::enable_if<mtl::traits::different_non_complex<T, U>, bool>::type
+inline operator!=(const complex<T>& x, const U& y)
+{
+    return real(x) != y || imag(x) != T(0);
+}
+
+template <typename T, typename U>
+typename boost::enable_if<mtl::traits::different_non_complex<T, U>, bool>::type
+inline operator!=(const T& x, const complex<U>& y)
+{
+    return x != real(y) || imag(y) != T(0);
+}
+
+template <typename T, typename U>
+typename boost::enable_if<mtl::traits::different_non_complex<T, U>, bool>::type
+inline operator!=(const complex<T>& x, const complex<U>& y)
+{
+    return real(x) != real(y) || imag(x) != imag(y);
+}
+
+} // namespace std
+
+#endif // MTL_OPERATION_EXTENDED_COMPLEX_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/fill.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_FILL_INCLUDE
+#define MTL_FILL_INCLUDE
+
+#include <algorithm>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+namespace mtl {
+
+// Sets all values of a matrix
+template <class Matrix>
+void fill(Matrix& ma, typename Matrix::value_type value)
+{
+    vampir_trace<3009> tracer;
+    std::fill(ma.elements(), ma.elements()+ma.num_elements(), value);
+}
+
+// Temporary solution
+// will be replaced by sequences and cursors generated by begin<all>(ma) and end<all>(ma)
+// Using segmented cursors, matrices with non-contigous element storing can be handled
+
+} // namespace mtl
+
+#endif // MTL_FILL_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/frobenius_norm.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_FROBENIUS_NORM_INCLUDE
+#define MTL_FROBENIUS_NORM_INCLUDE
+
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/utility/is_row_major.hpp>
+#include <boost/numeric/mtl/utility/tag.hpp>
+#include <boost/numeric/mtl/utility/range_generator.hpp>
+#include <boost/numeric/mtl/utility/property_map.hpp>
+#include <boost/numeric/mtl/operation/max_of_sums.hpp>
+#include <boost/numeric/mtl/operation/squared_abs.hpp>
+#include <boost/numeric/linear_algebra/identity.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+
+namespace mtl { namespace matrix {
+
+/// Frobenius norm, i.e. square root of sum of squares of all entries sqrt(sum_i sum_j(|a[i][j]|^2))
+template <typename Matrix>
+typename RealMagnitude<typename Collection<Matrix>::value_type>::type
+inline frobenius_norm(const Matrix& matrix)
+{
+    vampir_trace<3010> tracer;
+    using std::sqrt; using std::abs; using math::zero;
+    namespace traits = mtl::traits;
+    typename traits::const_value<Matrix>::type     value(matrix);
+
+    typedef typename Collection<Matrix>::value_type   value_type;
+    typedef typename RealMagnitude<value_type>::type  real_type;
+    real_type ref, sum= zero(ref);
+
+    typedef typename traits::range_generator<tag::major, Matrix>::type     cursor_type;
+    typedef typename traits::range_generator<tag::nz, cursor_type>::type   icursor_type;
+
+    for (cursor_type cursor = begin<tag::major>(matrix), cend = end<tag::major>(matrix); cursor != cend; ++cursor)
+	for (icursor_type icursor = begin<tag::nz>(cursor), icend = end<tag::nz>(cursor); icursor != icend; ++icursor)
+	    sum+= squared_abs(value(*icursor));
+    return sqrt(sum);
+}
+
+}} // namespace mtl::matrix
+
+#endif // MTL_FROBENIUS_NORM_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/fuse.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_FUSE_INCLUDE
+#define MTL_FUSE_INCLUDE
+
+#include <boost/utility/enable_if.hpp>
+#include <boost/mpl/and.hpp>
+
+#include <boost/numeric/mtl/utility/is_lazy.hpp>
+#include <boost/numeric/mtl/operation/fused_expr.hpp>
+
+namespace mtl {
+
+/// Fuse two lazy expressions (operator notation of fuse)
+template <typename T, typename U>
+typename boost::enable_if<boost::mpl::and_<traits::is_lazy<T>, traits::is_lazy<U> >, fused_expr<T, U> >::type
+operator||(const T& x, const U& y)
+{
+    return fused_expr<T, U>(const_cast<T&>(x), const_cast<U&>(y));
+}
+
+/// Fuse two lazy expressions
+template <typename T, typename U>
+typename boost::enable_if<boost::mpl::and_<traits::is_lazy<T>, traits::is_lazy<U> >, fused_expr<T, U> >::type
+fuse(const T& x, const U& y)
+{
+    return fused_expr<T, U>(const_cast<T&>(x), const_cast<U&>(y));
+}
+
+
+
+} // namespace mtl
+
+#endif // MTL_FUSE_INCLUDE

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

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef MTL_FUSED_EXPR_INCLUDE
+#define MTL_FUSED_EXPR_INCLUDE
+
+#include <boost/mpl/bool.hpp>
+#include <boost/mpl/and.hpp>
+#include <boost/numeric/mtl/operation/index_evaluator.hpp>
+#include <boost/numeric/mtl/utility/index_evaluatable.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+namespace mtl {
+
+/// Expression template for fusing other expression
+template <typename T, typename U>
+struct fused_expr
+{
+    typedef boost::mpl::and_<traits::forward_index_evaluatable<T>, traits::forward_index_evaluatable<U> >   forward;
+    typedef boost::mpl::and_<traits::backward_index_evaluatable<T>, traits::backward_index_evaluatable<U> > backward;
+
+    fused_expr(T& first, U& second)
+      : first(first), second(second), delayed_assign(false)
+    {
+	first.delay_assign(); second.delay_assign();
+    }
+
+    ~fused_expr() { if (!delayed_assign) eval(forward(), backward()); }
+
+    void delay_assign() const { delayed_assign= true; }
+
+    template <typename TT, typename UU>
+    void forward_eval_loop(const TT& const_first_eval, const UU& const_second_eval, boost::mpl::false_)
+    {
+	vampir_trace<6001> tracer;
+	// hope there is a more elegant way; copying the arguments causes errors due to double destructor evaluation
+	TT& first_eval= const_cast<TT&>(const_first_eval);
+	UU& second_eval= const_cast<UU&>(const_second_eval);
+	MTL_DEBUG_THROW_IF(mtl::vector::size(first_eval) != mtl::vector::size(second_eval), incompatible_size());
+
+	for (std::size_t i= 0, s= size(first_eval); i < s; i++) {
+	    first_eval(i); second_eval(i);
+	}
+    }
+
+    template <typename TT, typename UU>
+    void forward_eval_loop(const TT& const_first_eval, const UU& const_second_eval, boost::mpl::true_)
+    {
+	vampir_trace<6002> tracer;
+	// hope there is a more elegant way; copying the arguments causes errors due to double destructor evaluation
+	TT& first_eval= const_cast<TT&>(const_first_eval);
+	UU& second_eval= const_cast<UU&>(const_second_eval);
+	MTL_DEBUG_THROW_IF(mtl::vector::size(first_eval) != mtl::vector::size(second_eval), incompatible_size());
+
+	const std::size_t s= size(first_eval), sb= s >> 2 << 2;
+
+	for (std::size_t i= 0; i < sb; i+= 4) {
+	    first_eval.template at<0>(i); second_eval.template at<0>(i);
+	    first_eval.template at<1>(i); second_eval.template at<1>(i);
+	    first_eval.template at<2>(i); second_eval.template at<2>(i);
+	    first_eval.template at<3>(i); second_eval.template at<3>(i);
+	}
+
+	for (std::size_t i= sb; i < s; i++) {
+	    first_eval(i); second_eval(i);
+	}
+    }
+
+    // Forward evaluation dominates backward
+    template <bool B2>
+    void eval(boost::mpl::true_, boost::mpl::bool_<B2>)
+    {
+#ifdef MTL_LAZY_LOOP_WO_UNROLL
+	typedef boost::mpl::false_                                                                              to_unroll;
+#else
+	typedef boost::mpl::and_<traits::unrolled_index_evaluatable<T>, traits::unrolled_index_evaluatable<U> > to_unroll;
+#endif
+	// Currently lazy evaluation is only available on vector expressions, might change in the future
+	// std::cout << "Forward evaluation\n";
+	forward_eval_loop(index_evaluator(first), index_evaluator(second), to_unroll());
+    }
+
+    template <typename TT, typename UU>
+    void backward_eval_loop(const TT& const_first_eval, const UU& const_second_eval, boost::mpl::false_)
+    {
+	vampir_trace<6003> tracer;
+	// hope there is a more elegant way; copying the arguments causes errors due to double destructor evaluation
+	TT& first_eval= const_cast<TT&>(const_first_eval);
+	UU& second_eval= const_cast<UU&>(const_second_eval);
+	MTL_DEBUG_THROW_IF(mtl::vector::size(first_eval) != mtl::vector::size(second_eval), incompatible_size());
+
+	for (std::size_t i= size(first_eval); i-- > 0; ) {
+	    // std::cout << "i is " << i << "\n";
+	    first_eval(i); second_eval(i);
+	}
+    }
+
+    template <typename TT, typename UU>
+    void backward_eval_loop(const TT& const_first_eval, const UU& const_second_eval, boost::mpl::true_)
+    {
+	vampir_trace<6004> tracer;
+	// hope there is a more elegant way; copying the arguments causes errors due to double destructor evaluation
+	TT& first_eval= const_cast<TT&>(const_first_eval);
+	UU& second_eval= const_cast<UU&>(const_second_eval);
+	MTL_DEBUG_THROW_IF(mtl::vector::size(first_eval) != mtl::vector::size(second_eval), incompatible_size());
+
+	std::size_t s= size(first_eval), i= s-1, m= s % 4;
+	for (; m; m--) {
+	    // std::cout << "i is " << i << "\n";
+	    first_eval(i); second_eval(i--);
+	}
+	for(long j= i - 3; j >= 0; j-= 4) {
+	    // std::cout << "i is " << j+3 << ".." << j << "\n";
+	    first_eval.template at<3>(j); second_eval.template at<3>(j);
+	    first_eval.template at<2>(j); second_eval.template at<2>(j);
+	    first_eval.template at<1>(j); second_eval.template at<1>(j);
+	    first_eval.template at<0>(j); second_eval.template at<0>(j);
+	}
+    }
+
+    // Backward evaluation if forward isn't possible
+    void eval(boost::mpl::false_, boost::mpl::true_)
+    {
+	// std::cout << "Backward evaluation\n";
+	backward_eval_loop(index_evaluator(first), index_evaluator(second), boost::mpl::true_());
+    }
+
+    // Sequential evaluation
+    void eval(boost::mpl::false_, boost::mpl::false_)
+    {
+	// std::cout << "Non-fused evaluation\n";
+	evaluate_lazy(first); evaluate_lazy(second);
+    }
+
+    T& first;
+    U& second;
+    mutable bool                        delayed_assign;
+};
+
+
+} // namespace mtl
+
+#endif // MTL_FUSED_EXPR_INCLUDE

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

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG, www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also tools/license/license.mtl.txt in the distribution.
+
+#ifndef MTL_FUSED_INDEX_EVALUATOR_INCLUDE
+#define MTL_FUSED_INDEX_EVALUATOR_INCLUDE
+
+#include <boost/numeric/mtl/utility/index_evaluator.hpp>
+
+namespace mtl { namespace vector {
+
+template <typename T, typename U>
+struct fused_index_evaluator
+{
+    fused_index_evaluator(T& first, U& second)
+      : first(index_evaluator(first)), second(index_evaluator(second)) {}
+
+    template <unsigned Offset>
+    void at(std::size_t i)
+    {
+	first.template at<Offset>(i);
+	second.template at<Offset>(i);
+    }
+
+    void operator() (std::size_t i) { at<0>(i); }
+    void operator[] (std::size_t i) { at<0>(i); }
+
+    typename traits::index_evaluator<T>::type first;
+    typename traits::index_evaluator<U>::type second;
+};
+
+template <typename T, typename U>
+inline size_t size(const fused_index_evaluator<T, U>& expr) { return size(expr.first); }
+
+}} // namespace mtl::vector
+
+#endif // MTL_FUSED_INDEX_EVALUATOR_INCLUDE

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

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+// With contributions from Cornelius Steinhardt
+
+#ifndef MTL_MATRIX_GIVENS_INCLUDE
+#define MTL_MATRIX_GIVENS_INCLUDE
+
+#include <cmath>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/mtl/utility/irange.hpp>
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/operation/householder.hpp>
+#include <boost/numeric/mtl/operation/rank_one_update.hpp>
+#include <boost/numeric/mtl/operation/trans.hpp>
+
+namespace mtl { namespace matrix {
+
+/// Given's transformator
+/** Requires Hessenberg form, i.e. for transformations near the diagonal.
+    For general form use qr_givens.
+    \sa qr_givens. **/
+template <typename Matrix>
+class givens
+{
+    typedef typename Collection<Matrix>::value_type   value_type;
+    typedef typename Collection<Matrix>::size_type    size_type;
+
+  public:
+    /// Re-set the rotation parameters \p a and \p b
+    void set_rotation(value_type a, value_type b)
+    {
+	using std::abs;
+	value_type zero= math::zero(a), one= math::one(b), t;
+
+	if ( b == zero ) {
+	    c= one; s= zero;
+	} else if ( abs(b) > abs(a) ) {
+	    t= -a / b;
+	    s= one / sqrt(one + t*t);
+	    c= s * t;
+	} else {
+	    t= -b / a;
+	    c= one / sqrt(one + t*t);
+	    s= c * t;
+	}
+	G= c, s,
+	  -s, c;
+    }
+
+
+    /// Constructor takes %matrix \p H to be transformed and the rotation parameters \p a and \p b
+    givens(Matrix& H, value_type a, value_type b) : H(H), G(2, 2)
+    {	set_rotation(a, b);    }
+
+    /// Given's transformation of \p H with \p G regarding column \p k
+    Matrix& trafo(const Matrix& G, size_type k)
+    {
+	    irange r(k,k+2);
+	    // trans(H[r][ind])*= G; H[ind][r]*= G; // most compact form but does not work yet
+
+	    Matrix col_block(H[r][iall]), col_perm(trans(G) * col_block);
+	    H[r][iall]= col_perm;
+	    Matrix row_perm(H[iall][r] * G);
+	    H[iall][r]= row_perm;
+
+	    return H;
+    }
+
+    /// Given's transformation of \p H regarding column \p k
+    Matrix& trafo(size_type k)
+    {
+	return trafo(G, k);
+    }
+
+  private:
+    Matrix&    H, G;
+    value_type c, s;
+};
+
+}// namespace matrix
+
+
+namespace vector {
+
+/// Given's transformator on %vector (swap a*line(k) with b*line(k+1) )
+template <typename Vector>
+class givens
+{
+    typedef typename Collection<Vector>::value_type   value_type;
+    typedef typename Collection<Vector>::size_type    size_type;
+
+  public:
+    /// Constructor takes %vector \p v to be transformed and the rotation parameters \p a and \p b
+    givens(Vector& v, value_type a, value_type b) : v(v), a(a), b(b)
+    {  }
+
+    /// Given's transformation of \p v with \p a and \p b regarding column \p k
+    Vector& trafo(size_type k)
+    {
+	    value_type w1(0), w2(0);
+	    w1= a*v[k] - b*v[k+1]; //given's rotation on solution
+            w2= b*v[k] + a*v[k+1]; //rotation on vector
+            v[k]= w1;
+            v[k+1]= w2;
+
+	    return v;
+    }
+
+  private:
+    Vector&    v;
+    value_type a, b;
+};
+
+}// namespace vector
+
+
+} // namespace mtl
+
+#endif // MTL_MATRIX_GIVENS_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/hermitian.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_HERMITIAN_INCLUDE
+#define MTL_HERMITIAN_INCLUDE
+
+#include <boost/numeric/mtl/matrix/hermitian_view.hpp>
+#include <boost/numeric/mtl/matrix/view_ref.hpp>
+#include <boost/numeric/mtl/utility/enable_if.hpp>
+#include <boost/numeric/mtl/utility/view_code.hpp>
+#include <boost/numeric/mtl/utility/viewed_collection.hpp>
+#include <boost/numeric/mtl/utility/compose_view.hpp>
+
+namespace mtl {
+
+    // vector version to be done
+
+    namespace matrix {
+
+	namespace detail {
+
+	    template <typename Matrix>
+	    struct hermitian
+	    {
+		static const unsigned code= mtl::traits::view_toggle_hermitian<mtl::traits::view_code<Matrix> >::value;
+		typedef typename mtl::traits::compose_view<code, typename mtl::traits::viewed_collection<Matrix>::type>::type result_type;
+
+		static inline result_type apply(const Matrix& A)
+		{
+		    return result_type(view_ref(A));
+		}
+	    };
+
+	} // namespace detail
+
+	/// Return hermitian of matrix A
+	template <typename Matrix>
+	typename mtl::traits::enable_if_matrix<Matrix, typename detail::hermitian<Matrix>::result_type >::type
+	inline hermitian(const Matrix& A)
+	{
+	    return detail::hermitian<Matrix>::apply(A);
+	}
+    }
+
+    using matrix::hermitian;
+
+} // namespace mtl
+
+#endif // MTL_HERMITIAN_INCLUDE

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

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+// With contributions from Cornelius Steinhardt
+
+#ifndef MTL_MATRIX_HESSENBERG_INCLUDE
+#define MTL_MATRIX_HESSENBERG_INCLUDE
+
+#include <cmath>
+#include <boost/numeric/linear_algebra/identity.hpp>
+#include <boost/numeric/mtl/vector/parameter.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/mtl/utility/irange.hpp>
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/concept/magnitude.hpp>
+#include <boost/numeric/mtl/operation/householder.hpp>
+#include <boost/numeric/mtl/operation/rank_one_update.hpp>
+#include <boost/numeric/mtl/operation/trans.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+
+namespace mtl { namespace matrix {
+
+/// Hessenberg-Factorization of matrix A with householder-vectors
+/** Return Hessenberg matrix and tril(B,-2) are Householder-vectors **/
+template <typename Matrix>
+Matrix inline hessenberg_factors(const Matrix& A)
+{
+    vampir_trace<5014> tracer;
+    if (num_rows(A) < 3)
+	return A;
+
+    using mtl::imax;
+    typedef typename Collection<Matrix>::value_type   value_type;
+    typedef typename Collection<Matrix>::size_type    size_type;
+    size_type        ncols = num_cols(A), nrows = num_rows(A);
+    value_type       zero= math::zero(A[0][0]), beta;
+    Matrix           B(clone(A));
+
+    for(size_type i= 0; i < ncols-2; i++){
+	// mtl::vector::dense_vector<value_type>  v(B[irange(i+1, imax)][i]);
+		mtl::vector::dense_vector<value_type, vector::parameters<> >  v(nrows-i-1), w(nrows);
+        for (size_type j = 0; j < size(v); j++)
+            v[j]= B[j+i+1][i];
+        beta= householder(v).second;
+        v= householder(v).first;
+;
+	if( beta != zero){
+            w= beta * B[irange(0,imax)][irange(i+1,imax)] * v;
+	    //rank_one_update(B[irange(0,imax)][irange(i+1,imax)],-w,v);
+            for(size_type row = 0; row < nrows; row++){
+                for(size_type col = i+1; col < ncols; col++){
+                    B[row][col] -= w[row] * v[col-i-1];
+                }
+            }
+            //vector*Matrix
+            for(size_type k=0; k < size(w); k++){
+                w[k]= zero;
+                for(size_type j = 0; j < size(v); j++){
+                    w[k] += beta * v[j] * B[j+i+1][k];
+                }
+            }
+            //rank_one_update(A[irange(i+1,imax)][irange(0,imax)],-v,w);
+            for(size_type row = i+1; row < nrows; row++){
+                for(size_type col = 0; col < ncols; col++){
+                    B[row][col] -= v[row-i-1] * w[col];
+                }
+            }
+	    // B[irange(i+2, imax)][i]= v[irange(1, nrows-i-1)];
+            for(size_type row = i+2; row < nrows; row++){
+                B[row][i] = v[row-i-1];
+            }
+        }
+    }
+    return B;
+}
+
+
+/// Extract Householder vectors from Hessenberg factorization H of some A
+template <typename Matrix>
+Matrix inline extract_householder_hessenberg(const Matrix& H)
+{
+    vampir_trace<5015> tracer;
+    return Matrix(tril(H, -2));
+}
+
+/// Compute Householder vectors from Hessenberg factorization of A
+template <typename Matrix>
+Matrix inline householder_hessenberg(const Matrix& A)
+{
+    vampir_trace<5016> tracer;
+    return Matrix(tril(hessenberg_factors(A), -2));
+}
+
+
+/// Extract Hessenberg form from factorization H of some A
+template <typename Matrix>
+Matrix inline extract_hessenberg(const Matrix& H)
+{
+    vampir_trace<5017> tracer;
+    return Matrix(triu(H, -1));
+}
+
+/// Hessenberg form of A
+template <typename Matrix>
+Matrix inline hessenberg(const Matrix& A)
+{
+    vampir_trace<5018> tracer;
+    // return triu(hessenberg_factors(A), -1);
+    MTL_THROW_IF(num_rows(A) < 3, matrix_too_small());
+
+    typedef typename Collection<Matrix>::value_type   value_type;
+    // typedef typename Magnitude<value_type>::type      magnitude_type; // to multiply with 2 not 2+0i
+    typedef typename Collection<Matrix>::size_type    size_type;
+    size_type        ncols = num_cols(A), nrows = num_rows(A);
+    value_type       zero= math::zero(A[0][0]);
+    Matrix           H(nrows,ncols);
+
+    H= hessenberg_factors(A);
+
+    // H= bands(hessenberg_factors(A), -nrows, -1);
+    // set (doubly) strict lower triangle to zero
+    for(size_type row = 2; row < nrows; row++){
+        for(size_type col = 0; col < row-1; col++){
+            H[row][col]= zero;
+        }
+    }
+
+    return H;
+}
+
+
+/// Return Q where Q'*A*Q == hessenberg(A)
+template <typename Matrix>
+Matrix inline hessenberg_q(const Matrix& A)
+{
+    vampir_trace<5013> tracer;
+    using std::abs;
+    typedef typename Collection<Matrix>::value_type   value_type;
+    typedef typename Magnitude<value_type>::type      magnitude_type; // to multiply with 2 not 2+0i
+    typedef typename Collection<Matrix>::size_type    size_type;
+    size_type        ncols = num_cols(A), nrows = num_rows(A), mini;
+    value_type       zero= math::zero(A[0][0]), one= math::one(A[0][0]);
+    const magnitude_type two(2);
+    Matrix           Q(nrows,ncols);
+
+    MTL_THROW_IF(num_rows(A) < 3, matrix_too_small());
+
+    Q= one;
+
+    //Extract Q
+    for(size_type i = 0; i < nrows-2; i++){
+		mtl::vector::dense_vector<value_type, vector::parameters<> >   v(nrows-1), w(nrows);
+        v[0]= one;
+// 	std::cout<< "v=" << v << "\n";
+// 	std::cout<< "w=" << w << "\n";
+        for(size_type k = 1; k < size(v); k++)
+            v[k]= A[nrows-k][i];
+
+	magnitude_type beta= two / abs(dot(v, v)); // abs: x+0i -> x
+        if (beta != two) {
+            //trans(Vector)*Matrix
+            for(size_type k = 0; k < size(w); k++) {
+                w[k]= zero;
+                for(size_type j = 0; j < size(v); j++){
+// 		     std::cout<< "k=" << k << "  j=" << j << "\n";
+                    w[k]+= beta * v[j] * Q[j+i+1][k];
+		}
+            }
+            //rank_one_update(Q[irange(i+1,imax)][irange(0,imax)],-v,w);
+            for(size_type row = i+1; row < nrows; row++)
+                for(size_type col = 0; col < ncols; col++){
+// 		    std::cout<< "row=" << row << "  col=" << col << "\n";
+                    Q[row][col] -= v[row-1] * w[col];
+		}
+        }
+    }
+    return Q;
+}
+
+
+}} // namespace mtl::matrix
+
+#endif // MTL_MATRIX_HESSENBERG_INCLUDE
+

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

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+// With contributions from Cornelius Steinhardt
+
+#ifndef MTL_MATRIX_HOUSEHOLDER_INCLUDE
+#define MTL_MATRIX_HOUSEHOLDER_INCLUDE
+
+#include <cmath>
+#include <cassert>
+#include <boost/numeric/mtl/mtl_fwd.hpp>
+#include <boost/numeric/mtl/vector/parameter.hpp>
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+
+namespace mtl { namespace vector {
+
+
+/// Computes Householder vector v and scalar b for vector \p y
+/** such that identity_matrix(size(y))-b*v*v' projects the vector y
+    to a positive multiple of the first unit vector. **/
+template <typename Vector>
+std::pair<typename mtl::vector::dense_vector<typename Collection<Vector>::value_type, parameters<> >, typename Collection<Vector>::value_type>
+inline householder(Vector& y)
+{
+    vampir_trace<2004> tracer;
+    assert(size(y) > 0);
+    typedef typename  Collection<Vector>::value_type   value_type;
+    // typedef typename  Collection<Vector>::size_type    size_type;
+    const value_type  zero= math::zero(y[0]), one= math::one(y[0]);
+
+    Vector            v(y);
+    v[0]= one;
+    irange            tail(1, imax);
+    value_type        s( dot(v[tail], v[tail]) ), b, v0;
+
+    //evaluation of v and b
+    if (s == zero)
+        b= zero;
+    else {
+	value_type mu= sqrt(y[0] * y[0] + s);
+	v0= v[0]= y[0] <= zero ? y[0] - mu : -s / (y[0] + mu); // complex < zero????
+	b= 2 * v0 * v0 / (s + v0 * v0);                        // 2* complex???????
+	v/= v0;                                                // normalization of the first entry
+    }
+
+    return std::make_pair(v,b);
+}
+
+/// Computes Householder vector for vector \p y
+/** More stable Householder transformation, also for non-square matrices. **/
+template <typename Vector>
+typename mtl::vector::dense_vector<typename Collection<Vector>::value_type, parameters<> >
+inline householder_s(Vector& y)
+{
+    vampir_trace<2005> tracer;
+    typedef typename  Collection<Vector>::value_type   value_type;
+    const value_type  zero= math::zero(y[0]);
+
+    Vector            u(y);
+    value_type        nu(sqrt( dot(u, u) )), s;
+
+    if (nu != zero){
+	if(u[0] < 0){
+		s= -1;
+	} else {
+		s= 1;
+	}
+	u[0]= u[0] + s * nu;
+	u/= sqrt( dot(u, u) );
+    }
+
+    return u;
+}
+
+
+}} // namespace mtl::matrix
+
+#endif // MTL_MATRIX_HOUSEHOLDER_INCLUDE
+

AMDiS/lib/mtl4/boost/numeric/mtl/operation/imag.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_IMAG_INCLUDE
+#define MTL_IMAG_INCLUDE
+
+#include <complex>
+#include <boost/numeric/linear_algebra/identity.hpp>
+
+namespace mtl {
+
+namespace sfunctor {
+
+    template <typename Value>
+    struct imag
+    {
+	typedef Value result_type;
+
+	static inline result_type apply(const Value& v)
+	{
+	    using math::zero;
+	    return zero(v);
+	}
+	result_type operator()(const Value& v) const
+	{
+	    using math::zero;
+	    return zero(v);
+	}
+    };
+
+    template <typename Value>
+    struct imag<std::complex<Value> >
+    {
+	typedef Value result_type;
+
+	static inline result_type apply(const std::complex<Value>& v)
+	{
+	    return std::imag(v);
+	}
+	result_type operator()(const std::complex<Value>& v) const
+	{
+	    return std::imag(v);
+	}
+    };
+
+}
+
+/// Imaginary part of scalars (including non-complex)
+template <typename Value>
+typename mtl::traits::enable_if_scalar<Value, typename sfunctor::imag<Value>::result_type>::type
+inline imag(const Value& v)
+{
+    return sfunctor::imag<Value>::apply(v);
+}
+
+namespace vector {
+
+    /// Imaginary part of an vector
+    template <typename Vector>
+    typename mtl::traits::enable_if_vector<Vector, imag_view<Vector> >::type
+    inline imag(const Vector& v)
+    {
+	return imag_view<Vector>(v);
+    }
+}
+
+namespace matrix {
+
+    /// Imaginary part of an vector
+    template <typename Matrix>
+    typename mtl::traits::enable_if_matrix<Matrix, imag_view<Matrix> >::type
+    inline imag(const Matrix& v)
+    {
+	return imag_view<Matrix>(v);
+    }
+}
+
+
+} // namespace mtl
+
+#endif // MTL_IMAG_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/index_evaluator.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_INDEX_EVALUATOR_INCLUDE
+#define MTL_INDEX_EVALUATOR_INCLUDE
+
+#include <boost/utility/enable_if.hpp>
+#include <boost/mpl/and.hpp>
+
+#include <boost/numeric/mtl/mtl_fwd.hpp>
+#include <boost/numeric/mtl/operation/lazy_assign.hpp>
+#include <boost/numeric/mtl/utility/category.hpp>
+#include <boost/numeric/mtl/vector/vec_vec_aop_expr.hpp>
+#include <boost/numeric/mtl/vector/vec_scal_aop_expr.hpp>
+#include <boost/numeric/mtl/vector/reduction_index_evaluator.hpp>
+#include <boost/numeric/mtl/vector/dot_index_evaluator.hpp>
+#include <boost/numeric/mtl/vector/row_mat_cvec_index_evaluator.hpp>
+#include <boost/numeric/mtl/operation/fused_index_evaluator.hpp>
+#include <boost/numeric/mtl/operation/fused_expr.hpp>
+#include <boost/numeric/itl/pc/ic_0.hpp>
+#include <boost/numeric/itl/pc/ilu_0.hpp>
+
+namespace mtl {
+
+/// Overloaded function that maps from lazy expressions to the according index-wise evaluation classes
+template <typename T, typename U, typename Assign>
+typename boost::enable_if<boost::mpl::and_<mtl::traits::is_vector<T>, mtl::traits::is_vector<U> >,
+			  mtl::vector::vec_vec_aop_expr<T, U, Assign> >::type
+inline index_evaluator(lazy_assign<T, U, Assign>& lazy)
+{
+    return mtl::vector::vec_vec_aop_expr<T, U, Assign>(lazy.first, lazy.second, true);
+}
+
+template <typename T, typename U, typename Assign>
+typename boost::enable_if<boost::mpl::and_<mtl::traits::is_vector<T>, mtl::traits::is_scalar<U> >,
+			  mtl::vector::vec_scal_aop_expr<T, U, Assign> >::type
+inline index_evaluator(lazy_assign<T, U, Assign>& lazy)
+{
+    return mtl::vector::vec_scal_aop_expr<T, U, Assign>(lazy.first, lazy.second, true);
+}
+
+template <typename Scalar, typename Vector, typename Functor, typename Assign>
+mtl::vector::reduction_index_evaluator<Scalar, Vector, Functor, Assign>
+inline index_evaluator(lazy_assign<Scalar, mtl::vector::lazy_reduction<Vector, Functor>, Assign>& lazy)
+{
+    return mtl::vector::reduction_index_evaluator<Scalar, Vector, Functor, Assign>(lazy.first, lazy.second.v);
+}
+
+template <typename Scalar, unsigned long Unroll, typename Vector1,
+	  typename Vector2, typename ConjOpt, typename Assign>
+mtl::vector::dot_index_evaluator<Scalar, Vector1, Vector2, ConjOpt, Assign>
+inline index_evaluator(lazy_assign<Scalar, mtl::vector::dot_class<Unroll, Vector1, Vector2, ConjOpt>, Assign>& lazy)
+{
+    return mtl::vector::dot_index_evaluator<Scalar, Vector1, Vector2, ConjOpt, Assign>(lazy.first, lazy.second.v1, lazy.second.v2);
+}
+
+template <typename VectorOut, typename Matrix, typename VectorIn, typename Assign>
+mtl::vector::row_mat_cvec_index_evaluator<VectorOut, Matrix, VectorIn, Assign>
+inline index_evaluator(lazy_assign<VectorOut, mtl::mat_cvec_times_expr<Matrix, VectorIn>, Assign>& lazy)
+{
+    return mtl::vector::row_mat_cvec_index_evaluator<VectorOut, Matrix, VectorIn, Assign>(lazy.first, lazy.second.first, lazy.second.second);
+}
+
+template <typename V1, typename Matrix, typename Value, typename V2>
+itl::pc::ic_0_evaluator<V1, itl::pc::solver<itl::pc::ic_0<Matrix, Value>, V2, true> >
+inline index_evaluator(lazy_assign<V1, itl::pc::solver<itl::pc::ic_0<Matrix, Value>, V2, true>, assign::assign_sum>& lazy)
+{
+    return itl::pc::ic_0_evaluator<V1, itl::pc::solver<itl::pc::ic_0<Matrix, Value>, V2, true> >(lazy.first, lazy.second);
+}
+
+// reuse ic_0_evaluator with ilu_0 preconditioner since IC(0) and ILU(0) do the same at the upper triangle ;-)
+template <typename V1, typename Factorizer, typename Matrix, typename Value, typename V2>
+itl::pc::ic_0_evaluator<V1, itl::pc::solver<itl::pc::ilu<Matrix, Factorizer, Value>, V2, true> >
+inline index_evaluator(lazy_assign<V1, itl::pc::solver<itl::pc::ilu<Matrix, Factorizer, Value>, V2, true>, assign::assign_sum>& lazy)
+{
+    return itl::pc::ic_0_evaluator<V1, itl::pc::solver<itl::pc::ilu<Matrix, Factorizer, Value>, V2, true> >(lazy.first, lazy.second);
+}
+
+template <typename T, typename U>
+inline mtl::vector::fused_index_evaluator<T, U> index_evaluator(fused_expr<T, U>& expr)
+{  return mtl::vector::fused_index_evaluator<T, U>(expr.first, expr.second); }
+
+} // namespace mtl
+
+#endif // MTL_INDEX_EVALUATOR_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/infinity_norm.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_INFINITY_NORM_INCLUDE
+#define MTL_INFINITY_NORM_INCLUDE
+
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/concept/magnitude.hpp>
+#include <boost/numeric/mtl/utility/enable_if.hpp>
+#include <boost/numeric/mtl/utility/is_row_major.hpp>
+#include <boost/numeric/mtl/utility/category.hpp>
+#include <boost/numeric/mtl/utility/tag.hpp>
+#include <boost/numeric/mtl/utility/property_map.hpp>
+#include <boost/numeric/mtl/operation/max_of_sums.hpp>
+#include <boost/numeric/mtl/vector/lazy_reduction.hpp>
+#include <boost/numeric/mtl/vector/reduction.hpp>
+#include <boost/numeric/mtl/vector/reduction_functors.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+
+namespace mtl {
+
+    namespace vector {
+
+	template <unsigned long Unroll, typename Vector>
+	typename traits::enable_if_vector<Vector, typename RealMagnitude<typename Collection<Vector>::value_type>::type>::type
+	inline infinity_norm(const Vector& vector)
+	{
+	    vampir_trace<2006> tracer;
+	    typedef typename RealMagnitude<typename Collection<Vector>::value_type>::type result_type;
+	    return vector::reduction<Unroll, vector::infinity_norm_functor, result_type>::apply(vector);
+	}
+
+	/*! Infinity-norm for vectors: infinity_norm(x) \f$\rightarrow |x|_\infty\f$.
+	    \retval The magnitude type of the respective value type, see Magnitude.
+
+	    The norms are defined as \f$|v|_\infty=\max_i |v_i|\f$.
+
+	    Vector norms are unrolled 8-fold by default.
+	    An n-fold unrolling can be generated with infinity_norm<n>(x).
+	    The maximum for n is 8 (it might be increased later).
+	**/
+	template <typename Vector>
+	typename mtl::traits::enable_if_vector<Vector, typename RealMagnitude<typename Collection<Vector>::value_type>::type>::type
+	inline infinity_norm(const Vector& vector)
+	{
+	    return infinity_norm<8>(vector);
+	}
+
+	template <typename Vector>
+	lazy_reduction<Vector, infinity_norm_functor> inline lazy_infinity_norm(const Vector& v)
+	{  return lazy_reduction<Vector, infinity_norm_functor>(v); 	}
+    }
+
+    namespace matrix {
+
+	// Ignore unrolling for matrices
+	template <unsigned long Unroll, typename Matrix>
+	typename mtl::traits::enable_if_matrix<Matrix, typename RealMagnitude<typename Collection<Matrix>::value_type>::type>::type
+	inline infinity_norm(const Matrix& matrix)
+	{
+	    vampir_trace<3011> tracer;
+	    using mtl::impl::max_of_sums;
+	    typename mtl::traits::row<Matrix>::type                             row(matrix);
+	    return max_of_sums(matrix, mtl::traits::is_row_major<typename OrientedCollection<Matrix>::orientation>(),
+			       row, num_rows(matrix));
+	}
+
+	/*! Infinity-norm for matrices: infinity_norm(x) \f$\rightarrow |x|_\infty\f$.
+	    \retval The magnitude type of the respective value type, see Magnitude.
+
+	    The norms are defined as \f$|A|_\infty=\max_i\{\sum_j(|A_{ij}|)\}\f$.
+	    Matrix norms are not (yet) optimized by unrolling.
+	**/
+	template <typename Matrix>
+	typename mtl::traits::enable_if_matrix<Matrix, typename RealMagnitude<typename Collection<Matrix>::value_type>::type>::type
+	inline infinity_norm(const Matrix& matrix)
+	{
+	    return infinity_norm<8>(matrix);
+	}
+    }
+
+    using vector::infinity_norm;
+    using vector::lazy_infinity_norm;
+    using matrix::infinity_norm;
+
+} // namespace mtl
+
+#endif // MTL_INFINITY_NORM_INCLUDE

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

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef MTL_MATRIX_INV_INCLUDE
+#define MTL_MATRIX_INV_INCLUDE
+
+#include <boost/numeric/mtl/mtl_fwd.hpp>
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/matrix/identity.hpp>
+#include <boost/numeric/mtl/operation/upper_trisolve.hpp>
+#include <boost/numeric/mtl/operation/lu.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/mtl/utility/irange.hpp>
+#include <boost/numeric/mtl/utility/tag.hpp>
+#include <boost/numeric/mtl/vector/parameter.hpp>
+#include <boost/numeric/mtl/vector/dense_vector.hpp>
+#include <boost/numeric/mtl/vector/unit_vector.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+namespace mtl { namespace matrix {
+
+namespace traits {
+
+    /// Return type of inv(Matrix)
+    /** Might be specialized later for the sake of efficiency **/
+    template <typename Matrix>
+    struct inv
+    {
+	typedef typename Collection<Matrix>::value_type    value_type;
+	typedef ::mtl::matrix::dense2D<value_type>         type;
+    };
+
+} // traits
+
+/// Invert upper triangular matrix
+template <typename Matrix, typename MatrixOut>
+void inv_upper(Matrix const& A, MatrixOut& Inv)
+{
+    vampir_trace<5019> tracer;
+    typedef typename Collection<Matrix>::value_type    value_type;
+    typedef typename Collection<Matrix>::size_type     size_type;
+
+    const size_type N= num_rows(A);
+    MTL_DEBUG_THROW_IF(num_cols(A) != N, matrix_not_square());
+    MTL_DEBUG_THROW_IF(N != num_rows(Inv) || num_cols(A) != num_cols(Inv), incompatible_size());
+
+    Inv= math::zero(value_type());
+
+    for (size_type k= 0; k < N; ++k) {
+	irange r(k+1);
+	typename mtl::ColumnInMatrix<MatrixOut>::type col_k(Inv[r][k]);
+	upper_trisolve(A[r][r], vector::unit_vector<value_type>(k, k+1), col_k, mtl::tag::regular_diagonal());
+    }
+}
+
+
+/// Invert upper triangular matrix
+template <typename Matrix>
+inline typename traits::inv<Matrix>::type
+inv_upper(Matrix const& A)
+{
+    typedef typename Collection<Matrix>::size_type     size_type;
+    const size_type N= num_rows(A);
+    typename traits::inv<Matrix>::type Inv(N, N);
+    inv_upper(A, Inv);
+    return Inv;
+}
+
+#if 0
+/// Invert lower triangular matrix
+template <typename Matrix, typename MatrixOut>
+inline void inv_lower(Matrix const& A, MatrixOut& Inv)
+{
+    vampir_trace<5020> tracer;
+    typedef typename Collection<Matrix>::value_type    value_type;
+    typedef typename Collection<Matrix>::size_type     size_type;
+
+    const size_type N= num_rows(A);
+    MTL_DEBUG_THROW_IF(num_cols(A) != N, matrix_not_square());
+    MTL_DEBUG_THROW_IF(N != num_rows(Inv) || num_cols(A) != num_cols(Inv), incompatible_size());
+
+    Inv= math::zero(value_type());
+
+    for (size_type k= 0; k < N; ++k) {
+	irange r(k, N);
+	typename mtl::ColumnInMatrix<MatrixOut>::type col_k(Inv[r][k]);
+	lower_trisolve(A[r][r], vector::unit_vector<value_type>(0, N-k), col_k, mtl::tag::regular_diagonal());
+    }
+}
+
+template <typename Matrix>
+typename traits::inv<Matrix>::type
+inline inv_lower(Matrix const& A)
+{
+    typedef typename Collection<Matrix>::size_type     size_type;
+    const size_type N= num_rows(A);
+    typename traits::inv<Matrix>::type Inv(N, N);
+    inv_lower(A, Inv);
+    return Inv;
+}
+#endif
+
+#if 1
+/// Invert lower triangular matrix
+template <typename Matrix>
+typename traits::inv<Matrix>::type
+inline inv_lower(Matrix const& A)
+{
+    vampir_trace<5020> tracer;
+    Matrix T(trans(A)); // Shouldn't be needed
+    return typename traits::inv<Matrix>::type(trans(inv_upper(T)));
+}
+#endif
+
+
+
+
+/// Invert matrix
+/** Uses pivoting LU factorization and triangular inversion
+    \sa \ref lu, \ref inv_upper, \ref inv_lower **/
+template <typename Matrix, typename MatrixOut>
+inline void inv(Matrix const& A, MatrixOut& Inv)
+{
+    vampir_trace<5021> tracer;
+    typedef typename Collection<Matrix>::size_type     size_type;
+    typedef typename Collection<Matrix>::value_type    value_type;
+    typedef typename traits::inv<Matrix>::type         result_type;
+
+    const size_type N= num_rows(A);
+    MTL_THROW_IF(num_cols(A) != num_cols(A), matrix_not_square());
+    MTL_DEBUG_THROW_IF(N != num_rows(Inv) || num_cols(A) != num_cols(Inv), incompatible_size());
+
+    if (N == 1) {
+	Inv[0][0]= value_type(1) / A[0][0];
+	return;
+    }
+
+    result_type                    PLU(A);
+    mtl::vector::dense_vector<size_type, vector::parameters<> >   Pv(num_rows(A));
+
+    lu(PLU, Pv);
+    result_type  PU(upper(PLU)), PL(strict_lower(PLU));
+    for (size_type i= 0; i < num_rows(A); i++)
+	PL[i][i]= value_type(1);
+
+    Inv= inv_upper(PU) * inv_lower(PL) * permutation(Pv);
+}
+
+template <typename Matrix>
+typename traits::inv<Matrix>::type
+inline inv(Matrix const& A)
+{
+    typedef typename Collection<Matrix>::size_type     size_type;
+    const size_type N= num_rows(A);
+    typename traits::inv<Matrix>::type Inv(N, N);
+    inv(A, Inv);
+    return Inv;
+}
+
+}} // namespace mtl::matrix
+
+#endif // MTL_MATRIX_INV_INCLUDE

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

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef MTL_INVERT_DIAGONAL_INCLUDE
+#define MTL_INVERT_DIAGONAL_INCLUDE
+
+#include <boost/numeric/mtl/utility/tag.hpp>
+#include <boost/numeric/mtl/utility/range_generator.hpp>
+#include <boost/numeric/mtl/utility/property_map.hpp>
+#include <boost/numeric/mtl/matrix/compressed2D.hpp>
+#include <boost/numeric/linear_algebra/inverse.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+namespace mtl { namespace matrix {
+
+///Returns \p A with invert diagonal
+template <typename Matrix>
+inline void invert_diagonal(Matrix& A)
+{
+    vampir_trace<3012> tracer;
+    using math::reciprocal;
+    namespace traits = mtl::traits;
+
+    typename traits::row<Matrix>::type       row(A);
+    typename traits::col<Matrix>::type       col(A);
+    typename traits::value<Matrix>::type     value(A);
+
+    typedef typename traits::range_generator<tag::major, Matrix>::type     cursor_type;
+    typedef typename traits::range_generator<tag::nz, cursor_type>::type   icursor_type;
+
+    for (cursor_type cursor = begin<tag::major>(A), cend = end<tag::major>(A); cursor != cend; ++cursor)
+	for (icursor_type icursor = begin<tag::nz>(cursor), icend = end<tag::nz>(cursor); icursor != icend; ++icursor)
+	    if (row(*icursor) == col(*icursor))
+		value(*icursor, reciprocal(value(*icursor)));
+}
+
+// Specialization for compressed2D
+template <typename Value, typename Parameters>
+inline void invert_diagonal(compressed2D<Value, Parameters>& A)
+{
+    vampir_trace<3026> tracer;
+    using math::reciprocal;
+    for (typename Parameters::size_type i= 0, end= std::min(num_rows(A), num_cols(A)); i < end; i++)
+	A.lvalue(i, i)= reciprocal(A.lvalue(i, i));
+}
+
+
+}} // namespace mtl::matrix
+
+#endif // MTL_INVERT_DIAGONAL_INCLUDE

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

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef MTL_VECTOR_IOTA_INCLUDE
+#define MTL_VECTOR_IOTA_INCLUDE
+
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+namespace mtl { namespace vector {
+
+/// Assigns sequentially increasing values to %vector v
+template <typename Vector>
+void iota(Vector& v, const typename Collection<Vector>::value_type offset= 0)
+{
+    vampir_trace<3013> tracer;
+    for (typename Collection<Vector>::size_type i= 0; i < size(v); i++)
+	v[i]= i + offset;
+}
+
+}} // namespace mtl::vector
+
+#endif // MTL_VECTOR_IOTA_INCLUDE

AMDiS/lib/mtl4/boost/numeric/mtl/operation/is_negative.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_IS_NEGATIVE_INCLUDE
+#define MTL_IS_NEGATIVE_INCLUDE
+
+#include <boost/type_traits/is_unsigned.hpp>
+#include <boost/utility/enable_if.hpp>
+
+namespace mtl {
+
+template <typename T>
+typename boost::enable_if<boost::is_unsigned<T>, bool>::type
+inline is_negative(T) { return false; }
+
+template <typename T>
+typename boost::disable_if<boost::is_unsigned<T>, bool>::type
+inline is_negative(T x) { return x < 0; }
+
+} // namespace mtl
+
+#endif // MTL_IS_NEGATIVE_INCLUDE

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

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef MTL_LAZY_INCLUDE
+#define MTL_LAZY_INCLUDE
+
+#include <boost/numeric/mtl/operation/assign_mode.hpp>
+
+namespace mtl {
+
+/// Helper class for lazy evaluation
+template <typename T>
+struct lazy_t
+{
+    lazy_t(T& data) : data(data) {}
+
+    template <typename U>
+    lazy_assign<T, U, assign::assign_sum> operator=(const U& other)
+    { return lazy_assign<T, U, assign::assign_sum>(data, other); }
+
+    template <typename U>
+    lazy_assign<T, U, assign::plus_sum> operator+=(const U& other)
+    { return lazy_assign<T, U, assign::plus_sum>(data, other); }
+
+    template <typename U>
+    lazy_assign<T, U, assign::minus_sum> operator-=(const U& other)
+    { return lazy_assign<T, U, assign::minus_sum>(data, other); }
+
+    T& data;
+};
+
+/// Do not compute or assign x immediately but delay it until explicitly performed
+template <typename T>
+inline lazy_t<T> lazy(T& x)
+{ return lazy_t<T>(x); }
+
+/// Do not compute or assign x immediately but delay it until explicitly performed
+template <typename T>
+inline lazy_t<const T> lazy(const T& x)
+{ return lazy_t<const T>(x); }
+
+
+
+} // namespace mtl
+
+#endif // MTL_LAZY_INCLUDE