Commit 395f1b3c authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

cleanup for updated dune-* modules

parent 3c2194bd
Pipeline #142 failed with stage
// 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_BASIC_ITERATION_INCLUDE
#define ITL_BASIC_ITERATION_INCLUDE
#include <iostream>
#include <complex>
#include <string>
namespace itl {
template <class Real>
class basic_iteration
{
public:
typedef basic_iteration self;
typedef Real real;
template <class Vector>
basic_iteration(const Vector& r0, int max_iter_, Real t, Real a = Real(0))
: error(0), i(0), norm_r0(std::abs(two_norm(r0))),
max_iter(max_iter_), rtol_(t), atol_(a), is_finished(false), my_quite(false), my_suppress(false) { }
basic_iteration(Real nb, int max_iter_, Real t, Real a = Real(0))
: error(0), i(0), norm_r0(nb), max_iter(max_iter_), rtol_(t), atol_(a), is_finished(false),
my_quite(false), my_suppress(false) {}
virtual ~basic_iteration() {}
bool check_max()
{
if (i >= max_iter)
error= 1, is_finished= true, err_msg= "Too many iterations.";
return is_finished;
}
template <class Vector>
bool finished(const Vector& r)
{
if (converged(two_norm(r)))
return is_finished= true;
return check_max();
}
bool finished(const Real& r)
{
if (converged(r))
return is_finished= true;
return check_max();
}
template <typename T>
bool finished(const std::complex<T>& r)
{
if (converged(std::abs(r)))
return is_finished= true;
return check_max();
}
bool finished() const { return is_finished; }
template <class T>
int terminate(const T& r) { finished(r); return error; }
bool converged(const Real& r) { resid_= r; return converged(); }
bool converged() const
{
if (norm_r0 == 0)
return resid_ <= atol_; // ignore relative tolerance if |r0| is zero
return resid_ <= rtol_ * norm_r0 || resid_ <= atol_;
}
self& operator++() { ++i; return *this; }
self& operator+=(int n) { i+= n; return *this; }
bool first() const { return i <= 1; }
virtual operator int() const { return error; }
virtual int error_code() const { return error; }
bool is_converged() const { return is_finished && error == 0; }
int iterations() const { return i; }
int max_iterations() const { return max_iter; }
void set_max_iterations(int m) { max_iter= m; }
Real resid() const { return resid_; }
Real relresid() const { return resid_ / norm_r0; }
Real normb() const { return norm_r0; }
Real tol() const { return rtol_; }
Real atol() const { return atol_; }
int fail(int err_code) { error = err_code; return error_code(); }
int fail(int err_code, const std::string& msg)
{ error = err_code; err_msg = msg; return error_code(); }
void set(Real v) { norm_r0 = v; }
void set_quite(bool q) { my_quite= q; }
bool is_quite() const { return my_quite; }
void suppress_resume(bool s) { my_suppress= s; }
bool resume_suppressed() const { return my_suppress; }
void update_progress(const basic_iteration& that)
{
i= that.i;
resid_= that.resid_;
if (that.error > 1) { // copy error except too many iterations
error= that.error;
err_msg= that.err_msg;
is_finished= true;
} else
finished(resid_);
}
protected:
int error, i;
Real norm_r0;
int max_iter;
Real rtol_, atol_, resid_;
std::string err_msg;
bool is_finished, my_quite, my_suppress;
};
} // namespace itl
#endif // ITL_BASIC_ITERATION_INCLUDE
// 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 Library2
//
// See also license.mtl.txt in the distribution.
#ifndef ITL_CYCLIC_ITERATION_INCLUDE
#define ITL_CYCLIC_ITERATION_INCLUDE
#include <iostream>
#include <boost/numeric/itl/iteration/basic_iteration.hpp>
namespace itl {
/// Class for iteration control that cyclically prints residual
template <class Real, class OStream = std::ostream>
class cyclic_iteration : public basic_iteration<Real>
{
typedef basic_iteration<Real> super;
typedef cyclic_iteration self;
void print_resid()
{
if (!this->my_quite && this->i % cycle == 0)
if (multi_print || this->i != last_print) { // Avoid multiple print-outs in same iteration
out << "iteration " << this->i << ": resid " << this->resid()
// << " / " << this->norm_r0 << " = " << this->resid() / this->norm_r0 << " (rel. error)"
<< std::endl;
last_print= this->i;
}
}
public:
template <class Vector>
cyclic_iteration(const Vector& r0, int max_iter_, Real tol_, Real atol_ = Real(0), int cycle_ = 100,
OStream& out = std::cout)
: super(r0, max_iter_, tol_, atol_), cycle(cycle_), last_print(-1), multi_print(false), out(out)
{}
cyclic_iteration(Real r0, int max_iter_, Real tol_, Real atol_ = Real(0), int cycle_ = 100,
OStream& out = std::cout)
: super(r0, max_iter_, tol_, atol_), cycle(cycle_), last_print(-1), multi_print(false), out(out)
{}
bool finished() { return super::finished(); }
template <typename T>
bool finished(const T& r)
{
bool ret= super::finished(r);
print_resid();
return ret;
}
inline self& operator++() { ++this->i; return *this; }
inline self& operator+=(int n) { this->i+= n; return *this; }
operator int() const { return error_code(); }
/// Whether the residual is printed multiple times in iteration
bool is_multi_print() const { return multi_print; }
/// Set whether the residual is printed multiple times in iteration
void set_multi_print(bool m) { multi_print= m; }
int error_code() const
{
if (!this->my_suppress)
out << "finished! error code = " << this->error << '\n'
<< this->iterations() << " iterations\n"
<< this->resid() << " is actual final residual. \n"
<< this->relresid() << " is actual relative tolerance achieved. \n"
<< "Relative tol: " << this->rtol_ << " Absolute tol: " << this->atol_ << '\n'
<< "Convergence: " << pow(this->relresid(), 1.0 / double(this->iterations())) << std::endl;
return this->error;
}
protected:
int cycle, last_print;
bool multi_print;
OStream& out;
};
} // namespace itl
#endif // ITL_CYCLIC_ITERATION_INCLUDE
// 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_NOISY_ITERATION_INCLUDE
#define ITL_NOISY_ITERATION_INCLUDE
#include <iostream>
#include <boost/numeric/itl/iteration/cyclic_iteration.hpp>
namespace itl {
template <class Real, class OStream = std::ostream>
class noisy_iteration : public cyclic_iteration<Real, OStream>
{
public:
template <class Vector>
noisy_iteration(const Vector& r0, int max_iter_, Real tol_, Real atol_ = Real(0),
OStream& out = std::cout)
: cyclic_iteration<Real, OStream>(r0, max_iter_, tol_, atol_, 1, out)
{}
};
} // namespace itl
#endif // ITL_NOISY_ITERATION_INCLUDE
// 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_ITL_INCLUDE
#define ITL_ITL_INCLUDE
#include <boost/numeric/itl/iteration/basic_iteration.hpp>
#include <boost/numeric/itl/iteration/cyclic_iteration.hpp>
#include <boost/numeric/itl/iteration/noisy_iteration.hpp>
#include <boost/numeric/itl/krylov/cg.hpp>
#include <boost/numeric/itl/krylov/cgs.hpp>
#include <boost/numeric/itl/krylov/bicg.hpp>
#include <boost/numeric/itl/krylov/bicgstab.hpp>
#include <boost/numeric/itl/krylov/bicgstab_2.hpp>
#include <boost/numeric/itl/krylov/bicgstab_ell.hpp>
#include <boost/numeric/itl/krylov/fsm.hpp>
#include <boost/numeric/itl/krylov/idr_s.hpp>
#include <boost/numeric/itl/krylov/gmres.hpp>
#include <boost/numeric/itl/krylov/tfqmr.hpp>
#include <boost/numeric/itl/krylov/qmr.hpp>
#include <boost/numeric/itl/krylov/repeating_solver.hpp>
#include <boost/numeric/itl/minimization/quasi_newton.hpp>
#include <boost/numeric/itl/pc/identity.hpp>
#include <boost/numeric/itl/pc/is_identity.hpp>
#include <boost/numeric/itl/pc/diagonal.hpp>
#include <boost/numeric/itl/pc/ilu.hpp>
#include <boost/numeric/itl/pc/ilu_0.hpp>
#include <boost/numeric/itl/pc/ilut.hpp>
#include <boost/numeric/itl/pc/ic_0.hpp>
#include <boost/numeric/itl/pc/imf_preconditioner.hpp>
#include <boost/numeric/itl/pc/imf_algorithms.hpp>
#include <boost/numeric/itl/pc/sub_matrix_pc.hpp>
#include <boost/numeric/itl/pc/concat.hpp>
#include <boost/numeric/itl/smoother/gauss_seidel.hpp>
#include <boost/numeric/itl/stepper/armijo.hpp>
#include <boost/numeric/itl/stepper/wolf.hpp>
#include <boost/numeric/itl/updater/bfgs.hpp>
#include <boost/numeric/itl/updater/broyden.hpp>
#include <boost/numeric/itl/updater/dfp.hpp>
#include <boost/numeric/itl/updater/psb.hpp>
#include <boost/numeric/itl/updater/sr1.hpp>
#endif // ITL_ITL_INCLUDE
// 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_ITL_FWD_INCLUDE
#define ITL_ITL_FWD_INCLUDE
#include <boost/numeric/mtl/concept/collection.hpp>
namespace itl {
template <class Real> class basic_iteration;
template <class Real, class OStream> class cyclic_iteration;
template <class Real, class OStream> class noisy_iteration;
template <typename Solver, typename VectorIn, bool trans> class solver_proxy;
namespace pc {
template <typename PC, typename Vector, bool> struct solver;
template <typename Matrix, typename Value> class identity;
// template <typename Matrix, typename Value, typename Vector> Vector solve(const identity<Matrix>&, const Vector& x);
// template <typename Matrix, typename Vector> Vector adjoint_solve(const identity<Matrix>&, const Vector& x);
template <typename Matrix, typename Value> class diagonal;
// template <typename Matrix, typename Vector> Vector solve(const diagonal<Matrix>& P, const Vector& x);
// template <typename Matrix, typename Vector> Vector adjoint_solve(const diagonal<Matrix>& P, const Vector& x);
template <typename Matrix, typename Factorizer, typename Value> class ilu;
template <typename Matrix, typename Value> class ilu_0; // Maybe we should declare the default here???
template <typename Matrix, typename Value> class ilut; // Maybe we should declare the default here???
// template <typename Matrix, typename Vector> Vector solve(const ilu_0<Matrix>& P, const Vector& x);
// template <typename Matrix, typename Vector> Vector adjoint_solve(const ilu_0<Matrix>& P, const Vector& x);
template <typename Matrix, typename Value> class ic_0; // Maybe we should declare the default here???
// template <typename Matrix, typename Vector> Vector solve(const ic_0<Matrix>& P, const Vector& x);
// template <typename Matrix, typename Value, typename Vector> Vector adjoint_solve(const ic_0<Matrix, Value>& P, const Vector& x);
} // namespace pc
template < typename LinearOperator, typename HilbertSpaceX, typename HilbertSpaceB,
typename Preconditioner, typename Iteration >
int cg(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
const Preconditioner& M, Iteration& iter);
template < typename LinearOperator, typename Preconditioner= pc::identity<LinearOperator, double>,
typename RightPreconditioner= pc::identity<LinearOperator, double> >
class cg_solver;
template < typename LinearOperator, typename Vector,
typename Preconditioner, typename Iteration >
int bicg(const LinearOperator &A, Vector &x, const Vector &b,
const Preconditioner &M, Iteration& iter);
template < class LinearOperator, class HilbertSpaceX, class HilbertSpaceB,
class Preconditioner, class Iteration >
int bicgstab(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
const Preconditioner& M, Iteration& iter);
template < class LinearOperator, class HilbertSpaceX, class HilbertSpaceB,
class Preconditioner, class Iteration >
int bicgstab_2(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
const Preconditioner& M, Iteration& iter);
template < typename LinearOperator, typename Vector,
typename LeftPreconditioner, typename RightPreconditioner,
typename Iteration >
int bicgstab_ell(const LinearOperator &A, Vector &x, const Vector &b,
const LeftPreconditioner &L, const RightPreconditioner &R,
Iteration& iter, size_t l);
template < typename LinearOperator, typename Preconditioner= pc::identity<LinearOperator, double>,
typename RightPreconditioner= pc::identity<LinearOperator, double> >
class bicgstab_ell_solver;
template <typename Solver, unsigned N, bool Stored= false>
class repeating_solver;
} // namespace itl
#endif // ITL_ITL_FWD_INCLUDE
// 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_BICG_INCLUDE
#define ITL_BICG_INCLUDE
#include <complex>
#include <boost/numeric/itl/itl_fwd.hpp>
#include <boost/numeric/mtl/concept/collection.hpp>
#include <boost/numeric/mtl/operation/conj.hpp>
#include <boost/numeric/mtl/operation/resource.hpp>
#include <boost/numeric/mtl/interface/vpt.hpp>
namespace itl {
/// Bi-Conjugate Gradient
template < typename LinearOperator, typename Vector,
typename Preconditioner, typename Iteration >
int bicg(const LinearOperator &A, Vector &x, const Vector &b,
const Preconditioner &M, Iteration& iter)
{
mtl::vampir_trace<7003> tracer;
using mtl::conj;
typedef typename mtl::Collection<Vector>::value_type Scalar;
Scalar rho_1(0), rho_2(0), alpha(0), beta(0);
Vector r(b - A * x), z(resource(x)), p(resource(x)), q(resource(x)),
r_tilde(r), z_tilde(resource(x)), p_tilde(resource(x)), q_tilde(resource(x));
while ( ! iter.finished(r)) {
++iter;
z= solve(M, r);
z_tilde= adjoint_solve(M, r_tilde);
rho_1= dot(z_tilde, z);
if (rho_1 == 0.) return iter.fail(2, "bicg breakdown");
if (iter.first()) {
p= z;
p_tilde= z_tilde;
} else {
beta= rho_1 / rho_2;
p= z + beta * p;
p_tilde= z_tilde + conj(beta) * p_tilde;
}
q= A * p;
q_tilde= adjoint(A) * p_tilde;
alpha= rho_1 / dot(p_tilde, q);
x+= alpha * p;
r-= alpha * q;
r_tilde-= conj(alpha) * q_tilde;
rho_2= rho_1;
}
return iter;
}
/// Solver class for BiCG method; right preconditioner ignored (prints warning if not identity)
template < typename LinearOperator, typename Preconditioner= pc::identity<LinearOperator>,
typename RightPreconditioner= pc::identity<LinearOperator> >
class bicg_solver
{
public:
/// Construct solver from a linear operator; generate (left) preconditioner from it
explicit bicg_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
bicg_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 bicg(A, x, b, L, iter);
}
/// Perform one bicg 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 bicg(A, x, b, L, iter);
}
private:
const LinearOperator& A;
Preconditioner L;
};
} // namespace itl
#endif // ITL_BICG_INCLUDE
// 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_BICGSTAB_INCLUDE
#define ITL_BICGSTAB_INCLUDE
#include <boost/numeric/mtl/concept/collection.hpp>
#include <boost/numeric/mtl/utility/exception.hpp>
#include <boost/numeric/itl/utility/exception.hpp>
#include <boost/numeric/mtl/operation/resource.hpp>
#include <boost/numeric/mtl/interface/vpt.hpp>
namespace itl {
/// Bi-Conjugate Gradient Stabilized
template < class LinearOperator, class HilbertSpaceX, class HilbertSpaceB,
class Preconditioner, class Iteration >
int bicgstab(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
const Preconditioner& M, Iteration& iter)
{
typedef typename mtl::Collection<HilbertSpaceX>::value_type Scalar;
typedef HilbertSpaceX Vector;