AMDiS/lib/mtl4/boost/numeric/itl/krylov/cg.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef ITL_CG_INCLUDE
+#define ITL_CG_INCLUDE
+
+#include <cmath>
+#include <cassert>
+#include <iostream>
+#include <boost/mpl/bool.hpp>
+
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/itl/itl_fwd.hpp>
+#include <boost/numeric/itl/iteration/basic_iteration.hpp>
+#include <boost/numeric/itl/pc/identity.hpp>
+#include <boost/numeric/itl/pc/is_identity.hpp>
+#include <boost/numeric/mtl/operation/dot.hpp>
+#include <boost/numeric/mtl/operation/unary_dot.hpp>
+#include <boost/numeric/mtl/operation/conj.hpp>
+#include <boost/numeric/mtl/operation/resource.hpp>
+#include <boost/numeric/mtl/operation/lazy.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+namespace itl {
+
+/// Conjugate Gradients without preconditioning
+template < typename LinearOperator, typename HilbertSpaceX, typename HilbertSpaceB,
+	   typename Iteration >
+int cg(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
+       Iteration& iter)
+{
+    mtl::vampir_trace<7001> tracer;
+    using std::abs; using mtl::conj; using mtl::lazy;
+    typedef HilbertSpaceX Vector;
+    typedef typename mtl::Collection<HilbertSpaceX>::value_type Scalar;
+    typedef typename Iteration::real                            Real;
+
+    Scalar rho(0), rho_1(0), alpha(0), alpha_1(0);
+    Vector p(resource(x)), q(resource(x)), r(resource(x)), z(resource(x));
+
+    r = b - A*x;
+    rho = dot(r, r);
+    while (! iter.finished(Real(sqrt(abs(rho))))) {
+	++iter;
+	if (iter.first())
+	    p = r;
+	else
+	    p = r + (rho / rho_1) * p;
+
+	// q = A * p; alpha = rho / dot(p, q);
+	(lazy(q)= A * p) || (lazy(alpha_1)= lazy_dot(p, q));
+	alpha= rho / alpha_1;
+
+	x += alpha * p;
+	rho_1 = rho;
+	(lazy(r) -= alpha * q) || (lazy(rho) = lazy_unary_dot(r));
+    }
+
+    return iter;
+}
+
+/// Conjugate Gradients
+template < typename LinearOperator, typename HilbertSpaceX, typename HilbertSpaceB,
+	   typename Preconditioner, typename Iteration >
+int cg(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
+       const Preconditioner& L, Iteration& iter)
+{
+    using pc::is_identity;
+    if (is_identity(L))
+	return cg(A, x, b, iter);
+
+    mtl::vampir_trace<7002> tracer;
+    using std::abs; using mtl::conj; using mtl::lazy;
+    typedef HilbertSpaceX Vector;
+    typedef typename mtl::Collection<HilbertSpaceX>::value_type Scalar;
+    typedef typename Iteration::real                            Real;
+
+    Scalar rho(0), rho_1(0), rr, alpha(0), alpha_1;
+    Vector p(resource(x)), q(resource(x)), r(resource(x)), z(resource(x));
+
+    r = b - A*x;
+    rr = dot(r, r);
+    while (! iter.finished(Real(sqrt(abs(rr))))) {
+	++iter;
+	(lazy(z)= solve(L, r)) || (lazy(rho)= lazy_dot(r, z));
+
+	if (iter.first())
+	    p = z;
+	else
+	    p = z + (rho / rho_1) * p;
+
+	(lazy(q)= A * p) || (lazy(alpha_1)= lazy_dot(p, q));
+	alpha= rho / alpha_1;
+
+	x += alpha * p;
+	rho_1 = rho;
+	(lazy(r) -= alpha * q) || (lazy(rr) = lazy_unary_dot(r));
+    }
+    return iter;
+}
+
+/// Conjugate Gradients with ignored right preconditioner to unify interface
+template < typename LinearOperator, typename HilbertSpaceX, typename HilbertSpaceB,
+	   typename Preconditioner, typename RightPreconditioner, typename Iteration >
+int cg(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
+       const Preconditioner& L, const RightPreconditioner&, Iteration& iter)
+{
+    return cg(A, x, b, L, iter);
+}
+
+/// Solver class for CG method; right preconditioner ignored (prints warning if not identity)
+template < typename LinearOperator, typename Preconditioner,
+	   typename RightPreconditioner>
+class cg_solver
+{
+  public:
+    /// Construct solver from a linear operator; generate (left) preconditioner from it
+    explicit cg_solver(const LinearOperator& A) : A(A), L(A)
+    {
+	if (!pc::static_is_identity<RightPreconditioner>::value)
+	    std::cerr << "Right Preconditioner ignored!" << std::endl;
+    }
+
+    /// Construct solver from a linear operator and (left) preconditioner
+    cg_solver(const LinearOperator& A, const Preconditioner& L) : A(A), L(L)
+    {
+	if (!pc::static_is_identity<RightPreconditioner>::value)
+	    std::cerr << "Right Preconditioner ignored!" << std::endl;
+    }
+
+    /// Solve linear system approximately as specified by \p iter
+    template < typename HilbertSpaceB, typename HilbertSpaceX, typename Iteration >
+    int solve(const HilbertSpaceB& b, HilbertSpaceX& x, Iteration& iter) const
+    {
+	return cg(A, x, b, L, iter);
+    }
+
+    /// Perform one CG iteration on linear system
+    template < typename HilbertSpaceB, typename HilbertSpaceX >
+    int solve(const HilbertSpaceB& b, HilbertSpaceX& x) const
+    {
+	itl::basic_iteration<double> iter(x, 1, 0, 0);
+	return solve(b, x, iter);
+    }
+
+  private:
+    const LinearOperator& A;
+    Preconditioner        L;
+};
+
+
+} // namespace itl
+
+#endif // ITL_CG_INCLUDE

AMDiS/lib/mtl4/boost/numeric/itl/krylov/cgs.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef ITL_CGS_INCLUDE
+#define ITL_CGS_INCLUDE
+
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/operation/resource.hpp>
+#include <boost/numeric/mtl/operation/dot.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+namespace itl {
+
+/// Conjugate Gradient Squared
+template < typename LinearOperator, typename Vector,
+	   typename Preconditioner, typename Iteration >
+int cgs(const LinearOperator &A, Vector &x, const Vector &b,
+	const Preconditioner &M, Iteration& iter)
+{
+    mtl::vampir_trace<7007> tracer;
+    typedef typename mtl::Collection<Vector>::value_type Scalar;
+    Scalar     rho_1(0), rho_2(0), alpha(0), beta(0);
+    Vector     p(resource(x)), phat(resource(x)), q(resource(x)), qhat(resource(x)), vhat(resource(x)),
+	       u(resource(x)), uhat(resource(x)), r(b - A * x), rtilde= r;
+
+    while (! iter.finished(r)) {
+	++iter;
+	rho_1= dot(rtilde, r);
+
+	if (rho_1 == 0.) iter.fail(2, "cgs breakdown");
+
+	if (iter.first())
+	    p= u= r;
+	else {
+	    beta = rho_1 / rho_2;
+	    u= r + beta * q;
+	    p= u + beta * (q + beta * p);
+	}
+
+        vhat= A * Vector(solve(M, p));
+	alpha = rho_1 / dot(rtilde, vhat);
+	q= u - alpha * vhat;
+
+	u+= q;
+	uhat= solve(M, u);
+
+	x+= alpha * uhat;
+	qhat= A * uhat;
+	r-= alpha * qhat;
+
+	rho_2= rho_1;
+    }
+    return iter;
+}
+
+/// Solver class for CGS method; right preconditioner ignored (prints warning if not identity)
+template < typename LinearOperator, typename Preconditioner= pc::identity<LinearOperator>,
+	   typename RightPreconditioner= pc::identity<LinearOperator> >
+class cgs_solver
+{
+  public:
+    /// Construct solver from a linear operator; generate (left) preconditioner from it
+    explicit cgs_solver(const LinearOperator& A) : A(A), L(A)
+    {
+	if (!pc::static_is_identity<RightPreconditioner>::value)
+	    std::cerr << "Right Preconditioner ignored!" << std::endl;
+    }
+
+    /// Construct solver from a linear operator and (left) preconditioner
+    cgs_solver(const LinearOperator& A, const Preconditioner& L) : A(A), L(L)
+    {
+	if (!pc::static_is_identity<RightPreconditioner>::value)
+	    std::cerr << "Right Preconditioner ignored!" << std::endl;
+    }
+
+    /// Solve linear system approximately as specified by \p iter
+    template < typename HilbertSpaceB, typename HilbertSpaceX, typename Iteration >
+    int solve(const HilbertSpaceB& b, HilbertSpaceX& x, Iteration& iter) const
+    {
+	return cgs(A, x, b, L, iter);
+    }
+
+    /// Perform one CGS iteration on linear system
+    template < typename HilbertSpaceB, typename HilbertSpaceX >
+    int solve(const HilbertSpaceB& b, HilbertSpaceX& x) const
+    {
+	itl::basic_iteration<double> iter(x, 1, 0, 0);
+	return solve(b, x, iter);
+    }
+
+  private:
+    const LinearOperator& A;
+    Preconditioner        L;
+};
+
+} // namespace itl
+
+#endif // ITL_CGS_INCLUDE
+
+
+

AMDiS/lib/mtl4/boost/numeric/itl/krylov/fsm.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef ITL_FSM_INCLUDE
+#define ITL_FSM_INCLUDE
+
+#include <boost/numeric/mtl/operation/two_norm.hpp>
+
+namespace itl {
+
+/// Folded spectrum method
+/** Computed and named as in http://en.wikipedia.org/wiki/Folded_spectrum_method **/
+template < typename LinearOperator, typename VectorSpace, typename EigenValue,
+	   typename Damping, typename Iteration >
+int fsm(const LinearOperator& H, VectorSpace& phi, EigenValue eps, Damping alpha, Iteration& iter)
+{
+    VectorSpace v1(H * phi - eps * phi);
+    for (; !iter.finished(v1); ++iter) {
+	VectorSpace v2(H * v1 - eps * v1);
+	phi-= alpha * v2;
+	phi/= two_norm(phi);
+	v1= H * phi - eps * phi;
+    }
+    return iter;
+}
+
+
+} // namespace itl
+
+#endif // ITL_FSM_INCLUDE

AMDiS/lib/mtl4/boost/numeric/itl/krylov/gmres.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.
+
+// Written by Cornelius Steinhardt
+
+
+#ifndef ITL_GMRES_INCLUDE
+#define ITL_GMRES_INCLUDE
+
+#include <algorithm>
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/vector/dense_vector.hpp>
+#include <boost/numeric/mtl/matrix/dense2D.hpp>
+#include <boost/numeric/mtl/matrix/multi_vector.hpp>
+#include <boost/numeric/mtl/operation/givens.hpp>
+#include <boost/numeric/mtl/operation/two_norm.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/mtl/utility/irange.hpp>
+
+namespace itl {
+
+/// Generalized Minimal Residual method (without restart)
+/** It computes at most kmax_in iterations (or size(x) depending on what is smaller)
+    regardless on whether the termination criterion is reached or not.   **/
+template < typename Matrix, typename Vector, typename LeftPreconditioner, typename RightPreconditioner, typename Iteration >
+int gmres_full(const Matrix &A, Vector &x, const Vector &b,
+               LeftPreconditioner &L, RightPreconditioner &R, Iteration& iter)
+{
+    using mtl::size; using mtl::irange; using mtl::iall; using std::abs; using std::sqrt;
+    typedef typename mtl::Collection<Vector>::value_type Scalar;
+    typedef typename mtl::Collection<Vector>::size_type  Size;
+
+    if (size(b) == 0) throw mtl::logic_error("empty rhs vector");
+
+    const Scalar                zero= math::zero(Scalar());
+    Scalar                      rho, nu, hr;
+    Size                        k, kmax(std::min(size(x), Size(iter.max_iterations() - iter.iterations())));
+    Vector                      r0(b - A *x), r(solve(L,r0)), va(resource(x)), va0(resource(x)), va00(resource(x));
+    mtl::matrix::multi_vector<Vector>   V(Vector(resource(x), zero), kmax+1);
+    mtl::vector::dense_vector<Scalar>   s(kmax+1, zero), c(kmax+1, zero), g(kmax+1, zero), y(kmax, zero);  // replicated in distributed solvers
+    mtl::matrix::dense2D<Scalar>        H(kmax+1, kmax);                                             // dito
+    H= 0;
+
+    rho= g[0]= two_norm(r);
+    if (iter.finished(rho))
+	return iter;
+    V.vector(0)= r / rho;
+    H= zero;
+
+    // GMRES iteration
+    for (k= 0; k < kmax ; ++k, ++iter) {
+        va0= A * Vector(solve(R, V.vector(k)));
+        V.vector(k+1)= va= solve(L,va0);
+	// orth(V, V[k+1], false);
+        // modified Gram Schmidt method
+        for (Size j= 0; j < k+1; j++) {
+	    H[j][k]= dot(V.vector(j), V.vector(k+1));
+	    V.vector(k+1)-= H[j][k] * V.vector(j);
+        }
+        H[k+1][k]= two_norm(V.vector(k+1));
+        //reorthogonalize
+        for(Size j= 0; j < k+1; j++) {
+	    hr= dot(V.vector(k+1), V.vector(j));
+            H[j][k]+= hr;
+            V.vector(k+1)-= hr * V.vector(j);
+        }
+        H[k+1][k]= two_norm(V.vector(k+1));
+	if (H[k+1][k] != zero)                // watch for breakdown
+            V.vector(k+1)*= 1. / H[k+1][k];
+
+        // k Given's rotations
+	for(Size i= 0; i < k; i++)
+	    mtl::matrix::givens<mtl::matrix::dense2D<Scalar> >(H, H[i][k-1], H[i+1][k-1]).trafo(i);
+
+       nu= sqrt(H[k][k]*H[k][k]+H[k+1][k]*H[k+1][k]);
+       if(nu != zero){
+            c[k]=  H[k][k]/nu;
+            s[k]= -H[k+1][k]/nu;
+            H[k][k]=c[k]*H[k][k]-s[k]*H[k+1][k];
+            H[k+1][k]=0;
+ 	    mtl::vector::givens<mtl::vector::dense_vector<Scalar> >(g, c[k], s[k]).trafo(k);
+        }
+	rho= abs(g[k+1]);
+    }
+
+    //reduce k, to get regular matrix
+    while (k > 0 && abs(g[k-1]<= iter.atol())) k--;
+
+    // iteration is finished -> compute x: solve H*y=g as far as rank of H allows
+    irange                  range(k);
+    for (; !range.empty(); --range) {
+	try {
+	    y[range]= lu_solve(H[range][range], g[range]);
+	} catch (mtl::matrix_singular) { continue; } // if singular then try with sub-matrix
+	break;
+    }
+
+    if (range.finish() < k)
+  	std::cerr << "GMRES orhogonalized with " << k << " vectors but matrix singular, can only use "
+		  << range.finish() << " vectors!\n";
+    if (range.empty())
+        return iter.fail(2, "GMRES did not find any direction to correct x");
+    x+= Vector(solve(R, Vector(V.vector(range)*y[range])));
+
+    r= b - A*x;
+    return iter.terminate(r);
+}
+
+/// Generalized Minimal Residual method with restart
+template < typename Matrix, typename Vector, typename LeftPreconditioner,
+           typename RightPreconditioner, typename Iteration >
+int gmres(const Matrix &A, Vector &x, const Vector &b,
+          LeftPreconditioner &L, RightPreconditioner &R,
+	  Iteration& iter, typename mtl::Collection<Vector>::size_type restart)
+{
+     do {
+	 Iteration inner(iter);
+	 inner.set_max_iterations(std::min(int(iter.iterations()+restart), iter.max_iterations()));
+	 inner.suppress_resume(true);
+	 gmres_full(A, x, b, L, R, inner);
+	 iter.update_progress(inner);
+     } while (!iter.finished());
+
+     return iter;
+}
+
+/// Solver class for GMRES; right preconditioner ignored (prints warning if not identity)
+template < typename LinearOperator, typename Preconditioner= pc::identity<LinearOperator>,
+	   typename RightPreconditioner= pc::identity<LinearOperator> >
+class gmres_solver
+{
+  public:
+    /// Construct solver from a linear operator; generate (left) preconditioner from it
+    explicit gmres_solver(const LinearOperator& A, size_t restart= 8)
+      : A(A), restart(restart), L(A), R(A) {}
+
+    /// Construct solver from a linear operator and left preconditioner
+    gmres_solver(const LinearOperator& A, size_t restart, const Preconditioner& L)
+      : A(A), restart(restart), L(L), R(A) {}
+
+    /// Construct solver from a linear operator and left preconditioner
+    gmres_solver(const LinearOperator& A, size_t restart, const Preconditioner& L, const RightPreconditioner& R)
+      : A(A), restart(restart), L(L), R(R) {}
+
+    /// Solve linear system approximately as specified by \p iter
+    template < typename HilbertSpaceB, typename HilbertSpaceX, typename Iteration >
+    int solve(const HilbertSpaceB& b, HilbertSpaceX& x, Iteration& iter) const
+    {
+	return gmres(A, x, b, L, R, iter, restart);
+    }
+
+    /// Perform one GMRES iteration on linear system
+    template < typename HilbertSpaceB, typename HilbertSpaceX >
+    int solve(const HilbertSpaceB& b, HilbertSpaceX& x) const
+    {
+	itl::basic_iteration<double> iter(x, 1, 0, 0);
+	return solve(b, x, iter);
+    }
+
+  private:
+    const LinearOperator& A;
+    size_t                restart;
+    Preconditioner        L;
+    RightPreconditioner   R;
+};
+
+
+} // namespace itl
+
+#endif // ITL_GMRES_INCLUDE
+
+

