Loading...
Searching...
No Matches
solver.hh
Go to the documentation of this file.
// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_ISTL_SOLVER_HH
#define DUNE_ISTL_SOLVER_HH
#include <dune-istl-config.hh> // DUNE_ISTL_SUPPORT_OLD_CATEGORY_INTERFACE
#include <iomanip>
#include <ostream>
#include <string>
#include <functional>
#include <dune/common/exceptions.hh>
#include <dune/common/shared_ptr.hh>
#include <dune/common/simd/io.hh>
#include <dune/common/simd/simd.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/timer.hh>
#include "solvertype.hh"
#include "preconditioner.hh"
#include "operators.hh"
#include "scalarproducts.hh"
namespace Dune
{
struct InverseOperatorResult
{
InverseOperatorResult ()
{
clear();
}
void clear ()
{
iterations = 0;
reduction = 0;
converged = false;
conv_rate = 1;
elapsed = 0;
condition_estimate = -1;
}
int iterations;
double reduction;
bool converged;
double conv_rate;
double condition_estimate = -1;
double elapsed;
};
//=====================================================================
template<class X, class Y>
class InverseOperator {
public:
typedef X domain_type;
typedef Y range_type;
typedef typename X::field_type field_type;
typedef typename FieldTraits<field_type>::real_type real_type;
typedef Simd::Scalar<real_type> scalar_real_type;
virtual void apply (X& x, Y& b, InverseOperatorResult& res) = 0;
virtual void apply (X& x, Y& b, double reduction, InverseOperatorResult& res) = 0;
virtual SolverCategory::Category category() const
#ifdef DUNE_ISTL_SUPPORT_OLD_CATEGORY_INTERFACE
{
DUNE_THROW(Dune::Exception,"It is necessary to implement the category method in a derived classes, in the future this method will pure virtual.");
}
#else
= 0;
#endif
virtual ~InverseOperator () {}
protected:
// spacing values
enum { iterationSpacing = 5 , normSpacing = 16 };
void printHeader(std::ostream& s) const
{
s << std::setw(iterationSpacing) << " Iter";
s << std::setw(normSpacing) << "Defect";
s << std::setw(normSpacing) << "Rate" << std::endl;
}
template <typename CountType, typename DataType>
void printOutput(std::ostream& s,
const CountType& iter,
const DataType& norm,
const DataType& norm_old) const
{
const DataType rate = norm/norm_old;
s << std::setw(iterationSpacing) << iter << " ";
s << std::setw(normSpacing) << Simd::io(norm) << " ";
s << std::setw(normSpacing) << Simd::io(rate) << std::endl;
}
template <typename CountType, typename DataType>
void printOutput(std::ostream& s,
const CountType& iter,
const DataType& norm) const
{
s << std::setw(iterationSpacing) << iter << " ";
s << std::setw(normSpacing) << Simd::io(norm) << std::endl;
}
};
template<class X, class Y>
class IterativeSolver : public InverseOperator<X,Y>{
public:
using typename InverseOperator<X,Y>::domain_type;
using typename InverseOperator<X,Y>::range_type;
using typename InverseOperator<X,Y>::field_type;
using typename InverseOperator<X,Y>::real_type;
using typename InverseOperator<X,Y>::scalar_real_type;
IterativeSolver (const LinearOperator<X,Y>& op, Preconditioner<X,Y>& prec, scalar_real_type reduction, int maxit, int verbose) :
_op(stackobject_to_shared_ptr(op)),
_prec(stackobject_to_shared_ptr(prec)),
_sp(new SeqScalarProduct<X>),
_reduction(reduction), _maxit(maxit), _verbose(verbose), _category(SolverCategory::sequential)
{
if(SolverCategory::category(op) != SolverCategory::sequential)
DUNE_THROW(InvalidSolverCategory, "LinearOperator has to be sequential!");
if(SolverCategory::category(prec) != SolverCategory::sequential)
DUNE_THROW(InvalidSolverCategory, "Preconditioner has to be sequential!");
}
IterativeSolver (const LinearOperator<X,Y>& op, const ScalarProduct<X>& sp, Preconditioner<X,Y>& prec,
scalar_real_type reduction, int maxit, int verbose) :
_op(stackobject_to_shared_ptr(op)),
_prec(stackobject_to_shared_ptr(prec)),
_sp(stackobject_to_shared_ptr(sp)),
_reduction(reduction), _maxit(maxit), _verbose(verbose), _category(SolverCategory::category(op))
{
if(SolverCategory::category(op) != SolverCategory::category(prec))
DUNE_THROW(InvalidSolverCategory, "LinearOperator and Preconditioner must have the same SolverCategory!");
if(SolverCategory::category(op) != SolverCategory::category(sp))
DUNE_THROW(InvalidSolverCategory, "LinearOperator and ScalarProduct must have the same SolverCategory!");
}
IterativeSolver (std::shared_ptr<const LinearOperator<X,Y> > op, std::shared_ptr<Preconditioner<X,X> > prec, const ParameterTree& configuration) :
IterativeSolver(op,std::make_shared<SeqScalarProduct<X>>(),prec,
configuration.get<real_type>("reduction"),
configuration.get<int>("maxit"),
configuration.get<int>("verbose"))
{}
IterativeSolver (std::shared_ptr<const LinearOperator<X,Y> > op, std::shared_ptr<const ScalarProduct<X> > sp, std::shared_ptr<Preconditioner<X,X> > prec, const ParameterTree& configuration) :
IterativeSolver(op,sp,prec,
configuration.get<scalar_real_type>("reduction"),
configuration.get<int>("maxit"),
configuration.get<int>("verbose"))
{}
IterativeSolver (std::shared_ptr<const LinearOperator<X,Y>> op,
std::shared_ptr<const ScalarProduct<X>> sp,
std::shared_ptr<Preconditioner<X,Y>> prec,
scalar_real_type reduction, int maxit, int verbose) :
_op(op),
_prec(prec),
_sp(sp),
_reduction(reduction), _maxit(maxit), _verbose(verbose),
_category(SolverCategory::category(*op))
{
if(SolverCategory::category(*op) != SolverCategory::category(*prec))
DUNE_THROW(InvalidSolverCategory, "LinearOperator and Preconditioner must have the same SolverCategory!");
if(SolverCategory::category(*op) != SolverCategory::category(*sp))
DUNE_THROW(InvalidSolverCategory, "LinearOperator and ScalarProduct must have the same SolverCategory!");
}
// #warning actually we want to have this as the default and just implement the second one
// //! \copydoc InverseOperator::apply(X&,Y&,InverseOperatorResult&)
// virtual void apply (X& x, Y& b, InverseOperatorResult& res)
// {
// apply(x,b,_reduction,res);
// }
#ifndef DOXYGEN
// make sure the three-argument apply from the base class does not get shadowed
// by the redefined four-argument version below
using InverseOperator<X,Y>::apply;
#endif
void apply (X& x, X& b, double reduction, InverseOperatorResult& res) override
{
scalar_real_type saved_reduction = _reduction;
_reduction = reduction;
this->apply(x,b,res);
_reduction = saved_reduction;
}
SolverCategory::Category category() const override
{
return _category;
}
std::string name() const{
std::string name = className(*this);
return name.substr(0, name.find("<"));
}
template<class CountType = unsigned int>
class Iteration {
public:
Iteration(const IterativeSolver& parent, InverseOperatorResult& res)
: _i(0)
, _res(res)
, _parent(parent)
, _valid(true)
{
res.clear();
if(_parent._verbose>0){
std::cout << "=== " << parent.name() << std::endl;
if(_parent._verbose > 1)
_parent.printHeader(std::cout);
}
}
Iteration(const Iteration&) = delete;
Iteration(Iteration&& other)
: _def0(other._def0)
, _def(other._def)
, _i(other._i)
, _watch(other._watch)
, _res(other._res)
, _parent(other._parent)
, _valid(other._valid)
{
other._valid = false;
}
~Iteration(){
if(_valid)
finalize();
}
bool step(CountType i, real_type def){
if (!Simd::allTrue(isFinite(def))) // check for inf or NaN
{
if (_parent._verbose>0)
std::cout << "=== " << _parent.name() << ": abort due to infinite or NaN defect"
<< std::endl;
DUNE_THROW(SolverAbort,
_parent.name() << ": defect=" << Simd::io(def)
<< " is infinite or NaN");
}
if(i == 0)
_def0 = def;
if(_parent._verbose > 1){
if(i!=0)
_parent.printOutput(std::cout,i,def,_def);
else
_parent.printOutput(std::cout,i,def);
}
_def = def;
_i = i;
_res.converged = (Simd::allTrue(def<_def0*_parent._reduction || def<real_type(1E-30))); // convergence check
return _res.converged;
}
protected:
void finalize(){
_res.converged = (Simd::allTrue(_def<_def0*_parent._reduction || _def<real_type(1E-30)));
_res.iterations = _i;
_res.reduction = static_cast<double>(Simd::max(_def/_def0));
_res.conv_rate = pow(_res.reduction,1.0/_i);
_res.elapsed = _watch.elapsed();
if (_parent._verbose>0) // final print
{
std::cout << "=== rate=" << _res.conv_rate
<< ", T=" << _res.elapsed
<< ", TIT=" << _res.elapsed/_res.iterations
<< ", IT=" << _res.iterations << std::endl;
}
}
real_type _def0 = 0.0, _def = 0.0;
CountType _i;
Timer _watch;
InverseOperatorResult& _res;
const IterativeSolver& _parent;
bool _valid;
};
protected:
std::shared_ptr<const LinearOperator<X,Y>> _op;
std::shared_ptr<Preconditioner<X,Y>> _prec;
std::shared_ptr<const ScalarProduct<X>> _sp;
scalar_real_type _reduction;
int _maxit;
int _verbose;
SolverCategory::Category _category;
};
template <typename ISTLLinearSolver, typename BCRSMatrix>
class SolverHelper
{
public:
static void setMatrix (ISTLLinearSolver& solver,
const BCRSMatrix& matrix)
{
static const bool is_direct_solver
= Dune::IsDirectSolver<ISTLLinearSolver>::value;
SolverHelper<ISTLLinearSolver,BCRSMatrix>::
Implementation<is_direct_solver>::setMatrix(solver,matrix);
}
protected:
template <bool is_direct_solver, typename Dummy = void>
struct Implementation
{
static void setMatrix (ISTLLinearSolver&,
const BCRSMatrix&)
{}
};
template <typename Dummy>
struct Implementation<true,Dummy>
{
static void setMatrix (ISTLLinearSolver& solver,
const BCRSMatrix& matrix)
{
solver.setMatrix(matrix);
}
};
};
}
#endif
operators.hh
Define general, extensible interface for operators. The available implementation wraps a matrix.
solvertype.hh
Templates characterizing the type of a solver.
preconditioner.hh
scalarproducts.hh
Define base class for scalar product and norm.
std
STL namespace.
Dune
Definition allocator.hh:11
Dune::get
PropertyMapTypeSelector< Amg::VertexVisitedTag, Amg::PropertiesGraph< G, Amg::VertexProperties, EP, VM, EM > >::Type get(const Amg::VertexVisitedTag &tag, Amg::PropertiesGraph< G, Amg::VertexProperties, EP, VM, EM > &graph)
Definition dependency.hh:293
Dune::BCRSMatrix
A sparse block matrix with compressed row storage.
Definition bcrsmatrix.hh:467
Dune::SolverAbort
Thrown when a solver aborts due to some problem.
Definition istlexception.hh:46
Dune::LinearOperator
A linear operator.
Definition operators.hh:68
Dune::Preconditioner
Base class for matrix free definition of preconditioners.
Definition preconditioner.hh:33
Dune::ScalarProduct
Base class for scalar product and norm computation.
Definition scalarproducts.hh:52
Dune::SeqScalarProduct
Default implementation for the scalar case.
Definition scalarproducts.hh:168
Dune::InverseOperatorResult
Statistics about the application of an inverse operator.
Definition solver.hh:50
Dune::InverseOperatorResult::InverseOperatorResult
InverseOperatorResult()
Default constructor.
Definition solver.hh:52
Dune::InverseOperatorResult::condition_estimate
double condition_estimate
Estimate of condition number.
Definition solver.hh:81
Dune::InverseOperatorResult::elapsed
double elapsed
Elapsed time in seconds.
Definition solver.hh:84
Dune::InverseOperatorResult::iterations
int iterations
Number of iterations.
Definition solver.hh:69
Dune::InverseOperatorResult::reduction
double reduction
Reduction achieved: .
Definition solver.hh:72
Dune::InverseOperatorResult::clear
void clear()
Resets all data.
Definition solver.hh:58
Dune::InverseOperatorResult::conv_rate
double conv_rate
Convergence rate (average reduction per step).
Definition solver.hh:78
Dune::InverseOperatorResult::converged
bool converged
True if convergence criterion has been met.
Definition solver.hh:75
Dune::InverseOperator
Abstract base class for all solvers.
Definition solver.hh:101
Dune::InverseOperator::printHeader
void printHeader(std::ostream &s) const
helper function for printing header of solver output
Definition solver.hh:165
Dune::InverseOperator::~InverseOperator
virtual ~InverseOperator()
Destructor.
Definition solver.hh:158
Dune::InverseOperator::printOutput
void printOutput(std::ostream &s, const CountType &iter, const DataType &norm) const
helper function for printing solver output
Definition solver.hh:187
Dune::InverseOperator::printOutput
void printOutput(std::ostream &s, const CountType &iter, const DataType &norm, const DataType &norm_old) const
helper function for printing solver output
Definition solver.hh:174
Dune::InverseOperator::apply
virtual void apply(X &x, Y &b, double reduction, InverseOperatorResult &res)=0
apply inverse operator, with given convergence criteria.
Dune::InverseOperator< X, X >::scalar_real_type
Simd::Scalar< real_type > scalar_real_type
Definition solver.hh:116
Dune::InverseOperator< X, X >::range_type
X range_type
Definition solver.hh:107
Dune::InverseOperator< X, X >::normSpacing
@ normSpacing
Definition solver.hh:162
Dune::InverseOperator< X, X >::iterationSpacing
@ iterationSpacing
Definition solver.hh:162
Dune::InverseOperator< X, X >::domain_type
X domain_type
Definition solver.hh:104
Dune::InverseOperator::apply
virtual void apply(X &x, Y &b, InverseOperatorResult &res)=0
Apply inverse operator,.
Dune::InverseOperator< X, X >::field_type
X::field_type field_type
Definition solver.hh:110
Dune::InverseOperator< X, X >::real_type
FieldTraits< field_type >::real_type real_type
Definition solver.hh:113
Dune::InverseOperator::category
virtual SolverCategory::Category category() const =0
Category of the solver (see SolverCategory::Category).
Dune::IterativeSolver::IterativeSolver
IterativeSolver(std::shared_ptr< const LinearOperator< X, Y > > op, std::shared_ptr< const ScalarProduct< X > > sp, std::shared_ptr< Preconditioner< X, X > > prec, const ParameterTree &configuration)
Constructor.
Definition solver.hh:315
Dune::IterativeSolver::category
SolverCategory::Category category() const override
Category of the solver (see SolverCategory::Category).
Definition solver.hh:385
Dune::IterativeSolver::IterativeSolver
IterativeSolver(std::shared_ptr< const LinearOperator< X, Y > > op, std::shared_ptr< Preconditioner< X, X > > prec, const ParameterTree &configuration)
Constructor.
Definition solver.hh:292
Dune::IterativeSolver::_sp
std::shared_ptr< const ScalarProduct< X > > _sp
Definition solver.hh:508
Dune::IterativeSolver::apply
void apply(X &x, X &b, double reduction, InverseOperatorResult &res) override
Apply inverse operator with given reduction factor.
Definition solver.hh:376
Dune::IterativeSolver::IterativeSolver
IterativeSolver(std::shared_ptr< const LinearOperator< X, Y > > op, std::shared_ptr< const ScalarProduct< X > > sp, std::shared_ptr< Preconditioner< X, Y > > prec, scalar_real_type reduction, int maxit, int verbose)
General constructor to initialize an iterative solver.
Definition solver.hh:342
Dune::IterativeSolver::name
std::string name() const
Definition solver.hh:390
Dune::IterativeSolver::IterativeSolver
IterativeSolver(const LinearOperator< X, Y > &op, Preconditioner< X, Y > &prec, scalar_real_type reduction, int maxit, int verbose)
General constructor to initialize an iterative solver.
Definition solver.hh:232
Dune::IterativeSolver::_op
std::shared_ptr< const LinearOperator< X, Y > > _op
Definition solver.hh:506
Dune::IterativeSolver::_maxit
int _maxit
Definition solver.hh:510
Dune::IterativeSolver::_verbose
int _verbose
Definition solver.hh:511
Dune::IterativeSolver::_reduction
scalar_real_type _reduction
Definition solver.hh:509
Dune::IterativeSolver::IterativeSolver
IterativeSolver(const LinearOperator< X, Y > &op, const ScalarProduct< X > &sp, Preconditioner< X, Y > &prec, scalar_real_type reduction, int maxit, int verbose)
General constructor to initialize an iterative solver.
Definition solver.hh:264
Dune::IterativeSolver::_category
SolverCategory::Category _category
Definition solver.hh:512
Dune::IterativeSolver::_prec
std::shared_ptr< Preconditioner< X, Y > > _prec
Definition solver.hh:507
Dune::IterativeSolver::Iteration::Iteration
Iteration(const IterativeSolver &parent, InverseOperatorResult &res)
Definition solver.hh:415
Dune::IterativeSolver::Iteration::Iteration
Iteration(Iteration &&other)
Definition solver.hh:430
Dune::IterativeSolver::Iteration::_res
InverseOperatorResult & _res
Definition solver.hh:500
Dune::IterativeSolver::Iteration::_parent
const IterativeSolver & _parent
Definition solver.hh:501
Dune::IterativeSolver::Iteration::_watch
Timer _watch
Definition solver.hh:499
Dune::IterativeSolver::Iteration::Iteration
Iteration(const Iteration &)=delete
Dune::IterativeSolver::Iteration::_i
CountType _i
Definition solver.hh:498
Dune::IterativeSolver::Iteration::_def0
real_type _def0
Definition solver.hh:497
Dune::IterativeSolver::Iteration::step
bool step(CountType i, real_type def)
registers the iteration step, checks for invalid defect norm and convergence.
Definition solver.hh:457
Dune::IterativeSolver::Iteration::finalize
void finalize()
Definition solver.hh:482
Dune::IterativeSolver::Iteration::_def
real_type _def
Definition solver.hh:497
Dune::IterativeSolver::Iteration::~Iteration
~Iteration()
Definition solver.hh:442
Dune::IterativeSolver::Iteration::_valid
bool _valid
Definition solver.hh:502
Dune::SolverHelper
Helper class for notifying a DUNE-ISTL linear solver about a change of the iteration matrix object in...
Definition solver.hh:524
Dune::SolverHelper::setMatrix
static void setMatrix(ISTLLinearSolver &solver, const BCRSMatrix &matrix)
Definition solver.hh:526
Dune::SolverHelper::Implementation
Implementation that works together with iterative ISTL solvers, e.g. Dune::CGSolver or Dune::BiCGSTAB...
Definition solver.hh:542
Dune::SolverHelper::Implementation::setMatrix
static void setMatrix(ISTLLinearSolver &, const BCRSMatrix &)
Definition solver.hh:543
Dune::SolverHelper::Implementation< true, Dummy >::setMatrix
static void setMatrix(ISTLLinearSolver &solver, const BCRSMatrix &matrix)
Definition solver.hh:555
Dune::SolverCategory
Categories for the solvers.
Definition solvercategory.hh:22
Dune::SolverCategory::Category
Category
Definition solvercategory.hh:23
Dune::SolverCategory::sequential
@ sequential
Category for sequential solvers.
Definition solvercategory.hh:25
Dune::SolverCategory::category
static Category category(const OP &op, decltype(op.category()) *=nullptr)
Helperfunction to extract the solver category either from an enum, or from the newly introduced virtu...
Definition solvercategory.hh:34
Dune::InvalidSolverCategory
Definition solvercategory.hh:54
Dune::IsDirectSolver::value
@ value
Whether this is a direct solver.
Definition solvertype.hh:24