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;
303  scalar_real_type _reduction;
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
Base class for all implementations of iterative solvers.
Definition: solver.hh:195
static void setMatrix(ISTLLinearSolver &solver, const BCRSMatrix &matrix)
Definition: solver.hh:349
X::field_type field_type
The field type of the operator.
Definition: solver.hh:100
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
Implementation that works together with iterative ISTL solvers, e.g. Dune::CGSolver or Dune::BiCGSTAB...
Definition: solver.hh:335
Definition: allocator.hh:7
double reduction
Reduction achieved: .
Definition: solver.hh:62
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:422
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
SimdScalar< real_type > scalar_real_type
scalar type underlying the field_type
Definition: solver.hh:106
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
Categories for the solvers.
Definition: solvercategory.hh:19
std::shared_ptr< ScalarProduct< X > > _sp
Definition: solver.hh:302
void printHeader(std::ostream &s) const
helper function for printing header of solver output
Definition: solver.hh:155
A linear operator.
Definition: operators.hh:64
Define base class for scalar product and norm.
void printOutput(std::ostream &s, const CountType &iter, const DataType &norm) const
helper function for printing solver output
Definition: solver.hh:177
std::shared_ptr< LinearOperator< X, Y > > _op
Definition: solver.hh:300
scalar_real_type _reduction
Definition: solver.hh:303
SolverCategory::Category _category
Definition: solver.hh:306
double elapsed
Elapsed time in seconds.
Definition: solver.hh:74
virtual ~InverseOperator()
Destructor.
Definition: solver.hh:148
Define general, extensible interface for operators. The available implementation wraps a matrix...
static void setMatrix(ISTLLinearSolver &, const BCRSMatrix &)
Definition: solver.hh:337
static void setMatrix(ISTLLinearSolver &solver, const BCRSMatrix &matrix)
Definition: solver.hh:320
bool converged
True if convergence criterion has been met.
Definition: solver.hh:65
Abstract base class for all solvers.
Definition: solver.hh:91
Definition: solvercategory.hh:52
InverseOperatorResult()
Default constructor.
Definition: solver.hh:43
Default implementation for the scalar case.
Definition: scalarproducts.hh:84
int _maxit
Definition: solver.hh:304
Category
Definition: solvercategory.hh:21
Definition: solvertype.hh:13
X domain_type
Type of the domain of the operator to be inverted.
Definition: solver.hh:94
void clear()
Resets all data.
Definition: solver.hh:49
Base class for scalar product and norm computation.
Definition: scalarproducts.hh:46
virtual void apply(X &x, X &b, double reduction, InverseOperatorResult &res)
Apply inverse operator with given reduction factor.
Definition: solver.hh:285
Y range_type
Type of the range of the operator to be inverted.
Definition: solver.hh:97
virtual SolverCategory::Category category() const
Category of the solver (see SolverCategory::Category)
Definition: solver.hh:294
Helper class for notifying a DUNE-ISTL linear solver about a change of the iteration matrix object in...
Definition: solver.hh:317
Statistics about the application of an inverse operator.
Definition: solver.hh:40
Category for sequential solvers.
Definition: solvercategory.hh:23
Base class for matrix free definition of preconditioners.
Definition: preconditioner.hh:30
int _verbose
Definition: solver.hh:305
int iterations
Number of iterations.
Definition: solver.hh:59
std::shared_ptr< Preconditioner< X, Y > > _prec
Definition: solver.hh:301
Templates characterizing the type of a solver.
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
double conv_rate
Convergence rate (average reduction per step)
Definition: solver.hh:68
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
double condition_estimate
Estimate of condition number.
Definition: solver.hh:71