AMDiS/lib/mtl4/boost/numeric/itl/krylov/idr_s.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.
+
+// Peter Sonneveld and Martin B. van Gijzen, IDR(s): a family of simple and fast algorithms for solving large nonsymmetric linear systems.
+// SIAM J. Sci. Comput. Vol. 31, No. 2, pp. 1035-1062 (2008). (copyright SIAM)
+
+#ifndef ITL_IDR_S_INCLUDE
+#define ITL_IDR_S_INCLUDE
+
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/vector/dense_vector.hpp>
+#include <boost/numeric/mtl/operation/random.hpp>
+#include <boost/numeric/mtl/operation/orth.hpp>
+#include <boost/numeric/mtl/operation/resource.hpp>
+#include <boost/numeric/mtl/matrix/strict_upper.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/mtl/utility/irange.hpp>
+#include <boost/numeric/linear_algebra/identity.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+namespace itl {
+
+/// Induced Dimension Reduction on s dimensions (IDR(s))
+template < typename LinearOperator, typename Vector,
+	   typename LeftPreconditioner, typename RightPreconditioner,
+	   typename Iteration >
+int idr_s(const LinearOperator &A, Vector &x, const Vector &b,
+	  const LeftPreconditioner &, const RightPreconditioner &,
+	  Iteration& iter, size_t s)
+{
+    mtl::vampir_trace<7010> tracer;
+    using mtl::size; using mtl::iall; using mtl::matrix::strict_upper;
+    typedef typename mtl::Collection<Vector>::value_type Scalar;
+    typedef typename mtl::Collection<Vector>::size_type  Size;
+
+    if (size(b) == 0) throw mtl::logic_error("empty rhs vector");
+    if (s < 1) s= 1;
+
+    const Scalar                zero= math::zero(Scalar());
+    Scalar                      omega(zero);
+    Vector                      x0(x), y(resource(x)), v(resource(x)), t(resource(x)), q(resource(x)), r(b - A * x);
+    mtl::matrix::multi_vector<Vector>   dR(Vector(resource(x), zero), s), dX(Vector(resource(x), zero), s), P(Vector(resource(x), zero), s);
+    mtl::vector::dense_vector<Scalar>   m(s), c(s), dm(s);   // replicated in distributed solvers
+    mtl::matrix::dense2D<Scalar>        M(s, s);             // dito
+
+    random(P);
+    P.vector(0)= r;
+    orth(P);
+
+    for (size_t k= 0; k < s; k++) {
+	v= A * r;
+	omega= dot(v, r) / dot(v, v);
+	dX.vector(k)= omega * r;
+	dR.vector(k)= -omega * v;
+	x+= dX.vector(k);
+	r+= dR.vector(k);
+	if ((++iter).finished(r)) return iter;
+	M[iall][k]= trans(P) * dR.vector(k);
+    }
+
+    Size oldest= 0;
+    m= trans(P) * r;
+
+    while (! iter.finished(r)) {
+
+	for (size_t k= 0; k < s; k++) {
+	    c= lu_solve(M, m);
+	    q= dR * -c;
+	    v= r + q;
+	    if (k == 0) {
+		t= A * v;
+		omega= dot(t, v) / dot(t, t);
+		dR.vector(oldest)= q - omega * t;
+		dX.vector(oldest)= omega * v - dX * c;
+	    } else {
+		dX.vector(oldest)= omega * v - dX * c;
+		dR.vector(oldest)= A * -dX.vector(oldest);
+	    }
+	    r+= dR.vector(oldest);
+	    x+= dX.vector(oldest);
+
+	    if ((++iter).finished(r))
+		return iter;
+
+	    dm= trans(P) * dR.vector(oldest);
+	    M[iall][oldest]= dm;
+	    m+= dm;
+	    oldest= (oldest + 1) % s;
+	}
+    }
+    return iter;
+}
+
+/// Solver class for IDR(s) method; right preconditioner ignored (prints warning if not identity)
+template < typename LinearOperator, typename Preconditioner= pc::identity<LinearOperator>,
+	   typename RightPreconditioner= pc::identity<LinearOperator> >
+class idr_s_solver
+{
+  public:
+    /// Construct solver from a linear operator; generate (left) preconditioner from it
+    explicit idr_s_solver(const LinearOperator& A, size_t s= 8) : A(A), s(s), L(A), R(A) {}
+
+    /// Construct solver from a linear operator and left preconditioner
+    idr_s_solver(const LinearOperator& A, size_t s, const Preconditioner& L) : A(A), s(s), L(L), R(A) {}
+
+    /// Construct solver from a linear operator and left preconditioner
+    idr_s_solver(const LinearOperator& A, size_t s, const Preconditioner& L, const RightPreconditioner& R)
+      : A(A), s(s), L(L), R(R) {}
+
+    /// Solve linear system approximately as specified by \p iter
+    template < typename HilbertSpaceB, typename HilbertSpaceX, typename Iteration >
+    int solve(const HilbertSpaceB& b, HilbertSpaceX& x, Iteration& iter) const
+    {
+	return idr_s(A, x, b, L, R, iter, s);
+    }
+
+    /// Perform one IDR(s) iteration on linear system
+    template < typename HilbertSpaceB, typename HilbertSpaceX >
+    int solve(const HilbertSpaceB& b, HilbertSpaceX& x) const
+    {
+	itl::basic_iteration<double> iter(x, 1, 0, 0);
+	return solve(b, x, iter);
+    }
+
+  private:
+    const LinearOperator& A;
+    size_t                s;
+    Preconditioner        L;
+    RightPreconditioner   R;
+};
+
+
+} // namespace itl
+
+#endif // ITL_IDR_S_INCLUDE

AMDiS/lib/mtl4/boost/numeric/itl/krylov/qmr.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.
+
+// Written by Cornelius Steinhardt
+
+
+#ifndef ITL_QMR_INCLUDE
+#define ITL_QMR_INCLUDE
+
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/operation/trans.hpp>
+#include <boost/numeric/mtl/operation/resource.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+namespace itl {
+
+/// Quasi-Minimal Residual method
+template < typename Matrix, typename Vector,typename LeftPreconditioner,
+	   typename RightPreconditioner, typename Iteration >
+int qmr(const Matrix& A, Vector& x, const Vector& b, LeftPreconditioner& L,
+	const RightPreconditioner& R, Iteration& iter)
+{
+    mtl::vampir_trace<7008> tracer;
+    using mtl::size;
+    typedef typename mtl::Collection<Vector>::value_type Scalar;
+    if (size(b) == 0) throw mtl::logic_error("empty rhs vector");
+
+    const Scalar         zero= math::zero(Scalar()), one= math::one(Scalar());
+    Scalar               beta, gamma(one), gamma_1, delta, eta(-one), ep(one), rho_1, theta(zero), theta_1;
+    Vector               r(b - A * x), v_tld(r), y(solve(L, v_tld)), w_tld(r), z(adjoint_solve(R,w_tld)), v(resource(x)), w(resource(x)),
+                         y_tld(resource(x)), z_tld(resource(x)), p(resource(x)), q(resource(x)), p_tld(resource(x)), d(resource(x)), s(resource(x));
+
+    if (iter.finished(r))
+	return iter;
+
+    Scalar rho = two_norm(y), xi = two_norm(z);
+    while(! iter.finished(rho)) {
+	++iter;
+        if (rho == zero)
+	    return iter.fail(1, "qmr breakdown #1, rho=0");
+        if (xi == zero)
+            return iter.fail(2, "qmr breakdown #2, xi=0");
+
+        v= v_tld / rho;
+        y/= rho;
+        w= w_tld / xi;
+        z/= xi;
+
+        delta = dot(z,y);
+        if (delta == zero)
+            return iter.fail(3, "qmr breakdown, delta=0 #3");
+
+        y_tld = solve(R,y);
+        z_tld = adjoint_solve(L,z);
+
+	if (iter.first()) {
+            p = y_tld;
+            q = z_tld;
+	} else {
+            p = y_tld - ((xi * delta) / ep) * p;
+            q = z_tld - ((rho* delta) / ep) * q;
+        }
+
+        p_tld = A * p;
+        ep = dot(q, p_tld);
+        if (ep == zero)
+            return iter.fail(4, "qmr breakdown ep=0 #4");
+        beta= ep / delta;
+        if (beta == zero)
+            return iter.fail(5, "qmr breakdown beta=0 #5");
+        v_tld = p_tld - beta * v;
+        y = solve(L,v_tld);
+        rho_1 = rho;
+	rho = two_norm(y);
+        w_tld= trans(A)*q  - beta*w;
+        z = adjoint_solve(R, w_tld);
+        xi = two_norm(z);
+        gamma_1 = gamma;
+        theta_1 = theta;
+        theta = rho / (gamma_1 * beta);
+        gamma = one / (sqrt(one + theta * theta));
+
+        if (gamma == zero)
+            return iter.fail(6, "qmr breakdown gamma=0 #6");
+
+        eta= -eta * rho_1 * gamma * gamma / (beta * gamma_1 * gamma_1);
+	if (iter.first()) {
+           d= eta * p;
+	   s= eta * p_tld;
+	} else {
+            d= eta * p + (theta_1 * theta_1 * gamma * gamma) * d;
+            s= eta * p_tld + (theta_1 * theta_1 * gamma * gamma) * s;
+        }
+        x += d;
+        r -= s;
+    }
+    return iter;
+}
+
+/// Solver class for Quasi-minimal residual method; right preconditioner ignored (prints warning if not identity)
+template < typename LinearOperator, typename Preconditioner= pc::identity<LinearOperator>,
+	   typename RightPreconditioner= pc::identity<LinearOperator> >
+class qmr_solver
+{
+  public:
+    /// Construct solver from a linear operator; generate (left) preconditioner from it
+    explicit qmr_solver(const LinearOperator& A) : A(A), L(A), R(A) {}
+
+    /// Construct solver from a linear operator and left preconditioner
+    qmr_solver(const LinearOperator& A, const Preconditioner& L) : A(A), L(L), R(A) {}
+
+    /// Construct solver from a linear operator and left preconditioner
+    qmr_solver(const LinearOperator& A, const Preconditioner& L, const RightPreconditioner& R)
+      : A(A), L(L), R(R) {}
+
+    /// Solve linear system approximately as specified by \p iter
+    template < typename HilbertSpaceB, typename HilbertSpaceX, typename Iteration >
+    int solve(const HilbertSpaceB& b, HilbertSpaceX& x, Iteration& iter) const
+    {
+	return qmr(A, x, b, L, R, iter);
+    }
+
+    /// Perform one QMR iteration on linear system
+    template < typename HilbertSpaceB, typename HilbertSpaceX >
+    int solve(const HilbertSpaceB& b, HilbertSpaceX& x) const
+    {
+	itl::basic_iteration<double> iter(x, 1, 0, 0);
+	return solve(b, x, iter);
+    }
+
+  private:
+    const LinearOperator& A;
+    Preconditioner        L;
+    RightPreconditioner   R;
+};
+
+
+} // namespace itl
+
+#endif // ITL_QMR_INCLUDE
+

AMDiS/lib/mtl4/boost/numeric/itl/krylov/repeating_solver.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG, www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also tools/license/license.mtl.txt in the distribution.
+
+#ifndef ITL_REPEATING_SOLVER_INCLUDE
+#define ITL_REPEATING_SOLVER_INCLUDE
+
+#include <boost/mpl/if.hpp>
+#include <boost/static_assert.hpp>
+#include <boost/numeric/itl/itl_fwd.hpp>
+
+namespace itl {
+
+/// Class for calling \tparam N iterations of the given \tparam Solver
+/** If \tparam Stored is true then the \tparam Solver object is stored (i.e. possibly copied) here.
+    Otherwise it is only referred and passing temporary objects to the constructor will cause errors. **/
+template <typename Solver, unsigned N, bool Stored>
+class repeating_solver
+{
+    typedef typename boost::mpl::if_c<Stored, Solver, const Solver&>::type solver_type;
+  public:
+    explicit repeating_solver(const Solver& s) : s(s) {}
+
+    template <typename Matrix>
+    explicit repeating_solver(const Matrix& A) : s(A)
+    {
+	BOOST_STATIC_ASSERT((Stored)); // if matrix is passed class must own solver
+    }
+
+    /// Perform one GMRES iteration on linear system
+    template < typename VectorIn, typename VectorOut >
+    int solve(const VectorIn& b, VectorOut& x) const
+    {
+	itl::basic_iteration<double> iter(x, N, 0, 0);
+	return s.solve(b, x, iter);
+    }
+
+  private:
+    solver_type   s;
+};
+
+} // namespace itl
+
+#endif // ITL_REPEATING_SOLVER_INCLUDE

AMDiS/lib/mtl4/boost/numeric/itl/krylov/tfqmr.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.
+
+// Written by Cornelius Steinhardt
+
+
+#ifndef ITL_TFQMR_INCLUDE
+#define ITL_TFQMR_INCLUDE
+
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/linear_algebra/identity.hpp>
+#include <boost/numeric/linear_algebra/inverse.hpp>
+#include <boost/numeric/mtl/utility/irange.hpp>
+#include <boost/numeric/mtl/operation/resource.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+namespace itl {
+
+/// Transposed-free Quasi-minimal residual
+template < typename Matrix, typename Vector,
+	   typename LeftPreconditioner, typename RightPreconditioner, typename Iteration >
+int tfqmr(const Matrix &A, Vector &x, const Vector &b, const LeftPreconditioner &L,
+	  const RightPreconditioner &R, Iteration& iter)
+{
+    mtl::vampir_trace<7009> tracer;
+    using math::reciprocal; using mtl::size;
+    typedef typename mtl::Collection<Vector>::value_type Scalar;
+
+    if (size(b) == 0) throw mtl::logic_error("empty rhs vector");
+
+    const Scalar                zero= math::zero(Scalar()), one= math::one(Scalar());
+    Scalar                      theta(zero), eta(zero), tau, rho, rhon, sigma, alpha, beta, c;
+    Vector                      rt(b - A*Vector(solve(R, x))) /* shift x= R*x */, r(solve(L, rt)), u1(resource(x)), u2(resource(x)),
+                                y1(resource(x)), y2(resource(x)), w(resource(x)), d(resource(x), zero), v(resource(x));
+
+    if (iter.finished(rt))
+	return iter;
+    y1= w= r;
+    rt= A * Vector(solve(R, y1));
+    u1= v= solve(L,rt);
+    tau= two_norm(r);
+    rho= tau * tau;
+
+    // TFQMR iteration
+    while (! iter.finished(tau)) {
+	++iter;
+	sigma= dot(r,v);
+        if (sigma == zero)
+	    return iter.fail(1, "tfgmr breakdown, sigma=0 #1");
+        alpha= rho / sigma;
+
+        // inner loop
+        for(int j=1; j < 3; j++) {
+            if (j == 1) {
+                w-= alpha * u1;
+                d= y1+ (theta * theta * eta / alpha) * d;
+	    } else {
+                y2= y1 - alpha * v;
+                rt= A * Vector(solve(R, y2));
+                u2= solve(L, rt);
+                w-= alpha * u2;
+                d= y2 + (theta * theta * eta / alpha) * d;
+            }
+            theta= two_norm(w) / tau;
+            c= reciprocal(sqrt(one + theta*theta));
+            tau*= theta * c;
+            eta= c * c * alpha;
+            x+= eta * d;
+        } // end inner loop
+        if (rho == zero)
+            return iter.fail(1, "tfgmr breakdown, rho=0 #2");
+        rhon= dot(r,w);
+        beta= rhon/rho;
+        rho= rhon;
+        y1= w + beta*y2;
+        rt= A * Vector(solve(R, y1));
+        u1= solve(L, rt);
+        v= u1 + beta*(u2 + beta*v);
+        rt= A * x - b;
+    }
+    //shift back
+    x= solve(R, x);
+    return iter;
+}
+
+/// Solver class for Transposed-free quasi-minimal residual method; right preconditioner ignored (prints warning if not identity)
+template < typename LinearOperator, typename Preconditioner= pc::identity<LinearOperator>,
+	   typename RightPreconditioner= pc::identity<LinearOperator> >
+class tfqmr_solver
+{
+  public:
+    /// Construct solver from a linear operator; generate (left) preconditioner from it
+    explicit tfqmr_solver(const LinearOperator& A) : A(A), L(A), R(A) {}
+
+    /// Construct solver from a linear operator and left preconditioner
+    tfqmr_solver(const LinearOperator& A, const Preconditioner& L) : A(A), L(L), R(A) {}
+
+    /// Construct solver from a linear operator and left preconditioner
+    tfqmr_solver(const LinearOperator& A, const Preconditioner& L, const RightPreconditioner& R)
+      : A(A), L(L), R(R) {}
+
+    /// Solve linear system approximately as specified by \p iter
+    template < typename HilbertSpaceB, typename HilbertSpaceX, typename Iteration >
+    int solve(const HilbertSpaceB& b, HilbertSpaceX& x, Iteration& iter) const
+    {
+	return tfqmr(A, x, b, L, R, iter);
+    }
+
+    /// Perform one TFQMR iteration on linear system
+    template < typename HilbertSpaceB, typename HilbertSpaceX >
+    int solve(const HilbertSpaceB& b, HilbertSpaceX& x) const
+    {
+	itl::basic_iteration<double> iter(x, 1, 0, 0);
+	return solve(b, x, iter);
+    }
+
+  private:
+    const LinearOperator& A;
+    Preconditioner        L;
+    RightPreconditioner   R;
+};
+
+} // namespace itl
+
+#endif // ITL_TFQMR_INCLUDE

