dune-istl  2.7.0
solver.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 
4 #ifndef DUNE_ISTL_SOLVER_HH
5 #define DUNE_ISTL_SOLVER_HH
6 
7 #include <iomanip>
8 #include <ostream>
9 #include <string>
10 #include <functional>
11 
12 #include <dune/common/exceptions.hh>
13 #include <dune/common/shared_ptr.hh>
14 #include <dune/common/simd/io.hh>
15 #include <dune/common/simd/simd.hh>
16 #include <dune/common/parametertree.hh>
17 #include <dune/common/timer.hh>
18 
19 #include "solvertype.hh"
20 #include "preconditioner.hh"
21 #include "operators.hh"
22 #include "scalarproducts.hh"
23 
24 namespace Dune
25 {
46  {
49  {
50  clear();
51  }
52 
54  void clear ()
55  {
56  iterations = 0;
57  reduction = 0;
58  converged = false;
59  conv_rate = 1;
60  elapsed = 0;
61  condition_estimate = -1;
62  }
63 
66 
68  double reduction;
69 
71  bool converged;
72 
74  double conv_rate;
75 
77  double condition_estimate = -1;
78 
80  double elapsed;
81  };
82 
83 
84  //=====================================================================
96  template<class X, class Y>
98  public:
100  typedef X domain_type;
101 
103  typedef Y range_type;
104 
106  typedef typename X::field_type field_type;
107 
109  typedef typename FieldTraits<field_type>::real_type real_type;
110 
112  typedef Simd::Scalar<real_type> scalar_real_type;
113 
126  virtual void apply (X& x, Y& b, InverseOperatorResult& res) = 0;
127 
141  virtual void apply (X& x, Y& b, double reduction, InverseOperatorResult& res) = 0;
142 
144  virtual SolverCategory::Category category() const
145 #ifdef DUNE_ISTL_SUPPORT_OLD_CATEGORY_INTERFACE
146  {
147  DUNE_THROW(Dune::Exception,"It is necessary to implement the category method in a derived classes, in the future this method will pure virtual.");
148  }
149 #else
150  = 0;
151 #endif
152 
154  virtual ~InverseOperator () {}
155 
156  protected:
157  // spacing values
158  enum { iterationSpacing = 5 , normSpacing = 16 };
159 
161  void printHeader(std::ostream& s) const
162  {
163  s << std::setw(iterationSpacing) << " Iter";
164  s << std::setw(normSpacing) << "Defect";
165  s << std::setw(normSpacing) << "Rate" << std::endl;
166  }
167 
169  template <typename CountType, typename DataType>
170  void printOutput(std::ostream& s,
171  const CountType& iter,
172  const DataType& norm,
173  const DataType& norm_old) const
174  {
175  const DataType rate = norm/norm_old;
176  s << std::setw(iterationSpacing) << iter << " ";
177  s << std::setw(normSpacing) << Simd::io(norm) << " ";
178  s << std::setw(normSpacing) << Simd::io(rate) << std::endl;
179  }
180 
182  template <typename CountType, typename DataType>
183  void printOutput(std::ostream& s,
184  const CountType& iter,
185  const DataType& norm) const
186  {
187  s << std::setw(iterationSpacing) << iter << " ";
188  s << std::setw(normSpacing) << Simd::io(norm) << std::endl;
189  }
190  };
191 
200  template<class X, class Y>
201  class IterativeSolver : public InverseOperator<X,Y>{
202  public:
203  using typename InverseOperator<X,Y>::domain_type;
204  using typename InverseOperator<X,Y>::range_type;
205  using typename InverseOperator<X,Y>::field_type;
206  using typename InverseOperator<X,Y>::real_type;
208 
228  IterativeSolver (LinearOperator<X,Y>& op, Preconditioner<X,Y>& prec, scalar_real_type reduction, int maxit, int verbose) :
229  _op(stackobject_to_shared_ptr(op)),
230  _prec(stackobject_to_shared_ptr(prec)),
231  _sp(new SeqScalarProduct<X>),
232  _reduction(reduction), _maxit(maxit), _verbose(verbose), _category(SolverCategory::sequential)
233  {
235  DUNE_THROW(InvalidSolverCategory, "LinearOperator has to be sequential!");
237  DUNE_THROW(InvalidSolverCategory, "Preconditioner has to be sequential!");
238  }
239 
261  scalar_real_type reduction, int maxit, int verbose) :
262  _op(stackobject_to_shared_ptr(op)),
263  _prec(stackobject_to_shared_ptr(prec)),
264  _sp(stackobject_to_shared_ptr(sp)),
265  _reduction(reduction), _maxit(maxit), _verbose(verbose), _category(SolverCategory::category(op))
266  {
268  DUNE_THROW(InvalidSolverCategory, "LinearOperator and Preconditioner must have the same SolverCategory!");
270  DUNE_THROW(InvalidSolverCategory, "LinearOperator and ScalarProduct must have the same SolverCategory!");
271  }
272 
288  IterativeSolver (std::shared_ptr<LinearOperator<X,Y> > op, std::shared_ptr<Preconditioner<X,X> > prec, const ParameterTree& configuration) :
289  IterativeSolver(op,std::make_shared<SeqScalarProduct<X>>(),prec,
290  configuration.get<real_type>("reduction"),
291  configuration.get<int>("maxit"),
292  configuration.get<int>("verbose"))
293  {}
294 
311  IterativeSolver (std::shared_ptr<LinearOperator<X,Y> > op, std::shared_ptr<ScalarProduct<X> > sp, std::shared_ptr<Preconditioner<X,X> > prec, const ParameterTree& configuration) :
312  IterativeSolver(op,sp,prec,
313  configuration.get<scalar_real_type>("reduction"),
314  configuration.get<int>("maxit"),
315  configuration.get<int>("verbose"))
316  {}
317 
338  IterativeSolver (std::shared_ptr<LinearOperator<X,Y>> op,
339  std::shared_ptr<ScalarProduct<X>> sp,
340  std::shared_ptr<Preconditioner<X,Y>> prec,
341  scalar_real_type reduction, int maxit, int verbose) :
342  _op(op),
343  _prec(prec),
344  _sp(sp),
345  _reduction(reduction), _maxit(maxit), _verbose(verbose),
347  {
349  DUNE_THROW(InvalidSolverCategory, "LinearOperator and Preconditioner must have the same SolverCategory!");
351  DUNE_THROW(InvalidSolverCategory, "LinearOperator and ScalarProduct must have the same SolverCategory!");
352  }
353 
354  // #warning actually we want to have this as the default and just implement the second one
355  // //! \copydoc InverseOperator::apply(X&,Y&,InverseOperatorResult&)
356  // virtual void apply (X& x, Y& b, InverseOperatorResult& res)
357  // {
358  // apply(x,b,_reduction,res);
359  // }
360 
361 #ifndef DOXYGEN
362  // make sure the three-argument apply from the base class does not get shadowed
363  // by the redefined four-argument version below
365 #endif
366 
372  virtual void apply (X& x, X& b, double reduction, InverseOperatorResult& res)
373  {
374  scalar_real_type saved_reduction = _reduction;
375  _reduction = reduction;
376  this->apply(x,b,res);
377  _reduction = saved_reduction;
378  }
379 
382  {
383  return _category;
384  }
385 
386  std::string name() const{
387  std::string name = className(*this);
388  return name.substr(0, name.find("<"));
389  }
390 
408  template<class CountType = unsigned int>
409  class Iteration {
410  public:
412  : _i(0)
413  , _res(res)
414  , _parent(parent)
415  , _valid(true)
416  {
417  res.clear();
418  if(_parent._verbose>0){
419  std::cout << "=== " << parent.name() << std::endl;
420  if(_parent._verbose > 1)
421  _parent.printHeader(std::cout);
422  }
423  }
424 
425  Iteration(const Iteration&) = delete;
427  : _def0(other._def0)
428  , _def(other._def)
429  , _i(other._i)
430  , _watch(other._watch)
431  , _res(other._res)
432  , _parent(other._parent)
433  , _valid(other._valid)
434  {
435  other._valid = false;
436  }
437 
439  if(_valid)
440  finalize();
441  }
442 
453  bool step(CountType i, real_type def){
454  if (!Simd::allTrue(isFinite(def))) // check for inf or NaN
455  {
456  if (_parent._verbose>0)
457  std::cout << "=== " << _parent.name() << ": abort due to infinite or NaN defect"
458  << std::endl;
459  DUNE_THROW(SolverAbort,
460  _parent.name() << ": defect=" << Simd::io(def)
461  << " is infinite or NaN");
462  }
463  if(i == 0)
464  _def0 = def;
465  if(_parent._verbose > 1){
466  if(i!=0)
467  _parent.printOutput(std::cout,i,def,_def);
468  else
469  _parent.printOutput(std::cout,i,def);
470  }
471  _def = def;
472  _i = i;
473  _res.converged = (Simd::allTrue(def<_def0*_parent._reduction) || Simd::max(def)<1E-30); // convergence check
474  return _res.converged;
475  }
476 
477  protected:
478  void finalize(){
479  _res.converged = (Simd::allTrue(_def<_def0*_parent._reduction) || Simd::max(_def)<1E-30);
480  _res.iterations = _i;
481  _res.reduction = static_cast<double>(Simd::max(_def/_def0));
482  _res.conv_rate = pow(_res.reduction,1.0/_i);
483  _res.elapsed = _watch.elapsed();
484  if (_parent._verbose>0) // final print
485  {
486  std::cout << "=== rate=" << _res.conv_rate
487  << ", T=" << _res.elapsed
488  << ", TIT=" << _res.elapsed/_res.iterations
489  << ", IT=" << _res.iterations << std::endl;
490  }
491  }
492 
493  real_type _def0 = 0.0, _def = 0.0;
494  CountType _i;
495  Timer _watch;
498  bool _valid;
499  };
500 
501  protected:
502  std::shared_ptr<LinearOperator<X,Y>> _op;
503  std::shared_ptr<Preconditioner<X,Y>> _prec;
504  std::shared_ptr<ScalarProduct<X>> _sp;
506  int _maxit;
507  int _verbose;
509  };
510 
518  template <typename ISTLLinearSolver, typename BCRSMatrix>
520  {
521  public:
522  static void setMatrix (ISTLLinearSolver& solver,
523  const BCRSMatrix& matrix)
524  {
525  static const bool is_direct_solver
529  }
530 
531  protected:
536  template <bool is_direct_solver, typename Dummy = void>
538  {
539  static void setMatrix (ISTLLinearSolver&,
540  const BCRSMatrix&)
541  {}
542  };
543 
548  template <typename Dummy>
549  struct Implementation<true,Dummy>
550  {
551  static void setMatrix (ISTLLinearSolver& solver,
552  const BCRSMatrix& matrix)
553  {
554  solver.setMatrix(matrix);
555  }
556  };
557  };
558 
562 }
563 
564 #endif
Dune::IterativeSolver::_op
std::shared_ptr< LinearOperator< X, Y > > _op
Definition: solver.hh:502
Dune::IterativeSolver::Iteration::_res
InverseOperatorResult & _res
Definition: solver.hh:496
Dune::IterativeSolver::category
virtual SolverCategory::Category category() const
Category of the solver (see SolverCategory::Category)
Definition: solver.hh:381
Dune::IterativeSolver::_maxit
int _maxit
Definition: solver.hh:506
Dune::InverseOperator::domain_type
X domain_type
Type of the domain of the operator to be inverted.
Definition: solver.hh:100
Dune::SolverCategory::sequential
@ sequential
Category for sequential solvers.
Definition: solvercategory.hh:23
Dune::IterativeSolver::Iteration::_def
real_type _def
Definition: solver.hh:493
scalarproducts.hh
Define base class for scalar product and norm.
Dune::InverseOperatorResult::elapsed
double elapsed
Elapsed time in seconds.
Definition: solver.hh:80
Dune::SolverAbort
Thrown when a solver aborts due to some problem.
Definition: istlexception.hh:43
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:170
Dune::SolverHelper::Implementation
Implementation that works together with iterative ISTL solvers, e.g. Dune::CGSolver or Dune::BiCGSTAB...
Definition: solver.hh:537
Dune::InverseOperator
Abstract base class for all solvers.
Definition: solver.hh:97
Dune::IterativeSolver::Iteration::_i
CountType _i
Definition: solver.hh:494
Dune::InverseOperatorResult::condition_estimate
double condition_estimate
Estimate of condition number.
Definition: solver.hh:77
Dune::IterativeSolver::IterativeSolver
IterativeSolver(std::shared_ptr< LinearOperator< X, Y > > op, std::shared_ptr< ScalarProduct< X > > sp, std::shared_ptr< Preconditioner< X, X > > prec, const ParameterTree &configuration)
Constructor.
Definition: solver.hh:311
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:453
Dune::InverseOperatorResult::InverseOperatorResult
InverseOperatorResult()
Default constructor.
Definition: solver.hh:48
Dune::IterativeSolver::_prec
std::shared_ptr< Preconditioner< X, Y > > _prec
Definition: solver.hh:503
Dune::SolverCategory
Categories for the solvers.
Definition: solvercategory.hh:19
Dune::SolverHelper
Helper class for notifying a DUNE-ISTL linear solver about a change of the iteration matrix object in...
Definition: solver.hh:519
Dune::InverseOperatorResult::clear
void clear()
Resets all data.
Definition: solver.hh:54
Dune::SeqScalarProduct
Default implementation for the scalar case.
Definition: scalarproducts.hh:152
Dune::InverseOperator::range_type
Y range_type
Type of the range of the operator to be inverted.
Definition: solver.hh:103
Dune::IterativeSolver::_category
SolverCategory::Category _category
Definition: solver.hh:508
Dune::InverseOperator::scalar_real_type
Simd::Scalar< real_type > scalar_real_type
scalar type underlying the field_type
Definition: solver.hh:112
Dune::ScalarProduct
Base class for scalar product and norm computation.
Definition: scalarproducts.hh:48
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:292
Dune::IterativeSolver::_reduction
scalar_real_type _reduction
Definition: solver.hh:505
Dune::IterativeSolver::apply
virtual void apply(X &x, X &b, double reduction, InverseOperatorResult &res)
Apply inverse operator with given reduction factor.
Definition: solver.hh:372
Dune::IterativeSolver::Iteration::Iteration
Iteration(Iteration &&other)
Definition: solver.hh:426
Dune::InverseOperator::normSpacing
@ normSpacing
Definition: solver.hh:158
Dune::IterativeSolver::Iteration::finalize
void finalize()
Definition: solver.hh:478
Dune::LinearOperator
A linear operator.
Definition: operators.hh:65
Dune::SolverHelper::setMatrix
static void setMatrix(ISTLLinearSolver &solver, const BCRSMatrix &matrix)
Definition: solver.hh:522
Dune::IterativeSolver::Iteration::_def0
real_type _def0
Definition: solver.hh:493
Dune::InverseOperatorResult::converged
bool converged
True if convergence criterion has been met.
Definition: solver.hh:71
Dune::IterativeSolver::IterativeSolver
IterativeSolver(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:228
Dune::InverseOperatorResult
Statistics about the application of an inverse operator.
Definition: solver.hh:45
Dune::InverseOperator::printHeader
void printHeader(std::ostream &s) const
helper function for printing header of solver output
Definition: solver.hh:161
Dune::InverseOperatorResult::iterations
int iterations
Number of iterations.
Definition: solver.hh:65
Dune::BCRSMatrix
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:424
Dune::IterativeSolver::IterativeSolver
IterativeSolver(std::shared_ptr< LinearOperator< X, Y > > op, std::shared_ptr< Preconditioner< X, X > > prec, const ParameterTree &configuration)
Constructor.
Definition: solver.hh:288
Dune::IterativeSolver::IterativeSolver
IterativeSolver(std::shared_ptr< LinearOperator< X, Y >> op, std::shared_ptr< 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:338
Dune::InverseOperator::~InverseOperator
virtual ~InverseOperator()
Destructor.
Definition: solver.hh:154
Dune::InverseOperatorResult::conv_rate
double conv_rate
Convergence rate (average reduction per step)
Definition: solver.hh:74
Dune
Definition: allocator.hh:7
Dune::IterativeSolver
Base class for all implementations of iterative solvers.
Definition: solver.hh:201
Dune::SolverHelper::Implementation< true, Dummy >::setMatrix
static void setMatrix(ISTLLinearSolver &solver, const BCRSMatrix &matrix)
Definition: solver.hh:551
Dune::IsDirectSolver
Definition: solvertype.hh:13
Dune::IterativeSolver::Iteration::~Iteration
~Iteration()
Definition: solver.hh:438
preconditioner.hh
Dune::InverseOperator::real_type
FieldTraits< field_type >::real_type real_type
The real type of the field type (is the same if using real numbers, but differs for std::complex)
Definition: solver.hh:109
Dune::IterativeSolver::_sp
std::shared_ptr< ScalarProduct< X > > _sp
Definition: solver.hh:504
Dune::SolverHelper::Implementation::setMatrix
static void setMatrix(ISTLLinearSolver &, const BCRSMatrix &)
Definition: solver.hh:539
Dune::InverseOperator::field_type
X::field_type field_type
The field type of the operator.
Definition: solver.hh:106
Dune::IterativeSolver::Iteration::_valid
bool _valid
Definition: solver.hh:498
Dune::IterativeSolver::Iteration::Iteration
Iteration(const IterativeSolver &parent, InverseOperatorResult &res)
Definition: solver.hh:411
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:32
Dune::InverseOperator::iterationSpacing
@ iterationSpacing
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:183
Dune::IterativeSolver::name
std::string name() const
Definition: solver.hh:386
Dune::InverseOperator::apply
virtual void apply(X &x, Y &b, InverseOperatorResult &res)=0
Apply inverse operator,.
Dune::InverseOperator::category
virtual SolverCategory::Category category() const =0
Category of the solver (see SolverCategory::Category)
Dune::Preconditioner
Base class for matrix free definition of preconditioners.
Definition: preconditioner.hh:30
Dune::InvalidSolverCategory
Definition: solvercategory.hh:52
solvertype.hh
Templates characterizing the type of a solver.
Dune::IterativeSolver::IterativeSolver
IterativeSolver(LinearOperator< X, Y > &op, 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:260
Dune::IterativeSolver::Iteration::_watch
Timer _watch
Definition: solver.hh:495
Dune::SolverCategory::Category
Category
Definition: solvercategory.hh:21
Dune::IterativeSolver::Iteration::_parent
const IterativeSolver & _parent
Definition: solver.hh:497
Dune::IterativeSolver::Iteration
Class for controlling iterative methods.
Definition: solver.hh:409
Dune::IterativeSolver::_verbose
int _verbose
Definition: solver.hh:507
Dune::InverseOperatorResult::reduction
double reduction
Reduction achieved: .
Definition: solver.hh:68
operators.hh
Define general, extensible interface for operators. The available implementation wraps a matrix.