dune-istl  2.6-git
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 
10 #include <dune/common/exceptions.hh>
11 #include <dune/common/shared_ptr.hh>
12 #include <dune/common/simd.hh>
13 
14 #include "solvertype.hh"
15 #include "preconditioner.hh"
16 #include "operators.hh"
17 #include "scalarproducts.hh"
18 
19 namespace Dune
20 {
41  {
44  {
45  clear();
46  }
47 
49  void clear ()
50  {
51  iterations = 0;
52  reduction = 0;
53  converged = false;
54  conv_rate = 1;
55  elapsed = 0;
56  }
57 
60 
62  double reduction;
63 
65  bool converged;
66 
68  double conv_rate;
69 
71  double condition_estimate = -1;
72 
74  double elapsed;
75  };
76 
77 
78  //=====================================================================
90  template<class X, class Y>
92  public:
94  typedef X domain_type;
95 
97  typedef Y range_type;
98 
100  typedef typename X::field_type field_type;
101 
103  typedef typename FieldTraits<field_type>::real_type real_type;
104 
106  typedef SimdScalar<real_type> scalar_real_type;
107 
120  virtual void apply (X& x, Y& b, InverseOperatorResult& res) = 0;
121 
135  virtual void apply (X& x, Y& b, double reduction, InverseOperatorResult& res) = 0;
136 
138  virtual SolverCategory::Category category() const
139 #ifdef DUNE_ISTL_SUPPORT_OLD_CATEGORY_INTERFACE
140  {
141  DUNE_THROW(Dune::Exception,"It is necessary to implement the category method in a derived classes, in the future this method will pure virtual.");
142  }
143 #else
144  = 0;
145 #endif
146 
148  virtual ~InverseOperator () {}
149 
150  protected:
151  // spacing values
152  enum { iterationSpacing = 5 , normSpacing = 16 };
153 
155  void printHeader(std::ostream& s) const
156  {
157  s << std::setw(iterationSpacing) << " Iter";
158  s << std::setw(normSpacing) << "Defect";
159  s << std::setw(normSpacing) << "Rate" << std::endl;
160  }
161 
163  template <typename CountType, typename DataType>
164  void printOutput(std::ostream& s,
165  const CountType& iter,
166  const DataType& norm,
167  const DataType& norm_old) const
168  {
169  const DataType rate = norm/norm_old;
170  s << std::setw(iterationSpacing) << iter << " ";
171  s << std::setw(normSpacing) << norm << " ";
172  s << std::setw(normSpacing) << rate << std::endl;
173  }
174 
176  template <typename CountType, typename DataType>
177  void printOutput(std::ostream& s,
178  const CountType& iter,
179  const DataType& norm) const
180  {
181  s << std::setw(iterationSpacing) << iter << " ";
182  s << std::setw(normSpacing) << norm << std::endl;
183  }
184  };
185 
194  template<class X, class Y>
195  class IterativeSolver : public InverseOperator<X,Y>{
196  public:
197  using typename InverseOperator<X,Y>::domain_type;
198  using typename InverseOperator<X,Y>::range_type;
199  using typename InverseOperator<X,Y>::field_type;
200  using typename InverseOperator<X,Y>::real_type;
202 
222  IterativeSolver (LinearOperator<X,Y>& op, Preconditioner<X,Y>& prec, scalar_real_type reduction, int maxit, int verbose) :
223  _op(stackobject_to_shared_ptr(op)),
224  _prec(stackobject_to_shared_ptr(prec)),
225  _sp(new SeqScalarProduct<X>),
226  _reduction(reduction), _maxit(maxit), _verbose(verbose), _category(SolverCategory::sequential)
227  {
229  DUNE_THROW(InvalidSolverCategory, "LinearOperator has to be sequential!");
231  DUNE_THROW(InvalidSolverCategory, "Preconditioner has to be sequential!");
232  }
233 
255  scalar_real_type reduction, int maxit, int verbose) :
256  _op(stackobject_to_shared_ptr(op)),
257  _prec(stackobject_to_shared_ptr(prec)),
258  _sp(stackobject_to_shared_ptr(sp)),
259  _reduction(reduction), _maxit(maxit), _verbose(verbose), _category(SolverCategory::category(op))
260  {
262  DUNE_THROW(InvalidSolverCategory, "LinearOperator and Preconditioner must have the same SolverCategory!");
264  DUNE_THROW(InvalidSolverCategory, "LinearOperator and ScalarProduct must have the same SolverCategory!");
265  }
266 
267  // #warning actually we want to have this as the default and just implement the second one
268  // //! \copydoc InverseOperator::apply(X&,Y&,InverseOperatorResult&)
269  // virtual void apply (X& x, Y& b, InverseOperatorResult& res)
270  // {
271  // apply(x,b,_reduction,res);
272  // }
273 
274 #ifndef DOXYGEN
275  // make sure the three-argument apply from the base class does not get shadowed
276  // by the redefined four-argument version below
278 #endif
279 
285  virtual void apply (X& x, X& b, double reduction, InverseOperatorResult& res)
286  {
287  scalar_real_type saved_reduction = _reduction;
288  _reduction = reduction;
289  this->apply(x,b,res);
290  _reduction = saved_reduction;
291  }
292 
295  {
296  return _category;
297  }
298 
299  protected:
300  std::shared_ptr<LinearOperator<X,Y>> _op;
301  std::shared_ptr<Preconditioner<X,Y>> _prec;
302  std::shared_ptr<ScalarProduct<X>> _sp;
304  int _maxit;
305  int _verbose;
307  };
308 
316  template <typename ISTLLinearSolver, typename BCRSMatrix>
318  {
319  public:
320  static void setMatrix (ISTLLinearSolver& solver,
321  const BCRSMatrix& matrix)
322  {
323  static const bool is_direct_solver
327  }
328 
329  protected:
334  template <bool is_direct_solver, typename Dummy = void>
336  {
337  static void setMatrix (ISTLLinearSolver&,
338  const BCRSMatrix&)
339  {}
340  };
341 
346  template <typename Dummy>
347  struct Implementation<true,Dummy>
348  {
349  static void setMatrix (ISTLLinearSolver& solver,
350  const BCRSMatrix& matrix)
351  {
352  solver.setMatrix(matrix);
353  }
354  };
355  };
356 
360 }
361 
362 #endif
Dune::InvalidSolverCategory
Definition: solvercategory.hh:52
Dune::InverseOperatorResult::iterations
int iterations
Number of iterations.
Definition: solver.hh:59
preconditioner.hh
Dune::SolverCategory::Category
Category
Definition: solvercategory.hh:21
Dune::InverseOperatorResult::clear
void clear()
Resets all data.
Definition: solver.hh:49
Dune::InverseOperatorResult::reduction
double reduction
Reduction achieved: .
Definition: solver.hh:62
Dune::SolverCategory::sequential
@ sequential
Category for sequential solvers.
Definition: solvercategory.hh:23
Dune::IterativeSolver::_category
SolverCategory::Category _category
Definition: solver.hh:306
Dune::IterativeSolver
Base class for all implementations of iterative solvers.
Definition: solver.hh:195
Dune::InverseOperator::range_type
Y range_type
Type of the range of the operator to be inverted.
Definition: solver.hh:97
Dune::InverseOperator::field_type
X::field_type field_type
The field type of the operator.
Definition: solver.hh:100
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:103
Dune::InverseOperatorResult::elapsed
double elapsed
Elapsed time in seconds.
Definition: solver.hh:74
Dune::IterativeSolver::_maxit
int _maxit
Definition: solver.hh:304
Dune::IterativeSolver::apply
virtual void apply(X &x, X &b, double reduction, InverseOperatorResult &res)
Apply inverse operator with given reduction factor.
Definition: solver.hh:285
Dune::IterativeSolver::category
virtual SolverCategory::Category category() const
Category of the solver (see SolverCategory::Category)
Definition: solver.hh:294
Dune::InverseOperator::category
virtual SolverCategory::Category category() const =0
Category of the solver (see SolverCategory::Category)
Dune::SolverHelper::Implementation< true, Dummy >::setMatrix
static void setMatrix(ISTLLinearSolver &solver, const BCRSMatrix &matrix)
Definition: solver.hh:349
Dune::InverseOperator
Abstract base class for all solvers.
Definition: solver.hh:91
Dune::InverseOperator::domain_type
X domain_type
Type of the domain of the operator to be inverted.
Definition: solver.hh:94
Dune::InverseOperatorResult::converged
bool converged
True if convergence criterion has been met.
Definition: solver.hh:65
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::apply
virtual void apply(X &x, Y &b, InverseOperatorResult &res)=0
Apply inverse operator,.
Dune::InverseOperator::iterationSpacing
@ iterationSpacing
Definition: solver.hh:152
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:164
scalarproducts.hh
Define base class for scalar product and norm.
Dune::SeqScalarProduct
Default implementation for the scalar case.
Definition: scalarproducts.hh:84
operators.hh
Define general, extensible interface for operators. The available implementation wraps a matrix.
Dune::InverseOperatorResult
Statistics about the application of an inverse operator.
Definition: solver.hh:40
Dune::BCRSMatrix
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:422
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:222
Dune::IterativeSolver::_verbose
int _verbose
Definition: solver.hh:305
Dune::SolverHelper::Implementation::setMatrix
static void setMatrix(ISTLLinearSolver &, const BCRSMatrix &)
Definition: solver.hh:337
Dune::ScalarProduct
Base class for scalar product and norm computation.
Definition: scalarproducts.hh:46
Dune::InverseOperator::scalar_real_type
SimdScalar< real_type > scalar_real_type
scalar type underlying the field_type
Definition: solver.hh:106
Dune::InverseOperator::normSpacing
@ normSpacing
Definition: solver.hh:152
Dune::InverseOperatorResult::InverseOperatorResult
InverseOperatorResult()
Default constructor.
Definition: solver.hh:43
Dune::IsDirectSolver
Definition: solvertype.hh:13
Dune::IterativeSolver::_prec
std::shared_ptr< Preconditioner< X, Y > > _prec
Definition: solver.hh:301
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:254
Dune::InverseOperator::~InverseOperator
virtual ~InverseOperator()
Destructor.
Definition: solver.hh:148
Dune::SolverHelper::setMatrix
static void setMatrix(ISTLLinearSolver &solver, const BCRSMatrix &matrix)
Definition: solver.hh:320
Dune::InverseOperator::printOutput
void printOutput(std::ostream &s, const CountType &iter, const DataType &norm) const
helper function for printing solver output
Definition: solver.hh:177
Dune::SolverHelper::Implementation
Implementation that works together with iterative ISTL solvers, e.g. Dune::CGSolver or Dune::BiCGSTAB...
Definition: solver.hh:335
Dune::IterativeSolver::_sp
std::shared_ptr< ScalarProduct< X > > _sp
Definition: solver.hh:302
Dune
Definition: allocator.hh:7
solvertype.hh
Templates characterizing the type of a solver.
Dune::InverseOperatorResult::condition_estimate
double condition_estimate
Estimate of condition number.
Definition: solver.hh:71
Dune::IterativeSolver::_reduction
scalar_real_type _reduction
Definition: solver.hh:303
Dune::InverseOperatorResult::conv_rate
double conv_rate
Convergence rate (average reduction per step)
Definition: solver.hh:68
Dune::SolverCategory
Categories for the solvers.
Definition: solvercategory.hh:19
Dune::Preconditioner
Base class for matrix free definition of preconditioners.
Definition: preconditioner.hh:30
Dune::LinearOperator
A linear operator.
Definition: operators.hh:64
Dune::SolverHelper
Helper class for notifying a DUNE-ISTL linear solver about a change of the iteration matrix object in...
Definition: solver.hh:317
Dune::InverseOperator::printHeader
void printHeader(std::ostream &s) const
helper function for printing header of solver output
Definition: solver.hh:155
Dune::IterativeSolver::_op
std::shared_ptr< LinearOperator< X, Y > > _op
Definition: solver.hh:300