AMDiS/lib/mtl4/boost/numeric/itl/minimization/quasi_newton.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef ITL_QUASI_NEWTON_INCLUDE
+#define ITL_QUASI_NEWTON_INCLUDE
+
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/matrix/dense2D.hpp>
+#include <boost/numeric/mtl/matrix/operators.hpp>
+#include <boost/numeric/mtl/operation/operators.hpp>
+#include <boost/numeric/mtl/operation/trans.hpp>
+#include <boost/numeric/mtl/utility/gradient.hpp>
+
+// #include <iostream>
+
+namespace itl {
+
+/// Quasi-Newton method
+template <typename Matrix, typename Vector, typename F, typename Grad,
+	  typename Step, typename Update, typename Iter>
+Vector quasi_newton(Vector& x, F f, Grad grad_f, Step step, Update update, Iter& iter)
+{
+    typedef typename mtl::Collection<Vector>::value_type value_type;
+    Vector         d, y, x_k, s;
+    Matrix         H(size(x), size(x));
+
+    H= 1;
+    for (; !iter.finished(two_norm(grad_f(x))); ++iter) {
+	d= H * -grad_f(x);                                                  // std::cout << "d is " << d << '\n';
+	value_type alpha= step(x, d, f, grad_f); assert(alpha == alpha);
+	x_k= x + alpha * d;                                                 // std::cout << "x_k is " << x_k << '\n';
+	s= alpha * d;                                                       // std::cout << "alpha is " << alpha << '\n';
+	y= grad_f(x_k) - grad_f(x);
+	update(H, y, s);
+	x= x_k;
+    }
+    return x;
+}
+
+/// Quasi-Newton method
+template <typename Vector, typename F, typename Grad, typename Step, typename Update, typename Iter>
+Vector inline quasi_newton(Vector& x, F f, Grad grad_f, Step step, Update update, Iter& iter)
+{
+    typedef typename mtl::traits::gradient<Vector>::type hessian_type;
+    // typedef typename mtl::Collection<Vector>::value_type value_type;
+    return quasi_newton<hessian_type>(x, f, grad_f, step, update, iter);
+}
+
+
+} // namespace itl
+
+#endif // ITL_QUASI_NEWTON_INCLUDE

AMDiS/lib/mtl4/boost/numeric/itl/pc/binary_heap.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+//
+// Algorithm inspired by Nick Vannieuwenhoven, rewritten by Cornelius Steinhardt
+//
+//   Created on: Jan 10, 2010
+//       Author: heazk
+
+
+#ifndef MTL_BINARY_HEAP_INCLUDE
+#define MTL_BINARY_HEAP_INCLUDE
+
+namespace utils {
+
+/**
+ * An intrusive binary min-heap.
+ *
+ * KeyType: 		The type of the keys. Assumed to be a signed integer type.
+ * KeyComparator:	A functor to compare values of the key type. Should be a
+ * 					total order.
+ */
+template<
+	class DirectAccessIterator,
+	class KeyType,
+	class KeyComparator,
+	class ValueType,
+	class GetKey,
+	class GetParent,
+	class GetLeft,
+	class GetRight
+>
+class binary_heap {
+
+
+/*******************************************************************************
+ * Type Definitions
+ ******************************************************************************/
+
+
+public:
+	/**
+	 * The value type.
+	 */
+	typedef ValueType value_type;
+
+	/**
+	 * The key type.
+	 */
+	typedef KeyType key_type;
+
+	/**
+	 * The type of this heap.
+	 */
+	typedef binary_heap<
+		DirectAccessIterator,
+		KeyType,
+		KeyComparator,
+		ValueType,
+		GetKey,
+		GetParent,
+		GetLeft,
+		GetRight
+	> heap_type;
+
+
+/*******************************************************************************
+ * Constructors and Destructors
+ ******************************************************************************/
+
+
+public:
+	binary_heap(
+		DirectAccessIterator values_begin,
+		DirectAccessIterator values_end,
+		KeyComparator compare_keys,
+		GetKey get_key,
+		GetParent get_parent,
+		GetLeft get_left_child,
+		GetRight get_right_child
+	) :
+		m_root(values_begin),
+		m_compare_keys(compare_keys),
+		m_key(get_key),
+		m_parent(get_parent),
+		m_left_child(get_left_child),
+		m_right_child(get_right_child)
+	{
+		build_heap(values_begin, values_end);
+	}
+
+public:
+	binary_heap(
+		KeyComparator compare_keys,
+		GetKey get_key,
+		GetParent get_parent,
+		GetLeft get_left_child,
+		GetRight get_right_child
+	) :
+		m_root(0),
+		m_last(0),
+		m_compare_keys(compare_keys),
+		m_key(get_key),
+		m_parent(get_parent),
+		m_left_child(get_left_child),
+		m_right_child(get_right_child)
+	{
+	}
+
+	/**
+	 * Default destructor.
+	 */
+public:
+	~binary_heap() {
+		// Nothing here.
+	}
+
+	/**
+	 * Copying of heaps is disallowed.
+	 */
+private:
+	binary_heap(const heap_type&);
+	void operator=(const heap_type&);
+
+
+/*******************************************************************************
+ * Heap Construction and Maintaining
+ ******************************************************************************/
+
+
+	/**
+	 * Builds the min-heap.
+	 */
+private:
+	void build_heap(
+		DirectAccessIterator values_begin,
+		DirectAccessIterator values_end
+	) {
+		std::cout << "building heap ..." << std::endl;
+		// Construct heap connections.
+		const int size = values_end - values_begin;
+		m_parent(values_begin) = 0;
+		m_left_child(values_begin) = 	1 < size ? values_begin+1 : 0;
+		m_right_child(values_begin) = 	2 < size ? values_begin+2 : 0;
+		for(int i = 1; i < size; ++i) {
+			m_parent(values_begin+i) = values_begin+(((i+1)/2) - 1);
+			m_left_child(values_begin+i) =
+				2*(i+1)-1 < size ? values_begin+(2*(i+1)-1) : 0;
+			m_right_child(values_begin+i) =
+				2*(i+1) < size ? values_begin+(2*(i+1)) : 0;
+		}
+
+		m_root = size > 0 ? values_begin : 0;
+		m_last = size > 0 ? values_end-1 : 0;
+
+		build_heap_rec(m_root);
+
+	}
+
+	/**
+	 * Applies the min-heapify procedure at each node in postfix style.
+	 */
+private:
+	void build_heap_rec(value_type node) {
+		if(node) {
+//			dump(node);
+			build_heap_rec(m_left_child(node));
+			build_heap_rec(m_right_child(node));
+			min_heapify(node);
+		}
+	}
+
+	/**
+	 * The min-heapify operation to restore the min heap property at a given
+	 * node.
+	 */
+private:
+	void min_heapify(value_type node) {
+
+//		std::cout << "heapify: " << m_key(node) << std::endl;
+
+		value_type seek = node;
+		value_type smallest = seek;
+		do {
+			value_type left_child = m_left_child(seek);
+			value_type right_child = m_right_child(seek);
+			smallest = seek;
+
+			if(
+					left_child &&
+					m_compare_keys(left_child, seek)
+			) {
+				smallest = left_child;
+			}
+			if(
+					right_child &&
+					m_compare_keys(right_child, smallest)
+			) {
+				smallest = right_child;
+			}
+			if(smallest != seek) {
+				exchange(smallest, seek);
+			}
+		} while(smallest != seek);
+
+	}
+
+	/**
+	 * Exchanges two nodes.
+	 */
+private:
+	inline void exchange(value_type fst, value_type snd) {
+		if(!fst && !snd) {
+			return;
+		}
+		if(!fst || !snd) {
+			assert(false);
+		}
+		if(fst == snd) {
+			return;
+		}
+
+		// Update root and last pointers.
+		if(fst == m_root) {
+			m_root = snd;
+		} else if(snd == m_root) {
+			m_root = fst;
+		}
+		if(fst == m_last) {
+			m_last = snd;
+		} else if(snd == m_last) {
+			m_last = fst;
+		}
+
+//		std::cout << "fst pre" << std::endl;
+//		dump(fst);
+//		std::cout << "snd pre" << std::endl;
+//		dump(snd);
+
+		// Canonize the situation when fst and snd are directly connected such
+		// that fst is always the higher node.
+		if( m_parent(fst) == snd ) {
+			value_type tmp = fst;
+			fst = snd;
+			snd = tmp;
+		}
+
+		// Set the new parent pointers for the children.
+		value_type left_child = m_left_child(fst);
+		value_type right_child = m_right_child(fst);
+		if(left_child) {
+			m_parent(left_child) = snd;
+		}
+		if(right_child) {
+			m_parent(right_child) = snd;
+		}
+		left_child = m_left_child(snd);
+		right_child = m_right_child(snd);
+		if(left_child) {
+			m_parent(left_child) = fst;
+		}
+		if(right_child) {
+			m_parent(right_child) = fst;
+		}
+
+		// Set the new children pointers.
+		const value_type fst_left_child = m_left_child(fst);
+		const value_type fst_right_child = m_right_child(fst);
+		m_left_child(fst) = m_left_child(snd);
+		m_right_child(fst) = m_right_child(snd);
+		m_left_child(snd) = fst_left_child;
+		m_right_child(snd) = fst_right_child;
+
+		// Set the new children pointers of the parents.
+		if(m_parent(fst)) {
+			if( m_left_child(m_parent(fst)) == fst ) {
+				m_left_child(m_parent(fst)) = snd;
+			} else {
+				m_right_child(m_parent(fst)) = snd;
+			}
+		}
+		if(m_parent(snd)) {
+			if( m_left_child(m_parent(snd)) == snd ) {
+				m_left_child(m_parent(snd)) = fst;
+			} else {
+				m_right_child(m_parent(snd)) = fst;
+			}
+		}
+
+		// Set the new parent pointers.
+		const value_type fst_parent = m_parent(fst);
+		m_parent(fst) = m_parent(snd);
+		m_parent(snd) = fst_parent;
+
+//		std::cout << "fst post" << std::endl;
+//		dump(snd);
+//		std::cout << "snd post" << std::endl;
+//		dump(fst);
+
+//		if( m_parent(fst) == snd || m_parent(snd) == fst) {
+//			// Canonize the situation to this whereby fst is the higher node.
+//			if( m_parent(fst) == snd ) {
+//				value_type tmp = fst;
+//				fst = snd;
+//				snd = tmp;
+//			}
+//
+//			// Update root and last pointers.
+//			if(fst == m_root) {
+//				m_root = snd;
+//			}
+//			if(snd == m_last) {
+//				m_last = fst;
+//			}
+//
+//			// Get a hold on every value we need.
+//			value_type fst_parent = m_parent(fst);
+//			value_type fst_left_child = m_left_child(fst);
+//			value_type fst_right_child = m_right_child(fst);
+//			value_type snd_left_child = m_left_child(snd);
+//			value_type snd_right_child = m_right_child(snd);
+//
+//			// Set up fst's new pointers.
+//			m_left_child(fst) = snd_left_child;
+//			m_right_child(fst) = snd_right_child;
+//			m_parent(snd_left_child) = fst;
+//			m_parent(snd_right_child) = fst;
+//			m_parent(fst) = snd;
+//
+//			// Set up snd's new pointers.
+//			m_parent(snd) = fst_parent;
+//			if(fst_parent) {
+//				if( m_left_child(fst_parent) == fst ) {
+//					m_left_child(fst_parent) = snd;
+//				} else {
+//					m_right_child(fst_parent) = snd;
+//				}
+//			}
+//
+//			if(fst_left_child == snd) {
+//				m_parent(fst_right_child) = snd;
+//				m_right_child(snd) = fst_right_child;
+//				m_left_child(snd) = fst;
+//			} else {
+//				m_parent(fst_left_child) = snd;
+//				m_left_child(snd) = fst_left_child;
+//				m_right_child(snd) = fst;
+//			}
+//		}
+	}
+
+
+/*******************************************************************************
+ * Heap Operations
+ ******************************************************************************/
+
+	/**
+	 * Checks whether the heap is empty.
+	 */
+public:
+	bool empty() {
+		return m_root == 0;
+	}
+
+	/**
+	 * Returns a reference to the top of the heap.
+	 */
+public:
+	value_type& top() {
+		return m_root;
+	}
+
+	/**
+	 * Pops the top of the heap.
+	 */
+public:
+	void pop() {
+		assert(top());
+		assert(m_last);
+
+//		dump(m_root);
+//		dump(m_last);
+
+		value_type new_last = m_last;
+		if(m_root != m_last) {
+
+//			std::cout << "gonna	 exchange" << std::endl;
+			exchange(m_root, m_last);
+
+			assert(m_root == new_last);
+
+			new_last = m_last;
+//			std::cout << "donna	 exchange" << std::endl;
+
+//			dump(m_root);
+//			dump(m_last);
+
+			// Determine the new "last" node of the tree.
+			// The algorithm also works when the last node at the deepest level
+			// is removed. The algorithm will then continue from the rightmost
+			// node at the level one lower.
+			while(
+					m_parent(new_last) &&
+					m_left_child(m_parent(new_last)) == new_last
+			) {
+				new_last = m_parent(new_last);
+//				std::cout << "moving up ..." << std::endl;
+			}
+//			std::cout << "stopped moving up ..." << std::endl;
+			if(m_parent(new_last)) {
+				new_last = m_left_child(m_parent(new_last));
+//				std::cout << "moved left ..." << std::endl;
+			}
+			while(m_right_child(new_last)) {
+				new_last = m_right_child(new_last);
+//				std::cout << "moving right ..." << std::endl;
+			}
+//			std::cout << "stopped moving right ..." << std::endl;
+		}
+
+		if( m_parent(m_last) ) {
+			if( m_left_child(m_parent(m_last)) == m_last ) {
+				m_left_child(m_parent(m_last)) = 0;
+			} else {
+				m_right_child(m_parent(m_last)) = 0;
+			}
+		}
+		m_parent(m_last) = 0;
+		m_left_child(m_last) = 0;
+		m_right_child(m_last) = 0;
+//		m_key(m_last) = -1;
+
+		if(m_root != m_last) {
+			m_last = new_last;
+			min_heapify(m_root);
+		} else {
+			m_root = m_last = 0;
+		}
+	}
+
+	/**
+	 * Updates the key of a given node.
+	 */
+public:
+	void update_key(value_type node, key_type new_key) {
+		m_key(node) = new_key;
+
+		value_type seek = node;
+		while(
+				m_parent(seek) &&
+				m_compare_keys(seek, m_parent(seek))
+		) {
+			exchange(seek, m_parent(seek));
+		}
+	}
+
+	/**
+	 * Inserts a new value into the binary heap with the given key.
+	 */
+public:
+	void insert(value_type node, key_type key) {
+
+		// Special case handling an empty heap.
+		if(m_last == m_root && m_last == 0) {
+			m_last = m_root = node;
+			m_parent(node) = 0;
+			m_left_child(node) 	= 0;
+			m_right_child(node) = 0;
+			m_key(node) = key;
+			return;
+		}
+
+		// Search the new "last" position to insert the new node at.
+		value_type new_last = m_last;
+
+		// While the current node is right w.r.t. its parent, move up the tree.
+		while(
+				m_parent(new_last) &&
+				m_right_child(m_parent(new_last)) == new_last
+		) {
+			new_last = m_parent(new_last);
+		}
+		// If there is a right sibling, go to it and descent to the left.
+		if(m_parent(new_last) && m_right_child(m_parent(new_last))) {
+			new_last = m_right_child(m_parent(new_last));
+			while(m_left_child(new_last)) {
+				new_last = m_left_child(new_last);
+			}
+		}
+		// If there is a parent at this point, the right sibling does not exist,
+		// and the given node should be placed there.
+		if(m_parent(new_last)) {
+			// Add the node to the tree.
+			if(m_right_child(m_parent(new_last)) == 0) {
+				new_last = m_parent(new_last);
+				m_right_child(new_last) = node;
+			} else {
+				m_left_child(new_last) = node;
+			}
+		} else {
+			// We came from the very last node at the current level and hence
+			// must begin a new level. Descent to the left.
+			while( m_left_child(new_last) ) {
+				new_last = m_left_child(new_last);
+			}
+
+			// Add the node to the tree.
+			m_left_child(new_last) = node;
+		}
+
+		m_parent(node) = new_last;
+		m_left_child(node) 	= 0;
+		m_right_child(node) = 0;
+		m_key(node) = key;
+		m_last = node;
+
+		// Update the key and force the heap to rebalance itself.
+		update_key(node, key);
+	}
+
+	/**
+	 * Removes the node. The operation can safely be applied to an element no
+	 * longer in the heap.
+	 */
+public:
+	void remove(value_type node) {
+		if( m_parent(node) == 0 && m_root != node) {
+			return;
+		}
+
+		// Move the node up the tree.
+		value_type seek = node;
+		while(m_parent(seek)) {
+			exchange(seek, m_parent(seek));
+		}
+		// Remove it.
+		pop();
+	}
+
+
+public:
+	void dump(value_type node) {
+		std::cout << "seq: " << node->get_sequence_number() << std::endl;
+		std::cout << "key: " << m_key(node) << std::endl;
+		std::cout << "par: " << (m_parent(node) ? (m_parent(node))->get_sequence_number() : -1) << std::endl;
+		std::cout << "lef: " << (m_left_child(node) ? (m_left_child(node))->get_sequence_number() : -1) << std::endl;
+		std::cout << "rig: " << (m_right_child(node) ? (m_right_child(node))->get_sequence_number() : -1) << std::endl;
+
+	}
+
+/*******************************************************************************
+ * Data Members
+ ******************************************************************************/
+
+
+public:
+	/**
+	 * The root of the binary heap.
+	 */
+	ValueType m_root;
+
+	/**
+	 * The "last" node of the binary heap. The last node is the node that can
+	 * safely be removed such that the resulting tree is still a complete binary
+	 * tree, with the additional constraint that the nodes in the last level are
+	 * all in the "left" part of the tree.
+	 */
+	ValueType m_last;
+
+	/**
+	 * A functor to compare keys.
+	 */
+	KeyComparator m_compare_keys;
+
+	/**
+	 * A functor to obtain the key.
+	 */
+	GetKey m_key;
+
+	/**
+	 * A functor to obtain the parent.
+	 */
+	GetParent m_parent;
+
+	/**
+	 * A functor to obtain the left child.
+	 */
+	GetLeft m_left_child;
+
+	/**
+	 * A functor to obtain the right child.
+	 */
+	GetRight m_right_child;
+
+};
+
+}
+
+#endif // MTL_BINARY_HEAP_INCLUDE

AMDiS/lib/mtl4/boost/numeric/itl/pc/comparators.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+//
+// Algorithm inspired by Nick Vannieuwenhoven, written by Cornelius Steinhardt
+
+
+
+/*
+ *
+ *  Created on: Oct 14, 2010
+ *      Author: heazk
+ */
+
+
+#ifndef MTL_COMPARATORS_INCLUDE
+#define MTL_COMPARATORS_INCLUDE
+
+#define POINTER long
+
+namespace compare {
+
+/**
+ * Compares two given pointers based on the address to which they point.
+ */
+template<class Type>
+struct address_compare {
+	bool operator()(const Type *const a, const Type *const b) const {
+		return reinterpret_cast<POINTER>(a) < reinterpret_cast<POINTER>(b);
+	}
+};
+
+/**
+ * Compares two given pointers based on the address to which they point.
+ */
+template<class Type>
+struct address_compare_equal {
+	bool operator()(const Type *const a, const Type *const b) const {
+		return reinterpret_cast<POINTER>(a) == reinterpret_cast<POINTER>(b);
+	}
+};
+
+/**
+ * Hashes a given pointer based on its address.
+ */
+template<class Type>
+struct address_hasher {
+	POINTER operator()(const Type *const a) const {
+		return reinterpret_cast<POINTER>(a);
+	}
+};
+
+} // end namespace compare
+
+
+#endif // MTL_COMPARATORS_INCLUDE

AMDiS/lib/mtl4/boost/numeric/itl/pc/concat.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG, www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also tools/license/license.mtl.txt in the distribution.
+
+#ifndef ITL_PC_CONCAT_INCLUDE
+#define ITL_PC_CONCAT_INCLUDE
+
+#include <boost/mpl/if.hpp>
+#include <boost/numeric/mtl/operation/resource.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+#include <boost/numeric/itl/pc/solver.hpp>
+
+namespace itl { namespace pc {
+
+/// Class for concatenating \tparam PC1 and \tparam PC2
+template <typename PC1, typename PC2, typename Matrix, bool Store1= true, bool Store2= true>
+class concat
+{
+    typedef typename boost::mpl::if_c<Store1, PC1, const PC1&>::type pc1_type;
+    typedef typename boost::mpl::if_c<Store2, PC2, const PC2&>::type pc2_type;
+
+  public:
+    /// Construct both preconditioners from matrix \p A
+    explicit concat(const Matrix& A) : A(A), pc1(A), pc2(A)
+    {
+	BOOST_STATIC_ASSERT((Store1 && Store2));
+    }
+
+    /// Both preconditioners are already constructed and passed as arguments
+    /** If pc1 or pc2 is only constructed temporarily in the constructor call,
+	the according Store argument must be true; otherwise the preconditioner
+	will be a stale reference.
+	Conversely, if the preconditioner is already build outside the constructor call,
+	the according Store argument should be false for not storing the preconditioner twice. **/
+    concat(const Matrix& A, const PC1& pc1, const PC2& pc2) : A(A), pc1(pc1), pc2(pc2) {}
+
+  private:
+    template <typename VectorOut>
+    VectorOut& create_r(const VectorOut& y) const
+    {
+	static VectorOut  r(resource(y));
+	return r;
+    }
+
+    template <typename VectorOut>
+    VectorOut& create_d(const VectorOut& y) const
+    {
+	static VectorOut  d(resource(y));
+	return d;
+    }
+
+  public:
+    /// Concatenated preconditioning: pc2 is applied regularly and pc1 afterwards by defect correction
+    template <typename VectorIn, typename VectorOut>
+    void solve(const VectorIn& x, VectorOut& y) const
+    {
+	mtl::vampir_trace<5058> tracer;
+	y.checked_change_resource(x);
+	pc2.solve(x, y);
+
+	VectorOut &r= create_r(y), &d= create_d(y);
+	r= x;
+	r-= A * y;
+
+	pc1.solve(r, d);
+	y+= d;
+    }
+
+
+    /// Concatenated preconditioning: adjoint of pc1 is applied regularly and pc2's adjoint afterwards by defect correction
+    template <typename VectorIn, typename VectorOut>
+    void adjoint_solve(const VectorIn& x, VectorOut& y) const
+    {
+	mtl::vampir_trace<5059> tracer;
+	y.checked_change_resource(x);
+	pc1.adjoint_solve(x, y);
+
+	VectorOut &r= create_r(y), &d= create_d(y);
+	r= x;
+	r-= adjoint(A) * y;
+
+	pc2.adjoint_solve(r, d);
+	y+= d;
+    }
+
+   private:
+    const Matrix& A;
+    pc1_type 	  pc1;
+    pc2_type   	  pc2;
+};
+
+template <typename PC1, typename PC2, typename Matrix, bool Store1, bool Store2, typename Vector>
+solver<concat<PC1, PC2, Matrix, Store1, Store2>, Vector, false>
+inline solve(const concat<PC1, PC2, Matrix, Store1, Store2>& P, const Vector& x)
+{
+    return solver<concat<PC1, PC2, Matrix, Store1, Store2>, Vector, false>(P, x);
+}
+
+template <typename PC1, typename PC2, typename Matrix, bool Store1, bool Store2, typename Vector>
+solver<concat<PC1, PC2, Matrix, Store1, Store2>, Vector, true>
+inline adjoint_solve(const concat<PC1, PC2, Matrix, Store1, Store2>& P, const Vector& x)
+{
+    return solver<concat<PC1, PC2, Matrix, Store1, Store2>, Vector, true>(P, x);
+}
+
+
+}} // namespace itl::pc
+
+#endif // ITL_PC_CONCAT_INCLUDE

AMDiS/lib/mtl4/boost/numeric/itl/pc/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 ITL_PC_DIAGONAL_INCLUDE
+#define ITL_PC_DIAGONAL_INCLUDE
+
+#include <boost/numeric/linear_algebra/inverse.hpp>
+
+#include <boost/numeric/mtl/vector/dense_vector.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/operation/resource.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+#include <boost/numeric/itl/pc/solver.hpp>
+
+namespace itl { namespace pc {
+
+/// Diagonal Preconditioner
+template <typename Matrix, typename Value= typename mtl::Collection<Matrix>::value_type>
+class diagonal
+{
+  public:
+    typedef Value                                         value_type;
+    typedef typename mtl::Collection<Matrix>::size_type   size_type;
+    typedef diagonal                                      self;
+
+    /// Constructor takes matrix reference
+    explicit diagonal(const Matrix& A) : inv_diag(num_rows(A))
+    {
+	mtl::vampir_trace<5050> tracer;
+	MTL_THROW_IF(num_rows(A) != num_cols(A), mtl::matrix_not_square());
+	using math::reciprocal;
+
+	for (size_type i= 0; i < num_rows(A); ++i)
+	    inv_diag[i]= reciprocal(A[i][i]);
+    }
+
+    /// Member function solve, better use free function solve
+    template <typename Vector>
+    Vector solve(const Vector& x) const
+    {
+	Vector y(resource(x));
+	solve(x, y);
+	return y;
+    }
+
+    template <typename VectorIn, typename VectorOut>
+    void solve(const VectorIn& x, VectorOut& y) const
+    {
+	mtl::vampir_trace<5051> tracer;
+	y.checked_change_resource(x);
+	MTL_THROW_IF(size(x) != size(inv_diag), mtl::incompatible_size());
+	for (size_type i= 0; i < size(inv_diag); ++i)
+	    y[i]= inv_diag[i] * x[i];
+    }
+
+    /// Member function for solving adjoint problem, better use free function adjoint_solve
+    template <typename Vector>
+    Vector adjoint_solve(const Vector& x) const
+    {
+	Vector y(resource(x));
+	adjoint_solve(x, y);
+	return y;
+    }
+
+    template <typename VectorIn, typename VectorOut>
+    void adjoint_solve(const VectorIn& x, VectorOut& y) const
+    {
+	using mtl::conj;
+	y.checked_change_resource(x);
+	MTL_THROW_IF(size(x) != size(inv_diag), mtl::incompatible_size());
+	for (size_type i= 0; i < size(inv_diag); ++i)
+	    y[i]= conj(inv_diag[i]) * x[i];
+    }
+
+ protected:
+    mtl::vector::dense_vector<value_type>    inv_diag;
+};
+
+/// Solve approximately a sparse system in terms of inverse diagonal
+template <typename Matrix, typename Vector>
+solver<diagonal<Matrix>, Vector, false>
+inline solve(const diagonal<Matrix>& P, const Vector& x)
+{
+    return solver<diagonal<Matrix>, Vector, false>(P, x);
+}
+
+/// Solve approximately the adjoint of a sparse system in terms of inverse diagonal
+template <typename Matrix, typename Vector>
+solver<diagonal<Matrix>, Vector, true>
+inline adjoint_solve(const diagonal<Matrix>& P, const Vector& x)
+{
+    return solver<diagonal<Matrix>, Vector, true>(P, x);
+}
+
+
+}} // namespace itl::pc
+
+#endif // ITL_PC_DIAGONAL_INCLUDE

AMDiS/lib/mtl4/boost/numeric/itl/pc/ic_0.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef ITL_PC_IC_0_INCLUDE
+#define ITL_PC_IC_0_INCLUDE
+
+#include <boost/mpl/bool.hpp>
+
+#include <boost/numeric/linear_algebra/identity.hpp>
+#include <boost/numeric/linear_algebra/inverse.hpp>
+
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/utility/ashape.hpp>
+#include <boost/numeric/mtl/utility/tag.hpp>
+#include <boost/numeric/mtl/utility/category.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/mtl/operation/lower_trisolve.hpp>
+#include <boost/numeric/mtl/operation/upper_trisolve.hpp>
+#include <boost/numeric/mtl/matrix/upper.hpp>
+#include <boost/numeric/mtl/matrix/strict_lower.hpp>
+#include <boost/numeric/mtl/matrix/compressed2D.hpp>
+#include <boost/numeric/mtl/matrix/parameter.hpp>
+#include <boost/numeric/mtl/matrix/transposed_view.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+#include <boost/numeric/itl/pc/solver.hpp>
+
+
+namespace itl { namespace pc {
+
+template <typename Matrix, typename Value= typename mtl::Collection<Matrix>::value_type>
+class ic_0
+{
+  public:
+    typedef Value                                                    value_type;
+    typedef typename mtl::Collection<Matrix>::size_type              size_type;
+    typedef ic_0                                                     self;
+
+    typedef mtl::matrix::parameters<mtl::row_major, mtl::index::c_index, mtl::non_fixed::dimensions, false, size_type> para;
+    typedef mtl::matrix::compressed2D<value_type, para>                      U_type;
+#ifndef ITL_IC_0_ONE_MATRIX
+    typedef U_type                                                   L_type;
+#else
+    typedef typename mtl::matrix::transposed_view<U_type>            L_type;
+#endif
+    typedef mtl::matrix::detail::lower_trisolve_t<L_type, mtl::tag::inverse_diagonal, true> lower_solver_t;
+    typedef mtl::matrix::detail::upper_trisolve_t<U_type, mtl::tag::inverse_diagonal, true> upper_solver_t;
+
+    ic_0(const Matrix& A) : f(A, U), L(trans(U)), lower_solver(L), upper_solver(U) {}
+
+
+    // solve x = U^* U y --> y= U^{-1} U^{-*} x
+    template <typename Vector>
+    Vector solve(const Vector& x) const
+    {
+	mtl::vampir_trace<5036> tracer;
+	return inverse_upper_trisolve(U, inverse_lower_trisolve(adjoint(U), x));
+    }
+
+    // solve x = U^* y --> y0= U^{-*} x
+    template <typename VectorIn, typename VectorOut>
+    const VectorOut& solve_lower(const VectorIn& x, VectorOut&) const
+    {
+	static VectorOut y0;
+	y0.change_resource(resource(x));
+	lower_solver(x, y0);
+	return y0;
+    }
+
+    // solve x = U^* U y --> y= U^{-1} U^{-*} x
+    template <typename VectorIn, typename VectorOut>
+    void solve(const VectorIn& x, VectorOut& y) const
+    {
+	mtl::vampir_trace<5037> tracer;
+	const VectorOut& y0= solve_lower(x, y);
+
+	y.checked_change_resource(x);
+	upper_solver(y0, y);
+    }
+
+    // solve x = (LU)^* y --> y= L^{-*} U^{-*} x
+    template <typename Vector>
+    Vector adjoint_solve(const Vector& x) const
+    {
+	mtl::vampir_trace<5044> tracer;
+	return solve(x);
+    }
+
+    // solve x = (LU)^* y --> y= L^{-*} U^{-*} x
+    template <typename VectorIn, typename VectorOut>
+    void adjoint_solve(const VectorIn& x, VectorOut& y) const
+    {
+	mtl::vampir_trace<5044> tracer;
+	solve(x, y);
+    }
+
+
+    L_type get_L() { return L_type(L); }
+    U_type get_U() { return U; }
+
+  protected:
+    template <typename VectorOut, typename Solver> friend struct ic_0_evaluator;
+
+    // Dummy type to perform factorization in initializer list not in
+    struct factorizer
+    {
+	factorizer(const Matrix &A, U_type& U)
+	{   factorize(A, U, mtl::traits::is_sparse<Matrix>(), boost::is_same<Value, typename mtl::Collection<Matrix>::value_type>());  }
+
+	template <typename T>
+	void factorize(const Matrix&, U_type&, boost::mpl::false_, T)
+	{   MTL_THROW_IF(true, mtl::logic_error("IC(0) is not suited for dense matrices"));	}
+
+	// When we change the value_type then the factorization is still performed with that of A
+	template <typename UF>
+	void factorize(const Matrix& A, UF& U, boost::mpl::true_, boost::mpl::false_)
+	{
+	    typedef mtl::matrix::compressed2D<typename mtl::Collection<Matrix>::value_type, para> tmp_type;
+	    tmp_type U_tmp;
+	    factorize(A, U_tmp, boost::mpl::true_(), boost::mpl::true_());
+	    U= U_tmp;
+	}
+
+	// Factorization adapted from Saad
+	// Undefined (runtime) behavior if matrix is not symmetric
+	// UF is type for the factorization
+	template <typename UF>
+	void factorize(const Matrix& A, UF& U, boost::mpl::true_, boost::mpl::true_)
+	{
+	    using namespace mtl; using namespace mtl::tag;  using mtl::traits::range_generator;
+	    using math::reciprocal; using mtl::matrix::upper;
+	    mtl::vampir_trace<5035> tracer;
+
+	    // For the factorization we take still the value_type of A and later we copy it maybe to another value_type
+	    typedef typename mtl::Collection<Matrix>::value_type      value_type;
+	    typedef typename range_generator<row, UF>::type       cur_type;
+	    typedef typename range_generator<nz, cur_type>::type      icur_type;
+
+	    MTL_THROW_IF(num_rows(A) != num_cols(A), mtl::matrix_not_square());
+	    U= upper(A);
+
+	    typename mtl::traits::col<UF>::type                   col(U);
+	    typename mtl::traits::value<UF>::type                 value(U);
+
+	    cur_type kc= begin<row>(U), kend= end<row>(U);
+	    for (size_type k= 0; kc != kend; ++kc, ++k) {
+
+		icur_type ic= begin<nz>(kc), iend= end<nz>(kc);
+		MTL_DEBUG_THROW_IF(col(*ic) != k, mtl::missing_diagonal());
+
+		// U[k][k]= 1.0 / sqrt(U[k][k]);
+		value_type inv_dia= reciprocal(sqrt(value(*ic)));
+		value(*ic, inv_dia);
+		// icur_type jbegin=
+		++ic;
+		for (; ic != iend; ++ic) {
+		    // U[k][i] *= U[k][k]
+		    value_type d= value(*ic) * inv_dia;
+		    value(*ic, d);
+		    size_type i= col(*ic);
+
+		    // find non-zeros U[j][i] below U[k][i] for j in (k, i]
+		    // 1. Go to ith row in U (== ith column in U)
+		    cur_type irow(i, U); // = begin<row>(U); irow+= i;
+		    // 2. Find nonzeros with col() in (k, i]
+		    icur_type jc= begin<nz>(irow), jend= end<nz>(irow);
+		    while (col(*jc) <= k)  ++jc;
+		    while (col(*--jend) > i) ;
+		    ++jend;
+
+		    for (; jc != jend; ++jc) {
+			size_type j= col(*jc);
+			U.lvalue(j, i)-= d * U[k][j];
+		    }
+		    // std::cout << "U after eliminating U[" << i << "][" << k << "] =\n" << U;
+		}
+	    }
+	}
+    };
+
+    U_type                       U;
+    factorizer                   f;
+    L_type                       L;
+    lower_solver_t               lower_solver;
+    upper_solver_t               upper_solver;
+};
+
+#if 0
+template <typename Matrix, typename Value, typename Vector>
+struct ic_0_solver
+  : mtl::vector::assigner<ic_0_solver<Matrix, Value, Vector> >
+{
+    typedef ic_0<Matrix, Value> pc_type;
+
+    ic_0_solver(const ic_0<Matrix, Value>& P, const Vector& x) : P(P), x(x) {}
+
+    template <typename VectorOut>
+    void assign_to(VectorOut& y) const
+    {	P.solve(x, y);    }
+
+    const ic_0<Matrix, Value>& P;
+    const Vector&              x;
+};
+#endif
+
+template <typename VectorOut, typename Solver>
+struct ic_0_evaluator
+{
+    typedef typename Solver::pc_type                        pc_type;
+    typedef typename pc_type::size_type                     size_type;
+    typedef typename mtl::Collection<VectorOut>::value_type out_value_type;
+
+
+    ic_0_evaluator(VectorOut& y, const Solver& s)
+      : y(y), s(s), U(s.P.U), y0(s.P.solve_lower(s.x, y)) { MTL_DEBUG_ARG(lr= 99999999); }
+
+
+    void operator()(size_type i) { at<0>(i); }
+    void operator[](size_type i) { at<0>(i); }
+
+    template <unsigned Offset>
+    void at(size_type r)
+    {
+#ifndef NDEBUG
+	MTL_THROW_IF(r+Offset >= lr, mtl::logic_error("Traversal must be backward")); lr= r+Offset;
+#endif
+	size_type j0= U.ref_major()[r+Offset];
+	const size_type cj1= U.ref_major()[r+Offset+1];
+	MTL_DEBUG_THROW_IF(j0 == cj1 || U.ref_minor()[j0] != r+Offset, mtl::missing_diagonal());
+	out_value_type rr= y0[r+Offset], dia= U.data[j0++];
+	for (; j0 != cj1; ++j0) {
+	    MTL_DEBUG_THROW_IF(U.ref_minor()[j0] <= r+Offset, mtl::logic_error("Matrix entries must be sorted for this."));
+	    rr-= U.data[j0] * y[U.ref_minor()[j0]];
+	}
+	y[r+Offset]= rr * dia;
+    }
+
+    VectorOut&                               y;
+    const Solver&                            s;
+    const typename pc_type::U_type&          U;
+    const VectorOut&                         y0;
+    MTL_DEBUG_ARG(size_type                  lr;)
+};
+
+template <typename VectorOut, typename Solver>
+inline std::size_t size(const ic_0_evaluator<VectorOut, Solver>& eval)
+{   return size(eval.y); }
+
+template <typename Matrix, typename Value, typename Vector>
+solver<ic_0<Matrix, Value>, Vector, false>
+inline solve(const ic_0<Matrix, Value>& P, const Vector& x)
+{
+    return solver<ic_0<Matrix, Value>, Vector, false>(P, x);
+}
+
+template <typename Matrix, typename Value, typename Vector>
+solver<ic_0<Matrix, Value>, Vector, true>
+inline adjoint_solve(const ic_0<Matrix, Value>& P, const Vector& x)
+{
+    return solver<ic_0<Matrix, Value>, Vector, true>(P, x);
+}
+
+
+}} // namespace itl::pc
+
+namespace mtl { namespace vector {
+    using itl::pc::size;
+}} // namespace mtl::vector
+
+#endif // ITL_PC_IC_0_INCLUDE

AMDiS/lib/mtl4/boost/numeric/itl/pc/identity.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef ITL_PC_IDENTITY_INCLUDE
+#define ITL_PC_IDENTITY_INCLUDE
+
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+#include <boost/numeric/itl/pc/solver.hpp>
+
+namespace itl { namespace pc {
+
+/// Identity preconditioner, i.e. no preconditioning and vector is just copied
+/** Second template is just for a uniform interface with other preconditioners. **/
+template <typename Matrix, typename Value= double>
+class identity
+{
+  public:
+    typedef typename mtl::Collection<Matrix>::value_type  value_type;
+    typedef typename mtl::Collection<Matrix>::size_type   size_type;
+    typedef identity                                      self;
+
+    identity(const Matrix&) {}
+
+    template <typename Vector>
+    Vector solve(const Vector& x) const
+    {
+	mtl::vampir_trace<5032> tracer;
+	return x;
+    }
+
+    template <typename VectorIn, typename VectorOut>
+    void solve(const VectorIn& b, VectorOut& x) const
+    {
+	mtl::vampir_trace<5032> tracer;
+	x= b;
+    }
+
+    template <typename Vector>
+    Vector adjoint_solve(const Vector& x) const
+    {
+	mtl::vampir_trace<5034> tracer;
+	return x;
+    }
+
+    template <typename VectorIn, typename VectorOut>
+    void adjoint_solve(const VectorIn& b, VectorOut& x) const
+    {
+	mtl::vampir_trace<5034> tracer;
+	x= b;
+    }
+};
+
+template <typename Matrix, typename Vector>
+solver<identity<Matrix>, Vector, false>
+inline solve(const identity<Matrix>& P, const Vector& x)
+{   return solver<identity<Matrix>, Vector, false>(P, x); }
+
+template <typename Matrix, typename Vector>
+solver<identity<Matrix>, Vector, true>
+inline adjoint_solve(const identity<Matrix>& P, const Vector& x)
+{   return solver<identity<Matrix>, Vector, true>(P, x); }
+
+
+}} // namespace itl::pc
+
+#endif // ITL_PC_IDENTITY_INCLUDE

AMDiS/lib/mtl4/boost/numeric/itl/pc/ilu.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef ITL_PC_ILU_INCLUDE
+#define ITL_PC_ILU_INCLUDE
+
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/utility/tag.hpp>
+#include <boost/numeric/mtl/operation/adjoint.hpp>
+#include <boost/numeric/mtl/operation/lower_trisolve.hpp>
+#include <boost/numeric/mtl/operation/upper_trisolve.hpp>
+#include <boost/numeric/mtl/operation/lu.hpp>
+#include <boost/numeric/mtl/matrix/parameter.hpp>
+#include <boost/numeric/mtl/matrix/compressed2D.hpp>
+#include <boost/numeric/mtl/matrix/dense2D.hpp>
+#include <boost/numeric/mtl/vector/dense_vector.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+#include <boost/numeric/itl/pc/solver.hpp>
+
+namespace itl { namespace pc {
+
+/// Incomplete LU factorization of a \p Matrix into matrices of type \p Value using a \p Factorizer (e.g. ILU(0) or ILUT)
+template <typename Matrix, typename Factorizer, typename Value= typename mtl::Collection<Matrix>::value_type>
+class ilu
+{
+  public:
+    typedef Value                                         value_type;
+    typedef typename mtl::Collection<Matrix>::size_type   size_type;
+    typedef ilu                                           self;
+    typedef Factorizer                                    factorizer_type;
+
+    typedef mtl::matrix::parameters<mtl::row_major, mtl::index::c_index, mtl::non_fixed::dimensions, false, size_type> para;
+    typedef mtl::matrix::compressed2D<value_type, para>                     L_type;
+    typedef mtl::matrix::compressed2D<value_type, para>                     U_type;
+    typedef typename mtl::matrix::traits::adjoint<L_type>::type     adjoint_L_type;
+    typedef typename mtl::matrix::traits::adjoint<U_type>::type     adjoint_U_type;
+
+    typedef mtl::matrix::detail::lower_trisolve_t<L_type, mtl::tag::unit_diagonal, true>    lower_solver_t;
+    typedef mtl::matrix::detail::upper_trisolve_t<U_type, mtl::tag::inverse_diagonal, true> upper_solver_t;
+
+    typedef mtl::matrix::detail::lower_trisolve_t<adjoint_U_type, mtl::tag::inverse_diagonal, true> adjoint_lower_solver_t;
+    typedef mtl::matrix::detail::upper_trisolve_t<adjoint_L_type, mtl::tag::unit_diagonal, true>    adjoint_upper_solver_t;
+
+    /// Factorization adapted from Saad
+    explicit ilu(const Matrix& A)
+      : f(A, L, U), lower_solver(L), upper_solver(U), adjoint_L(adjoint(L)), adjoint_U(adjoint(U)),
+	adjoint_lower_solver(adjoint_U), adjoint_upper_solver(adjoint_L) {}
+
+    template <typename FactPara>
+    ilu(const Matrix& A, const FactPara& p)
+      : f(A, p, L, U), lower_solver(L), upper_solver(U), adjoint_L(adjoint(L)), adjoint_U(adjoint(U)),
+	adjoint_lower_solver(adjoint_U), adjoint_upper_solver(adjoint_L) {}
+
+    /// Solve  LU y = x --> y= U^{-1} L^{-1} x
+    template <typename Vector>
+    Vector solve(const Vector& x) const
+    {
+	Vector y;
+	solve(x, y);
+	return y;
+    }
+
+    // solve x = L y --> y0= L^{-1} x
+    template <typename VectorIn, typename VectorOut>
+    const VectorOut& solve_lower(const VectorIn& x, VectorOut&) const
+    {
+	static VectorOut y0;
+	y0.change_resource(resource(x));
+	lower_solver(x, y0);
+	return y0;
+    }
+
+    // Solve  LU y = x --> y= U^{-1} L^{-1} x
+    template <typename VectorIn, typename VectorOut>
+    void solve(const VectorIn& x, VectorOut& y) const
+    {
+	mtl::vampir_trace<5039> tracer;
+	const VectorOut& y0= solve_lower(x, y);
+
+	y.checked_change_resource(x);
+	upper_solver(y0, y);
+    }
+
+
+    /// Solve (LU)^H y = x --> y= L^{-H} U^{-H} x
+    template <typename Vector>
+    Vector adjoint_solve(const Vector& x) const
+    {
+	mtl::vampir_trace<5040> tracer;
+	Vector y(resource(x));
+	adjoint_solve(x, y);
+	return y;
+    }
+
+    /// Solve (LU)^H y = x --> y= L^{-H} U^{-H} x
+    template <typename VectorIn, typename VectorOut>
+    void adjoint_solve(const VectorIn& x, VectorOut& y) const
+    {
+	mtl::vampir_trace<5040> tracer;
+	y.checked_change_resource(x);
+	// y= unit_upper_trisolve(adjoint(L), inverse_lower_trisolve(adjoint(U), x));
+	static VectorOut y0;
+	y0.change_resource(resource(x));
+	adjoint_lower_solver(x, y0);
+	adjoint_upper_solver(y0, y);
+    }
+
+
+    L_type get_L() { return L; }
+    U_type get_U() { return U; }
+
+  public:
+    L_type                      L;
+    U_type                      U;
+  private:
+    Factorizer                  f;
+    lower_solver_t              lower_solver;
+    upper_solver_t              upper_solver;
+    adjoint_L_type              adjoint_L;
+    adjoint_U_type              adjoint_U;
+    adjoint_lower_solver_t      adjoint_lower_solver;
+    adjoint_upper_solver_t      adjoint_upper_solver;
+
+};
+
+template <typename Value, typename Factorizer, typename V2>
+class ilu<mtl::matrix::dense2D<Value, mtl::matrix::parameters<> >, Factorizer, V2> // last 2 arguments are dummies
+{
+  public:
+    typedef mtl::matrix::dense2D<Value, mtl::matrix::parameters<> >    Matrix;
+    typedef typename mtl::Collection<Matrix>::value_type  value_type;
+    typedef typename mtl::Collection<Matrix>::size_type   size_type;
+    typedef ilu                                           self;
+    typedef Matrix                                        LU_type;
+
+    ilu(const Matrix& A) : LU(A) { lu(LU, P); std::cout << "LU is\n" << LU << "P is " << P << "\n"; }
+
+    // Solve  P^{-1}LU x = b --> x= U^{-1} L^{-1} P b
+    template <typename Vector>
+    Vector solve(const Vector& b) const { return lu_apply(LU, P, b); }
+
+    // Solve  P^{-1}LU x = b --> x= U^{-1} L^{-1} P b
+    template <typename VectorIn, typename VectorOut>
+    void solve(const VectorIn& b, VectorOut& x) const { x= lu_apply(LU, P, b); }
+
+    // Solve (P^{-1}LU)^H x = b --> x= P^{-1}L^{-H} U^{-H} b // P^{-1}^{-1}^H = P^{-1})
+    template <typename Vector>
+    Vector adjoint_solve(const Vector& b) const { return lu_adjoint_apply(LU, P, b); }
+
+    // Solve (P^{-1}LU)^H x = b --> x= P^{-1}L^{-H} U^{-H} b // P^{-1}^{-1}^H = P^{-1})
+    template <typename VectorIn, typename VectorOut>
+    void adjoint_solve(const VectorIn& b, VectorOut& x) const { x= lu_adjoint_apply(LU, P, b); }
+
+  private:
+    LU_type                        LU;
+    mtl::vector::dense_vector<size_type, mtl::vector::parameters<> >   P;
+};
+
+#if 0
+template <typename Matrix, typename Factorizer, typename Value, typename Vector>
+struct ilu_solver
+  : mtl::vector::assigner<ilu_solver<Matrix, Factorizer, Value, Vector> >
+{
+    typedef ilu<Matrix, Factorizer, Value> pc_type;
+
+    ilu_solver(const pc_type& P, const Vector& x) : P(P), x(x) {}
+
+    template <typename VectorOut>
+    void assign_to(VectorOut& y) const
+    {	P.solve(x, y);    }
+
+    const pc_type&       P;
+    const Vector&        x;
+};
+#endif
+
+
+
+/// Solve LU x = b --> x= U^{-1} L^{-1} b
+template <typename Matrix, typename Factorizer, typename Value, typename Vector>
+// ilu_solver<Matrix, Factorizer, Value, Vector>
+solver<ilu<Matrix, Factorizer, Value>, Vector, false>
+inline solve(const ilu<Matrix, Factorizer, Value>& P, const Vector& x)
+{
+    return solver<ilu<Matrix, Factorizer, Value>, Vector, false>(P, x);
+}
+
+
+/// Solve (LU)^H x = b --> x= L^{-H} U^{-H} b
+template <typename Matrix, typename Factorizer, typename Value, typename Vector>
+// Vector
+solver<ilu<Matrix, Factorizer, Value>, Vector, true>
+inline adjoint_solve(const ilu<Matrix, Factorizer, Value>& P, const Vector& b)
+{
+    return solver<ilu<Matrix, Factorizer, Value>, Vector, true>(P, b);
+}
+
+// ic_0_evaluator not needed IC(0) and ILU(0) do the same at the upper triangle ;-)
+
+}} // namespace itl::pc
+
+
+#endif // ITL_PC_ILU_INCLUDE

AMDiS/lib/mtl4/boost/numeric/itl/pc/ilu_0.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef ITL_PC_ILU_0_INCLUDE
+#define ITL_PC_ILU_0_INCLUDE
+
+#include <boost/mpl/bool.hpp>
+#include <boost/numeric/linear_algebra/identity.hpp>
+
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/operation/invert_diagonal.hpp>
+#include <boost/numeric/mtl/utility/tag.hpp>
+#include <boost/numeric/mtl/utility/category.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/mtl/matrix/strict_lower.hpp>
+#include <boost/numeric/mtl/matrix/upper.hpp>
+#include <boost/numeric/mtl/matrix/compressed2D.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+#include <boost/numeric/itl/pc/ilu.hpp>
+
+namespace itl { namespace pc {
+
+// Dummy type to perform factorization in initializer list not in
+struct ilu_0_factorizer
+{
+    template <typename Matrix, typename L_type, typename U_type>
+    ilu_0_factorizer(const Matrix &A, L_type& L, U_type& U)
+    {   factorize(A, L, U, mtl::traits::is_sparse<Matrix>());  }
+
+    template <typename Matrix, typename L_type, typename U_type>
+    void factorize(const Matrix&, L_type&, U_type&, boost::mpl::false_)
+    {  MTL_THROW_IF(true, mtl::logic_error("ILU is not intended for dense matrices")); }
+
+    template <typename Matrix, typename L_type, typename U_type>
+    void factorize(const Matrix& A, L_type& L, U_type& U, boost::mpl::true_)
+    {
+	using namespace mtl; using namespace mtl::tag;  using mtl::traits::range_generator;
+	using math::reciprocal;
+	MTL_THROW_IF(num_rows(A) != num_cols(A), mtl::matrix_not_square());
+	mtl::vampir_trace<5038> tracer;
+
+	typedef typename mtl::Collection<Matrix>::value_type      value_type;
+	typedef typename mtl::Collection<Matrix>::size_type       size_type;
+	typedef mtl::matrix::parameters<mtl::row_major, mtl::index::c_index, mtl::non_fixed::dimensions, false, size_type> para;
+	typedef mtl::matrix::compressed2D<value_type, para>  LU_type;
+	LU_type LU(A);
+
+	typedef typename range_generator<row, LU_type>::type      cur_type;
+	typedef typename range_generator<nz, cur_type>::type      icur_type;
+	typename mtl::traits::col<LU_type>::type                  col(LU);
+	typename mtl::traits::value<LU_type>::type                value(LU);
+	mtl::vector::dense_vector<value_type, mtl::vector::parameters<> > inv_dia(num_rows(A));
+	cur_type ic= begin<row>(LU), iend= end<row>(LU);
+	for (size_type i= 0; ic != iend; ++ic, ++i) {
+
+	    for (icur_type kc= begin<nz>(ic), kend= end<nz>(ic); kc != kend; ++kc) {
+		size_type k= col(*kc);
+		if (k == i) break;
+
+		value_type aik= value(*kc) * inv_dia[k];
+		value(*kc, aik);
+
+		for (icur_type jc= kc + 1; jc != kend; ++jc)
+		    value(*jc, value(*jc) - aik * LU[k][col(*jc)]);
+		// std::cout << "LU after eliminating A[" << i << "][" << k << "] =\n" << LU;
+	    }
+	    inv_dia[i]= reciprocal(LU[i][i]);
+	}
+	invert_diagonal(LU);
+	L= strict_lower(LU);
+	U= upper(LU);
+    }
+};
+
+template <typename Matrix, typename Value= typename mtl::Collection<Matrix>::value_type>
+class ilu_0
+  : public ilu<Matrix, ilu_0_factorizer, Value>
+{
+    typedef ilu<Matrix, ilu_0_factorizer, Value> base;
+  public:
+    ilu_0(const Matrix& A) : base(A) {}
+};
+
+// ic_0_evaluator not needed IC(0) and ILU(0) do the same at the upper triangle ;-)
+
+}} // namespace itl::pc
+
+
+#endif // ITL_PC_ILU_0_INCLUDE

AMDiS/lib/mtl4/boost/numeric/itl/pc/ilut.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef ITL_PC_ILUT_INCLUDE
+#define ITL_PC_ILUT_INCLUDE
+
+#include <boost/numeric/mtl/vector/sparse_vector.hpp>
+#include <boost/numeric/mtl/operation/invert_diagonal.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+namespace itl { namespace pc {
+
+
+struct ilut_factorizer
+{
+
+    template <typename Matrix, typename Para, typename L_type, typename U_type>
+    ilut_factorizer(const Matrix &A, const Para& p, L_type& L, U_type& U)
+    { factorize(A, p, L, U, mtl::traits::is_row_major<Matrix>()); }
+
+    // column-major matrices are copied first
+    template <typename Matrix, typename Para, typename L_type, typename U_type, bool B>
+    void factorize(const Matrix &A, const Para& p, L_type& L, U_type& U, boost::mpl::bool_<B>)
+    {
+	typedef typename mtl::Collection<Matrix>::value_type      value_type;
+	typedef typename mtl::Collection<Matrix>::size_type       size_type;
+	typedef mtl::matrix::parameters<mtl::row_major, mtl::index::c_index, mtl::non_fixed::dimensions, false, size_type> para;
+	typedef mtl::matrix::compressed2D<value_type, para>  LU_type;
+	LU_type LU(A);
+	factorize(LU, p, L, U, boost::mpl::true_());
+    }
+
+    // According Yousef Saad: ILUT, NLAA, Vol 1(4), 387-402 (1994)
+#if 0
+    template <typename Value, typename MPara, typename Para, typename L_type, typename U_type>
+    factorize(const mtl::matrix::compressed2D<Value, MPara>& A, const Para& p, L_type& L, U_type& U, boost::mpl::true_)
+#endif
+
+    template <typename Matrix, typename Para, typename L_type, typename U_type>
+    void factorize(const Matrix& A, const Para& p, L_type& L, U_type& U, boost::mpl::true_)
+
+    {
+	mtl::vampir_trace<5049> tracer;
+	using std::abs; using mtl::traits::range_generator; using mtl::begin; using mtl::end;
+	using namespace mtl::tag;
+	MTL_THROW_IF(num_rows(A) != num_cols(A), mtl::matrix_not_square());
+
+	typedef typename mtl::Collection<Matrix>::value_type      value_type;
+	typedef typename mtl::Collection<Matrix>::size_type       size_type;
+	typedef typename range_generator<row, Matrix>::type       cur_type;
+	typedef typename range_generator<nz, cur_type>::type      icur_type;
+	typename mtl::traits::col<Matrix>::type                   col(A);
+	typename mtl::traits::const_value<Matrix>::type           value(A);
+
+	size_type n= num_rows(A);
+	L.change_dim(n, n);
+	U.change_dim(n, n);
+	{
+	    mtl::matrix::inserter<L_type> L_ins(L, p.first);
+	    mtl::matrix::inserter<U_type> U_ins(U, p.first + 1); // plus one for diagonal
+
+	    mtl::vector::sparse_vector<value_type> vec(n); // corr. row in paper
+	    cur_type ic= begin<row>(A); // , iend= end<row>(A);
+	    for (size_type i= 0; i < n; ++i, ++ic) {
+
+		for (icur_type kc= begin<nz>(ic), kend= end<nz>(ic); kc != kend; ++kc) // row= A[i][*]
+		    vec.insert(col(*kc), value(*kc));
+		// std::cerr << "vec_" << i << " = " << vec << std::endl;
+		value_type tau_i= p.second * two_norm(vec); // threshold for i-th row
+		// loop over non-zeros in vec; changes in vec considered
+		for (size_type j= 0; j < vec.nnz() && vec.index(j) < i; j++) {
+		    size_type k= vec.index(j);
+		    value_type ukk= U_ins.value(k, k);
+		    MTL_DEBUG_THROW_IF(ukk == value_type(0), mtl::missing_diagonal());
+		    value_type vec_k= vec.value(j)/= ukk;
+		    // std::cout << "vec after updating from U[" << k << "][" << k << "] is " << vec << '\n';
+		    for (size_type j0= U_ins.ref_major()[k], j1= U_ins.ref_slot_ends()[k]; j0 < j1; j0++) { // U[k][k+1:n]
+			size_type k1= U_ins.ref_minor()[j0];
+			if (k1 > k)
+			    vec[k1]-= vec_k * U_ins.ref_elements()[j0];
+			// std::cout << "vec after updating from U[" << k << "][" << k1 << "] is " << vec << '\n';
+		    }
+		    // if (i > 1000 && i < 1010) std::cout << "vec before crop in row " << i << ", updating from row " << k << ": \n" << vec << "\n";
+		    vec.crop(tau_i);
+		    // if (i > 1000 && i < 1010) std::cout << "vec after crop: \n" << vec << "\n";
+		}
+		// std::cerr << "vec_" << i << " = " << vec << std::endl;
+		vec.sort_on_data();
+		// std::cerr << "vec_" << i << " sorted on data = " << vec << std::endl;
+
+		// std::cout << "vec at " << i << " is " << vec << '\n';
+		// mtl::vampir_trace<9904> tracer2;
+		bool diag_found= false;
+		for (size_type cntu= 0, cntl= 0, j= 0; j < vec.nnz() && (cntu < p.first || cntl < p.first); j++) {
+		    size_type k= vec.index(j);
+		    value_type v= vec.value(j);
+		    // if (abs(v) < tau_i) break;
+		    if (i == k) {
+			U_ins[i][i] << v; diag_found= true;
+		    } else if (i < k) {
+			if (cntu++ < p.first)
+			    U_ins[i][k] << v;
+		    } else // i > k
+			if (cntl++ < p.first)
+			    L_ins[i][k] << v;
+		}
+		if (!diag_found) std::cerr << "Deleted diagonal!!!!\n";
+		vec.make_empty();
+	    }
+	} // destroy inserters
+	invert_diagonal(U);
+    }
+};
+
+// Not usable yet !!!!!
+template <typename Matrix, typename Value= typename mtl::Collection<Matrix>::value_type>
+class ilut
+  : public ilu<Matrix, ilut_factorizer, Value>
+{
+    typedef ilu<Matrix, ilut_factorizer, Value> base;
+  public:
+    ilut(const Matrix& A, std::size_t p, typename mtl::Collection<Matrix>::value_type tau)
+      : base(A, std::make_pair(p, tau)) {}
+};
+
+}} // namespace itl::pc
+
+#endif // ITL_PC_ILUT_INCLUDE

AMDiS/lib/mtl4/boost/numeric/itl/pc/imf_algorithms.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+//
+// Algorithm inspired by Nick Vannieuwenhoven, written by Cornelius Steinhardt
+
+/*
+ * Implementation of the IMF factorization and application routines.
+ */
+
+#ifndef MTL_IMF_ALGORITHMS_INCLUDE
+#define MTL_IMF_ALGORITHMS_INCLUDE
+
+#include <iostream>
+#include <string.h>
+#include <vector>
+#include <map>
+#include <set>
+#include <complex>
+#include <limits>
+
+#include <boost/numeric/mtl/interface/vpt.hpp>
+#include <boost/numeric/mtl/matrix/compressed2D.hpp>
+#include <boost/numeric/mtl/matrix/coordinate2D.hpp>
+#include <boost/numeric/mtl/operation/clone.hpp>
+#include <boost/numeric/mtl/utility/irange.hpp>
+#include <boost/numeric/mtl/utility/make_copy_or_reference.hpp>
+
+#include <boost/numeric/itl/pc/binary_heap.hpp>
+#include <boost/numeric/itl/pc/sorting.hpp>
+
+#include "boost/unordered_map.hpp"
+#include "boost/unordered_set.hpp"
+
+
+namespace itl {   namespace pc {
+////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+//
+// INTRUSIVE HEAP HELPER FUNCTORS
+//
+////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+
+namespace heap {
+/**
+ * An information extending structure for elements to support the reduction
+ * process.
+ */
+template<class Element, class Key>
+struct HeapInformation {
+	Key m_key;
+	Element* m_parent;
+	Element* m_left_child;
+	Element* m_right_child;
+
+	HeapInformation() :
+		m_key(0.0), m_parent(0), m_left_child(0), m_right_child(0) {}
+};
+
+/**
+ * A functor returning a reference to the degree of an element.
+ */
+template<class Element, class Key>
+struct GetKey {
+	Key& operator()(Element* element) const {
+		return static_cast<HeapInformation<Element,Key>* >(
+			element->get_extra_pointer()
+		)->m_key;
+	}
+};
+
+/**
+ * A functor returning a reference to the parent of an element (in a heap).
+ */
+template<class Element, class Key>
+struct GetParent {
+	Element*& operator()(Element* element) const {
+		return static_cast<HeapInformation<Element,Key>* >(
+			element->get_extra_pointer()
+		)->m_parent;
+	}
+};
+
+/**
+ * A functor returning a reference to the left child of an element (in a heap).
+ */
+template<class Element, class Key>
+struct GetLeftChild {
+	Element*& operator()(Element* element) const {
+		return static_cast<HeapInformation<Element,Key>* >(
+			element->get_extra_pointer()
+		)->m_left_child;
+	}
+};
+
+/**
+ * A functor returning a reference to the right child of an element (in a heap).
+ */
+template<class Element, class Key>
+struct GetRightChild {
+	Element*& operator()(Element* element) const {
+		return static_cast<HeapInformation<Element,Key>* >(
+			element->get_extra_pointer()
+		)->m_right_child;
+	}
+};
+
+/**
+ * Compares two elements based on their degrees.
+ */
+template<class Element, class Key>
+struct DegreeCompare {
+	/**
+	 * Compares two elements based on their degree. Ties are broken by the
+	 * sequence numbers of the elements.
+	 */
+	inline bool operator()(Element* first, Element* second) const {
+		const int seq_fst = first->get_id();
+		const int seq_snd = second->get_id();
+
+		const Key key_fst = static_cast<HeapInformation<Element,Key>* >(
+				first->get_extra_pointer()
+			)->m_key;
+		const Key key_snd = static_cast<HeapInformation<Element,Key>* >(
+				second->get_extra_pointer()
+			)->m_key;
+
+		if (key_fst == key_snd) {
+			return seq_fst < seq_snd;
+		}
+		return key_fst < key_snd;
+	}
+};
+
+} // end namespace heap
+
+
+enum Status { UNMARKED, DIAGONAL, MARKED_CURRENT, REMOVED, NON_DIAGONAL };
+
+
+////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+//
+// DIAGONAL BLOCK SELECTION PRIORITY ESTIMATORS
+//
+////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+
+/**
+ * Estimates the priority based on the number of nodes the element is connected
+ * to for the purpose of the local Schur complement computation. A higher
+ * priority is given if the element is connected to fewer nodes.
+ */
+template< class Element, class NodeStatusVector, bool UseStatus = false >
+struct MinConnectedNodesEstimation {
+	inline int operator()(
+			const Element& el,
+			const NodeStatusVector& status
+	) const {
+		typedef typename Element::neighbor_collection_type neigh_type;
+
+		// Determine set of all nodes.
+		std::vector<int> nodes;
+		const neigh_type& neighs = el.get_neighbors();
+		for(int i = 0; i < el.get_nb_neighbors(); ++i) {
+			nodes.insert(
+				nodes.end(),
+				neighs[i]->get_indices().begin(),
+				neighs[i]->get_indices().end()
+			);
+		}
+
+//		radix_sort( &nodes[0], nodes.size() );   // TODO INCLUDE RADIX_SORT
+		std::sort( nodes.begin(), nodes.end() );
+
+		int degree = -el.nb_vars()+1;
+		for(unsigned int i = 1; i < nodes.size(); ++i)
+			if( nodes[i-1] != nodes[i] )
+			    if (!UseStatus || status[nodes[i]] == UNMARKED)
+				++degree;
+		return degree;
+	}
+};
+
+
+
+////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+//
+// IMF Factorization
+//
+////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+
+
+template<class Type>
+struct AddressHasher {
+	long operator()(const Type*& object) const {
+		return reinterpret_cast<long>( object );
+	}
+};
+
+template<class Element, class StatusVector>
+struct IsRemoved {
+	const StatusVector& m_status;
+	IsRemoved(const StatusVector& status) : m_status(status) { }
+
+	bool operator()(const Element *const el) const {
+		return m_status[ el->get_id() ] == REMOVED;
+	}
+};
+
+/**
+ * Constructs the EBE-ML-ILU preconditioner from the given mesh. The mesh is
+ * altered in the process.
+ */
+template< typename ValType >
+template<  class Mesh >
+void itl::pc::imf_preconditioner<ValType>::factor(const Mesh& mesh , const int maxlofi
+) {
+	mtl::vampir_trace<5052> tracer;
+	typedef typename Mesh::element_type element_type;
+	typedef typename Mesh::element_iterator element_iterator;
+	typedef typename element_type::value_type value_type;
+	typedef typename element_type::neighbor_collection_type neigh_coll_type;
+	typedef typename element_type::neighbor_iterator neigh_iterator;
+	typedef typename element_type::neighbor_set_type neigh_set_type;
+	typedef typename element_type::neighbor_set_iterator_type neigh_set_iterator;
+	// typedef typename neigh_coll_type::const_iterator const_neigh_iterator;
+	typedef typename element_type::index_type index_type;
+	typedef typename element_type::matrix_type matrix_type;
+
+	typedef typename mtl::matrix::coordinate2D<value_type> coo_sparse_type_upper;
+	typedef typename mtl::matrix::coordinate2D<value_type> coo_sparse_type_lower;
+
+
+	typedef std::map<int, int> cmap;
+	// typedef typename cmap::iterator cmap_iterator;
+
+	typedef value_type key_type;
+	typedef utils::binary_heap<
+		element_iterator,
+		key_type,
+		heap::DegreeCompare<element_type, key_type>,
+		element_type*,
+		heap::GetKey<element_type, key_type>,
+		heap::GetParent<element_type, key_type>,
+		heap::GetLeftChild<element_type, key_type>,
+		heap::GetRightChild<element_type, key_type>
+	> my_heap;
+
+
+	// Constants.
+	const int nb_elements = mesh.get_total_elements();
+	const int nb_vars = mesh.get_total_vars();
+	const int UNPERMUTED = -1;
+	const value_type zero(0);
+	// A DS tracking the set of elements during the construction.
+	//	Invariant: elements[i]->get_sequence_number() == i
+	//	Invariant: elements[i] == 0  iff  element[i] is removed
+	std::vector<element_type*> elements;
+	elements.reserve( nb_elements + (nb_elements >> 4) );
+	element_iterator it = mesh.element_begin();
+	for (int i = 0; i < nb_elements; ++i) {
+	  elements.push_back( *&it );
+	  ++it;
+	}
+
+
+	// Data structures for the preconditioner.
+	std::vector<element_type*> block_diagonal;
+	std::vector<coo_sparse_type_upper*> upper_matrices;
+	std::vector<coo_sparse_type_lower*> lower_matrices;
+	std::vector<int> diagonal_offsets;
+	diagonal_offsets.push_back(0);
+	m_ordering = UNPERMUTED;
+
+	////////////////////////////////////////////////////////////////////////////
+	// Auxiliary data structures.
+	////////////////////////////////////////////////////////////////////////////
+
+		// Binary heap containing the dense elements that should still be
+		// considered for selection on the diagonal of the current level.
+		heap::DegreeCompare<element_type, key_type> key_compare;
+		heap::GetKey<element_type, key_type> get_key;
+		heap::GetParent<element_type, key_type> get_parent;
+		heap::GetLeftChild<element_type, key_type> get_left;
+		heap::GetRightChild<element_type, key_type> get_right;
+		my_heap unmarked_elements(
+				key_compare,
+				get_key,
+				get_parent,
+				get_left,
+				get_right
+		);
+		std::vector< heap::HeapInformation<element_type, key_type>* >	red_info( nb_elements );
+
+		// Another data structure for the elements that should be considered
+		// for the diagonal
+		std::vector<element_type*> unmarked_elements_srtd;
+		std::vector<int          > unmarked_elements_degr;
+		unmarked_elements_srtd.reserve( nb_elements );
+		unmarked_elements_degr.reserve( nb_elements );
+
+		// A DS to keep track of whether the element is on the diagonal, not
+		// on the diagonal or not yet marked.
+		std::vector<Status> el_status(nb_elements, UNMARKED);
+
+		// A DS to keep track of the variables that have been forced into the
+		// reduced system.
+		std::vector<Status> in_reduced( mesh.get_total_vars(), UNMARKED );
+
+	////////////////////////////////////////////////////////////////////////////
+	////////////////////////////////////////////////////////////////////////////
+
+	// Degree estimators.
+	const MinConnectedNodesEstimation<element_type, std::vector<Status> > min_nodes_estimate =
+	      MinConnectedNodesEstimation<element_type, std::vector<Status> >();
+
+	std::cout << mesh.get_total_vars() << " nodes remaining" << std::endl;
+
+	int level = 1;
+	int perm_off = 0;
+	int perm_high = 0;
+	int block_diag_low = 0;
+	int max_sequence_number = nb_elements;
+
+ 	do {
+	    mtl::vampir_trace<9901> tb1;
+
+		block_diag_low = block_diagonal.size();
+
+/***************************************************************************
+ * Phase 0: Determine priority of unmarked elements
+ **************************************************************************/
+
+		// Make sure the set of unmarked elements is empty again.
+		assert( unmarked_elements.empty() );
+		unmarked_elements_degr.clear();
+		unmarked_elements_srtd.clear();
+		// Construct the set of available diagonal elements along with their
+		// degrees.
+		for( unsigned int i = 0; i < elements.size(); ++i ) {
+			// Skip the element if it no longer exists, or it is a diagonal
+			// element.
+			if( (el_status[i] == REMOVED) || (el_status[i] == DIAGONAL) ) {
+				continue;
+			}
+
+			element_type& el = *(elements[i]);
+			// Reset the status.
+			el_status[i] = UNMARKED;
+
+			int degree = min_nodes_estimate(el, in_reduced);
+		        unmarked_elements_degr.push_back( degree );
+			unmarked_elements_srtd.push_back( &el );
+		}
+
+/***************************************************************************
+ * Phase 1: Select the block diagonal elements
+ **************************************************************************/
+/********************************SKIP UPDATE***********************************/
+
+		// Sort the candidate diagonal elements by their degree in ascending
+		// order.
+		sort_along<int, element_type*>(
+			&*(unmarked_elements_degr.begin()),
+			&*(unmarked_elements_srtd.begin()),
+			unmarked_elements_degr.size()
+		);
+
+		// For each of the candidate diagonal elements ...
+		for(unsigned int i = 0; i < unmarked_elements_srtd.size(); ++i) {
+			// Select the minimum element.
+			element_type& el = *(unmarked_elements_srtd[i]);
+			// If the element is already marked, skip it.
+			if( el_status[el.get_id()] != UNMARKED ) {
+				continue;
+			}
+
+			// If the element is unmarked, add it to the set of diagonal
+			// elements
+			block_diagonal.push_back(&el);
+			// Mark the element.
+			el_status[el.get_id()] = DIAGONAL;
+
+			// Update permutation vector
+			for(int i = 0; i < el.nb_vars(); ++i) {
+				m_ordering( el.get_indices()(i) ) = perm_off;
+				++perm_off;
+				assert(perm_off <= mesh.get_total_vars());
+			}
+
+			// Determine level-2 neighbors.
+			neigh_set_type lvl2_neighs = el.get_level_neighbors( 2 );
+
+			// Mark all level-2 neighbors.
+			for(neigh_set_iterator neigh_it = lvl2_neighs.begin(); neigh_it != lvl2_neighs.end(); ++neigh_it) {
+				element_type& neigh = **neigh_it;
+				assert( el_status[neigh.get_id()] != DIAGONAL );
+				el_status[ neigh.get_id() ] = MARKED_CURRENT;
+			}
+		}
+
+
+/*******************************UPDATE DEGREE**********************************/
+		// Update diagonal block offset.
+		diagonal_offsets.push_back( perm_off );
+		//save upperbound for number of L and U entrys
+		unsigned int upperbound(0);
+		for(unsigned int i=0;i< block_diagonal.size();i++){
+			mtl::vector::dense_vector<int> involve_node(block_diagonal[i]->get_indices());
+			for(unsigned int j=0;j< block_diagonal[i]->get_neighbors().size();j++){
+			  mtl::vector::dense_vector<int> involve_neigh(block_diagonal[i]->get_neighbors()[j]->get_indices());
+			  unsigned int c(0);
+			  for(unsigned int a= 0; a < size(involve_node); a++){
+			      for(unsigned int b= 0; b < size(involve_neigh); b++){
+				if(involve_node[a] == involve_neigh[b])
+				  c++;
+			      }
+			  }
+			  upperbound+= c*(size(involve_neigh)-c);
+			}
+		}
+
+		// Update permutation offsets.
+		perm_high = perm_off;
+
+/***************************************************************************
+ * Phase 2: Compute the update matrices
+ **************************************************************************/
+		mtl::vampir_trace<9902> tb2;
+
+
+		coo_sparse_type_lower* L=new coo_sparse_type_lower(nb_vars, nb_vars, upperbound);
+		coo_sparse_type_upper* U=new coo_sparse_type_lower(nb_vars, nb_vars, upperbound);
+		lower_matrices.push_back(L);
+		upper_matrices.push_back(U);
+		// For each diagonal block element ... (in parallel)
+		unsigned int ku= 0, kl= 0;
+		for(std::size_t b_i = block_diag_low; b_i < block_diagonal.size();++b_i ) {
+			element_type& diag_el = *block_diagonal[b_i];
+			// Copy the level-1 neighbors.
+			neigh_coll_type& diag_neighs = diag_el.get_neighbors();
+
+			// Determine the set of incident nodes.
+			boost::unordered_set<int> diag_incident_nodes =
+				diag_el.get_incident_nodes();
+			index_type& p = diag_el.get_indices();
+			index_type q(
+				diag_incident_nodes.size() == 0 ?
+					1 : diag_incident_nodes.size()
+			);
+			typename boost::unordered_set<int>::const_iterator it =
+				diag_incident_nodes.begin();
+			for(unsigned int i = 0; i < diag_incident_nodes.size(); ++i) {
+				q(i) = *it;
+				++it;
+			}
+			const int n1 = size(p);
+			const int n2 = diag_incident_nodes.size();
+			sort(q);
+			assert(n1 > 0);
+
+			// Construct a mapping from global to local node numbers.
+			cmap to_local;
+			for(int i = 0; i < n1; ++i) {
+				to_local[p(i)] = i;
+			}
+
+			for(int i = 0; i < n2; ++i) {
+				to_local[q(i)] = mtl::size(p) + i;
+			}
+			// Construct the frontal matrix.
+			block_type frontal( n1+n2, n1+n2 );
+			frontal = zero;
+
+			// For each connected neighbor, add their values to the frontal
+			// matrix.
+			for(neigh_iterator neigh_it = diag_neighs.begin(); neigh_it != diag_neighs.end(); ++neigh_it) {
+				element_type& neigh = **neigh_it;
+				assert( el_status[neigh.get_id()] != DIAGONAL );
+				assert( el_status[neigh.get_id()] != REMOVED );
+
+				// Remap indices.
+				index_type local_idx( neigh.get_indices() );
+
+				for(int i = 0; i < neigh.nb_vars(); ++i) {
+					local_idx(i) = to_local[ local_idx(i) ];
+				}
+				//insert connectet neighbor
+				{
+				  mtl::matrix::inserter<mtl::matrix::dense2D<value_type>, mtl::operations::update_plus<value_type> > ins(frontal);
+				  ins << element_matrix(neigh.get_values(), local_idx, local_idx);
+				}
+				neigh.get_values()= zero;
+			}
+			// Add the values of the diagonal element to the frontal matrix.
+			{
+				// Remap indices.
+				index_type local_idx( diag_el.get_indices() );
+				for(int i = 0; i < diag_el.nb_vars(); ++i) {
+					local_idx(i) = to_local[ local_idx(i) ];
+				}
+				//insert the diagonal element
+				{
+				  mtl::matrix::inserter<matrix_type, mtl::operations::update_plus<value_type> > ins(frontal);
+				  ins << element_matrix(
+				  diag_el.get_values(), local_idx, local_idx);
+				}
+				diag_el.get_values()= zero;
+			}
+			//
+			// Store the L and U part.
+			//
+			for(int i = 0; i < n1; ++i) { // p-part
+				for(int j = n1; j < n1+n2; ++j) { // q-part
+				  // Assumption: structural symmetry.
+					if(frontal(i,j) != zero) {
+						U->insert(p(i),q(j-n1),frontal(i,j));
+						ku++;
+					}
+					if(frontal(j,i) != zero){
+						L->insert(q(j-n1),p(i),frontal(j,i));
+						kl++;
+					}
+				}
+			}
+
+			mtl::irange  n0(0,n1);
+    			diag_el.get_values() = inv(mtl::clone(frontal[n0][n0]));
+			frontal[n0][n0] = diag_el.get_values();
+
+ 			if( (level <= maxlofi) && (n2 > 0) ) {
+//------------------------------------------------------------------------------
+//  					ELEMENT COALESCING APPROACH
+//------------------------------------------------------------------------------
+				// The lofi is below the user-requested level. Use Algorithm 2
+				// in [1].
+				// Compute the Schur complement.
+				mtl::irange n1r(0,n1), n2r(n1, mtl::imax);
+				frontal[n2r][n2r]-= frontal[n2r][n1r] * frontal[n1r][n1r] * frontal[n1r][n2r];
+				mtl::irange nz(n1,n1+n2);
+				matrix_type Z( mtl::clone(frontal[nz][nz]) );
+				// Construct the new element.
+				element_type* fill = new element_type(max_sequence_number, q, Z);
+				// Add processing information to the required DSs.
+				elements.push_back( fill );
+				el_status.push_back( UNMARKED );
+
+
+				// Determine the level-2 neighbors.
+				neigh_set_type lvl2_neighs = diag_el.get_level_neighbors( 2 );
+				// Remove the level-1 neighbors of the diagonal element.
+				for(neigh_iterator it = diag_neighs.begin(); it != diag_neighs.end(); ++it) {
+					element_type& neigh = **it;
+					neigh.clear();
+					el_status[neigh.get_id()] = REMOVED;
+				}
+
+				// Update the neighborhood of the level-2 neighbors.
+				const IsRemoved<element_type, std::vector<Status> > is_removed( el_status );
+				for(neigh_set_iterator it = lvl2_neighs.begin(); it != lvl2_neighs.end(); ++it) {
+					element_type& neigh = **it;
+
+					// Skip the level-1 neighbors (removed) and the diagonal
+					// element.
+					if((el_status[neigh.get_id()] == DIAGONAL) || (el_status[neigh.get_id()] == REMOVED)) {
+						continue;
+					}
+					// The element is in the strict level-2 neighborhood.
+					// Remove the level-1 neighbors from its set of neighbors.
+					neigh_coll_type& neigh_neighs = neigh.get_neighbors();
+					neigh_iterator new_end = std::remove_if(neigh_neighs.begin(), neigh_neighs.end(), is_removed);
+					neigh_neighs.erase(new_end, neigh_neighs.end());
+
+					// Add the newly generated element as neighbor, and vice
+					// versa.
+					neigh_neighs.push_back( fill );
+					fill->get_neighbors().push_back( &neigh );
+
+				}
+				// Update sequence number.
+				++max_sequence_number;
+			} else if (n2 > 0) {
+
+//------------------------------------------------------------------------------
+//  					ELEMENT DISTRIBUTION APPROACH
+//------------------------------------------------------------------------------
+				// The lofi is above the user-requested level. Use Algorithm 3
+				// in [1] to compute the approximate Schur complement element-
+				// wise.
+
+				// Distribute the values of the generalized element across the
+				// level-2 neighbors of the diagonal element.
+
+				// Remove the nodes of the diagonal element from its
+				// neighbors.
+				for(neigh_iterator neigh_it = diag_neighs.begin(); neigh_it != diag_neighs.end(); ++neigh_it)
+				{
+					element_type& neigh = **neigh_it;
+					neigh.remove_nodes( diag_el.get_indices(), diag_el );
+
+					assert( el_status[neigh.get_id()] != DIAGONAL );
+					assert( el_status[neigh.get_id()] != REMOVED );
+
+					// If the element is entirely removed, mark it as such.
+					if( neigh.nb_vars() == 0 ) {
+						el_status[neigh.get_id()] =	REMOVED;
+					}
+				}
+
+				// Compute the Schur complement.
+				mtl::irange n1r(0,n1), n2r(n1, mtl::imax);
+				frontal[n2r][n2r]-= frontal[n2r][n1r] * frontal[n1r][n1r] * frontal[n1r][n2r];
+				// Distribute the values of the (modified) update matrix
+				// over the level-1 neighbors.
+				mtl::irange nz(n1,n1+n2);
+  				matrix_type S( frontal[nz][nz] );
+				for(neigh_iterator neigh_it = diag_neighs.begin(); neigh_it != diag_neighs.end();++neigh_it) {
+					element_type& el = **neigh_it;
+					if( el_status[el.get_id()] != REMOVED ) {
+						el.absorb(S, q);
+					}
+				}
+				// Do not distribute over the entire level-2 neighborhood.
+				// This is expensive, while not adding much in terms of quality
+				// of the approximation.
+
+			} else {
+				// There are no more nodes in the Schur complement, but it could
+				// be that some elements other than the diagonal element overlap
+				// completely with said element.
+
+				// In this case, simply remove all neighbors.
+				for(neigh_iterator neigh_it = diag_neighs.begin(); neigh_it != diag_neighs.end(); ++neigh_it) {
+					element_type& neigh = **neigh_it;
+					neigh.remove_nodes( diag_el.get_indices(), diag_el );
+					assert( el_status[neigh.get_id()] != DIAGONAL );
+					assert( el_status[neigh.get_id()] != REMOVED );
+					// The element is now entirely removed, mark it as such.
+					assert( neigh.nb_vars() == 0 );
+					el_status[neigh.get_id()] = REMOVED;
+				}
+
+			} // END ELEMENT DISTRIBUTION
+			// Clear the neighborhood of the diagonal element.
+			diag_el.get_neighbors().clear();
+		}
+
+		std::cout << "level " << level << " complete: ";
+		std::cout << mesh.get_total_vars()-perm_high << " nodes remaining."<< "\n";
+
+		++level;
+	} while( (mesh.get_total_vars() - perm_high > 0) );	// There are still variables to cover.
+
+
+	diagonal_offsets.push_back( mesh.get_total_vars() );
+	assert( mesh.get_total_vars() - perm_high == 0 );
+	assert( lower_matrices.size()+2 == diagonal_offsets.size() );
+	/***************************************************************************
+	 * Phase 3: Apply permutation vector to lower and upper matrices
+	 **************************************************************************/
+	mtl::vampir_trace<9903> tb3;
+
+
+	mtl::matrix::traits::permutation<>::type P(permutation(m_ordering));
+	typedef typename coo_sparse_type_lower::size_type size_type;
+
+	for( std::size_t k = 0; k < lower_matrices.size(); ++k ) {
+ 		coo_sparse_type_lower& L = *(lower_matrices[k]);
+ 		coo_sparse_type_upper& U = *(upper_matrices[k]);
+ 		if (nnz(L)> 0) {
+		    std::vector<size_type> &L_row(L.row_index_array()),
+			                      &L_col(L.column_index_array()),
+					      &U_row(U.row_index_array()),
+					      &U_col(U.column_index_array());
+		    const int off_low = diagonal_offsets[k];
+		    for(unsigned int i = 0; i < nnz(L); ++i) {
+			L_row[i] = m_ordering( L_row[i] );
+			L_col[i] = m_ordering( L_col[i] ) - off_low;
+		    };
+		    for(unsigned int i = 0; i < nnz(U); ++i) {
+			  U_row[i] = m_ordering( U_row[i] ) - off_low;
+			  U_col[i] = m_ordering( U_col[i] );
+		    }
+		}
+
+		m_lower.push_back(mtl::matrix::compressed2D<value_type>(L));
+		m_upper.push_back(mtl::matrix::compressed2D<value_type>(U));
+	}
+	/***************************************************************************
+	 * Phase 4: Construct the IMF preconditioner
+	 **************************************************************************/
+	// Copy the block diagonal values into a consecutive array.
+//	mtl::vampir_trace<9904> tb4;
+
+	m_diagonal = new matrix_type[ block_diagonal.size() ];
+	for(std::size_t i = 0; i < block_diagonal.size(); ++i) {
+		assert( el_status[block_diagonal[i]->get_id()] == DIAGONAL );
+ 		unsigned int diarows(num_rows(block_diagonal[i]->get_values()));
+  		m_diagonal[i].change_dim(diarows,diarows);
+		m_diagonal[i] = block_diagonal[i]->get_values();
+	}
+	// Copy the diagonal offsets to a consecutive array.
+	int* diagonal_index = new int[ diagonal_offsets.size() ];
+	memcpy(
+		diagonal_index,
+		&(diagonal_offsets[0]),
+		sizeof(int)*diagonal_offsets.size()
+	);
+	assert( diagonal_offsets.back() == mesh.get_total_vars() );
+	assert( lower_matrices.size()+2 == diagonal_offsets.size() );
+
+	m_levels = lower_matrices.size();
+	m_nb_blocks = block_diagonal.size();
+	m_diagonal_index = diagonal_index;
+
+	// Todo: create element_structure from all elements
+	// Delete only elements that we generated; those at the beginning are just referred
+	for (unsigned i= nb_elements; i < elements.size(); i++)
+	  delete elements[i];
+
+	for (unsigned i= 0; i < lower_matrices.size(); i++)
+	  delete lower_matrices[i];
+	for (unsigned i= 0; i < upper_matrices.size(); i++)
+	  delete upper_matrices[i];
+}
+
+////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+//
+// IMF APPLICATION
+//
+////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+
+
+/**
+ * Applies the IMF preconditioner to a matrix. The contents of the rhs vector
+ * are overwritten.
+ *
+ * The algorithm is implemented without recursion, as was already suggested by
+ * Y. Saad in [2].
+ */
+
+template< class ValType >
+template< class Vector >
+Vector imf_preconditioner<ValType>::imf_apply(const Vector& rhs) const
+{
+	using namespace mtl;
+	Vector res(rhs);
+	ValType zero(0);
+	// Forward elimination.
+	{
+          mtl::vampir_trace<9901> tb1;
+	  int b_off = 0;
+	  for(int level = 0; level < m_levels; ++level) {
+		const int off_low = m_diagonal_index[level];
+		const int off_high = m_diagonal_index[level+1];
+		const int n1 = off_high - off_low;
+		assert(off_low <= off_high);
+		// Compute dy = inv(D)*y
+		vector_type dy(n1);
+		dy = zero;
+		for(int off = off_low; off < off_high; ++b_off ) { //parallel
+			const int block_size = num_rows( m_diagonal[b_off] );
+			dy[mtl::irange(off-off_low, off-off_low+block_size)] = m_diagonal[b_off] * res[mtl::irange(off, off + block_size)];
+			off += block_size;
+		}
+		// Compute x = x - E*dy
+		Vector big(num_cols(m_lower[level]), ValType(0));
+		big[mtl::irange(0,n1)] = dy;
+ 		res -= m_lower[level] * big;
+	  }
+	}
+
+	// Backward elimination.
+	{
+	  mtl::vampir_trace<9902> tb2;
+	  int b_off = m_nb_blocks-1;
+	  for(int level = m_levels-1; level >= 0; --level) {
+
+		const int off_low = m_diagonal_index[level];
+		const int off_high = m_diagonal_index[level+1];
+		// y' = y - Fx
+//  		assert( m_upper[level] );
+
+		vector_type yp(m_upper[level] * res);
+		res[mtl::irange(off_low, off_high) ] -= yp[mtl::irange(0, off_high-off_low) ];
+		// y = inv(D)*y'
+		for(int off = off_high; off > off_low; --b_off ) {
+			const int block_size = num_rows(m_diagonal[b_off]);
+			assert(b_off >= 0);
+			assert(off-block_size >= off_low);
+
+			vector_type dy(m_diagonal[b_off] * res[irange(off-block_size, off)] );
+			res[irange(off-block_size, off)] = dy;
+			off -= block_size;
+		}
+	}
+	assert(b_off == -1);
+	}
+    return res;
+}
+
+
+} // end namespace pc
+} // end namespace itl
+
+
+#endif // MTL_IMF_ALGORITHMS_INCLUDE

AMDiS/lib/mtl4/boost/numeric/itl/pc/imf_preconditioner.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+//
+// Algorithm inspired by Nick Vannieuwenhoven, written by Cornelius Steinhardt
+
+/*
+ * The IMF preconditioner.
+ *
+ * References
+ * 	[1] N. Vannieuwenhoven and K. Meerbergen, IMF: An incomplete multifron-
+ *		tal LU-factorization for element-structured sparse linear systems,
+ *		Tech. Rep. TW581, Department of Computer Science, KULeuven,
+ *		December 2010.
+ */
+
+#ifndef MTL_IMF_PRECONDITIONER_INCLUDE
+#define MTL_IMF_PRECONDITIONER_INCLUDE
+
+#include <boost/numeric/itl/pc/matrix_algorithms.hpp>
+#include <boost/numeric/itl/pc/solver.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+#include <boost/numeric/mtl/matrix/compressed2D.hpp>
+#include <boost/numeric/mtl/matrix/coordinate2D.hpp>
+#include <boost/numeric/mtl/matrix/dense2D.hpp>
+#include <boost/numeric/mtl/vector/dense_vector.hpp>
+
+
+namespace itl {   namespace pc {
+
+/// The IMF preconditioner, as described in [1].
+template<class ValType>
+class imf_preconditioner {
+
+public:
+	/// The type of the values.
+	typedef ValType value_type;
+
+	/// The type of the sparse data structure for the upper matrices.
+	typedef mtl::matrix::coordinate2D<value_type> sparse_type_upper;
+
+	/// The type of the sparse data structure for the lower matrices.
+	typedef mtl::matrix::coordinate2D<value_type> sparse_type_lower;
+
+	/// The type of the vectors.
+	typedef mtl::vector::dense_vector<value_type> vector_type;
+
+	///The type of the permutation vector.
+	typedef mtl::vector::dense_vector<int>  index_type;
+
+	/// The type of the matrices on the block diagonal.
+	typedef mtl::matrix::dense2D<value_type>  block_type;
+
+	/// The type of the sequence of lower matrices.
+	typedef std::vector<mtl::matrix::compressed2D<value_type> > lower_matrix_coll_type;
+
+	/// The type of the sequence of upper matrices.
+	typedef std::vector<mtl::matrix::compressed2D<value_type> > upper_matrix_coll_type;
+
+	/// Constructor
+	template< class ElementStructure >
+	imf_preconditioner(
+		const ElementStructure& element_structure ,
+  		const int maxlofi=0,
+		const bool copy_on=true
+	)
+	: 	m_nb_vars( element_structure.get_total_vars() ),
+	  	m_ordering( element_structure.get_total_vars() ),
+	  	m_diagonal_index(0),
+	  	m_diagonal(0)
+	{
+ 		mtl::vampir_trace<5053> tracer;
+		if(copy_on){
+		  ElementStructure es(element_structure);
+		  factor(es, maxlofi);
+		} else {
+		  factor(element_structure, maxlofi);
+		}
+		P= permutation(m_ordering);
+	}
+
+	/// Destructor
+	~imf_preconditioner()
+        {
+	  delete[] m_diagonal_index;
+	  delete[] m_diagonal;
+	}
+
+private:
+	/// Disallow the copy constructor and assignment.
+	imf_preconditioner();
+	imf_preconditioner(const imf_preconditioner& other);
+	void operator=(const imf_preconditioner& other);
+
+        /// Constructs the IMF preconditioner.Forward declaration
+	template< class ElementStructure >
+  	void factor(const ElementStructure&, const int);
+
+public:
+
+	/// Returns the number of levels (equals the number of lower and upper matrices.
+	int get_nb_levels() const { return m_levels; }
+
+	/// Returns the number of blocks on the diagonal.
+	int get_nb_blocks() const { return m_nb_blocks; }
+
+
+        /// Applies the preconditioner to the given matrix.
+        template <typename VectorIn, typename VectorOut>
+        void solve(const VectorIn& b, VectorOut& x) const
+        {
+	    VectorIn m(trans(P)*b), m_tmp(imf_apply(m));
+	    x= P * m_tmp;
+        }
+
+private:
+        /// Applies the preconditioner   prototype
+	template< class Vector >
+	Vector imf_apply(const Vector&) const;
+
+	/// The number of variables (also the size of the m_ordering vector).
+	unsigned int m_nb_vars;
+
+	/// The number of blocks on the diagonal.
+	unsigned int m_nb_blocks;
+
+	/// A vector containing the renumbering of IMF.
+	index_type m_ordering;
+
+	/// A matrix containing the renumbering of IMF.
+        mtl::matrix::traits::permutation<>::type P;
+
+	/// The number of levels (equals the number of entries in the diagonal index array minus one).
+	int m_levels;
+
+	 /** The index array for the matrices on the block diagonal. The i^th entry
+	 * indicates where the i^th level of block diagonal matrices starts in the
+	 * right hand side vector. */
+	int* m_diagonal_index;
+
+        /// The matrices on the block diagonal.
+	block_type* m_diagonal;
+
+	/// The sparse lower matrices of each level, sorted by level.
+	lower_matrix_coll_type m_lower;
+
+	/// The sparse upper matrices of each level, sorted by level.
+	upper_matrix_coll_type m_upper;
+};
+
+/// Solve
+template <typename Matrix, typename Vector>
+solver<imf_preconditioner<Matrix>, Vector, false>
+inline solve(const imf_preconditioner<Matrix>& P, const Vector& b)
+{
+    mtl::vpt::vampir_trace<5054> tracer;
+    return solver<imf_preconditioner<Matrix>, Vector, false>(P, b);
+}
+}//namespace pc
+}//namespace itl
+#endif // MTL_IMF_PRECONDITIONER_INCLUDE