DUNE PDELab (git)

solver.hh
Go to the documentation of this file.
1// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
2// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
3// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
4// vi: set et ts=4 sw=2 sts=2:
5
6#ifndef DUNE_ISTL_SOLVER_HH
7#define DUNE_ISTL_SOLVER_HH
8
9#include <dune-istl-config.hh> // DUNE_ISTL_SUPPORT_OLD_CATEGORY_INTERFACE
10
11#include <iomanip>
12#include <ostream>
13#include <string>
14#include <functional>
15
21#include <dune/common/timer.hh>
22
23#include "solvertype.hh"
24#include "preconditioner.hh"
25#include "operators.hh"
26#include "scalarproducts.hh"
27
28namespace Dune
29{
50 {
53 {
54 clear();
55 }
56
58 void clear ()
59 {
60 iterations = 0;
61 reduction = 0;
62 converged = false;
63 conv_rate = 1;
64 elapsed = 0;
66 }
67
70
72 double reduction;
73
76
78 double conv_rate;
79
81 double condition_estimate = -1;
82
84 double elapsed;
85 };
86
87
88 //=====================================================================
100 template<class X, class Y>
102 public:
104 typedef X domain_type;
105
107 typedef Y range_type;
108
110 typedef typename X::field_type field_type;
111
113 typedef typename FieldTraits<field_type>::real_type real_type;
114
117
130 virtual void apply (X& x, Y& b, InverseOperatorResult& res) = 0;
131
145 virtual void apply (X& x, Y& b, double reduction, InverseOperatorResult& res) = 0;
146
149#ifdef DUNE_ISTL_SUPPORT_OLD_CATEGORY_INTERFACE
150 {
151 DUNE_THROW(Dune::Exception,"It is necessary to implement the category method in a derived classes, in the future this method will pure virtual.");
152 }
153#else
154 = 0;
155#endif
156
158 virtual ~InverseOperator () {}
159
160 protected:
161 // spacing values
162 enum { iterationSpacing = 5 , normSpacing = 16 };
163
165 void printHeader(std::ostream& s) const
166 {
167 s << std::setw(iterationSpacing) << " Iter";
168 s << std::setw(normSpacing) << "Defect";
169 s << std::setw(normSpacing) << "Rate" << std::endl;
170 }
171
173 template <typename CountType, typename DataType>
174 void printOutput(std::ostream& s,
175 const CountType& iter,
176 const DataType& norm,
177 const DataType& norm_old) const
178 {
179 const DataType rate = norm/norm_old;
180 s << std::setw(iterationSpacing) << iter << " ";
181 s << std::setw(normSpacing) << Simd::io(norm) << " ";
182 s << std::setw(normSpacing) << Simd::io(rate) << std::endl;
183 }
184
186 template <typename CountType, typename DataType>
187 void printOutput(std::ostream& s,
188 const CountType& iter,
189 const DataType& norm) const
190 {
191 s << std::setw(iterationSpacing) << iter << " ";
192 s << std::setw(normSpacing) << Simd::io(norm) << std::endl;
193 }
194 };
195
204 template<class X, class Y>
205 class IterativeSolver : public InverseOperator<X,Y>{
206 public:
210 using typename InverseOperator<X,Y>::real_type;
212
232 IterativeSolver (const LinearOperator<X,Y>& op, Preconditioner<X,Y>& prec, scalar_real_type reduction, int maxit, int verbose) :
234 _prec(stackobject_to_shared_ptr(prec)),
235 _sp(new SeqScalarProduct<X>),
236 _reduction(reduction), _maxit(maxit), _verbose(verbose), _category(SolverCategory::sequential)
237 {
239 DUNE_THROW(InvalidSolverCategory, "LinearOperator has to be sequential!");
241 DUNE_THROW(InvalidSolverCategory, "Preconditioner has to be sequential!");
242 }
243
265 scalar_real_type reduction, int maxit, int verbose) :
267 _prec(stackobject_to_shared_ptr(prec)),
269 _reduction(reduction), _maxit(maxit), _verbose(verbose), _category(SolverCategory::category(op))
270 {
272 DUNE_THROW(InvalidSolverCategory, "LinearOperator and Preconditioner must have the same SolverCategory!");
274 DUNE_THROW(InvalidSolverCategory, "LinearOperator and ScalarProduct must have the same SolverCategory!");
275 }
276
292 IterativeSolver (std::shared_ptr<const LinearOperator<X,Y> > op, std::shared_ptr<Preconditioner<X,X> > prec, const ParameterTree& configuration) :
293 IterativeSolver(op,std::make_shared<SeqScalarProduct<X>>(),prec,
294 configuration.get<real_type>("reduction"),
295 configuration.get<int>("maxit"),
296 configuration.get<int>("verbose"))
297 {}
298
315 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) :
316 IterativeSolver(op,sp,prec,
317 configuration.get<scalar_real_type>("reduction"),
318 configuration.get<int>("maxit"),
319 configuration.get<int>("verbose"))
320 {}
321
342 IterativeSolver (std::shared_ptr<const LinearOperator<X,Y>> op,
343 std::shared_ptr<const ScalarProduct<X>> sp,
344 std::shared_ptr<Preconditioner<X,Y>> prec,
345 scalar_real_type reduction, int maxit, int verbose) :
346 _op(op),
347 _prec(prec),
348 _sp(sp),
349 _reduction(reduction), _maxit(maxit), _verbose(verbose),
350 _category(SolverCategory::category(*op))
351 {
353 DUNE_THROW(InvalidSolverCategory, "LinearOperator and Preconditioner must have the same SolverCategory!");
355 DUNE_THROW(InvalidSolverCategory, "LinearOperator and ScalarProduct must have the same SolverCategory!");
356 }
357
358 // #warning actually we want to have this as the default and just implement the second one
359 // //! \copydoc InverseOperator::apply(X&,Y&,InverseOperatorResult&)
360 // virtual void apply (X& x, Y& b, InverseOperatorResult& res)
361 // {
362 // apply(x,b,_reduction,res);
363 // }
364
365#ifndef DOXYGEN
366 // make sure the three-argument apply from the base class does not get shadowed
367 // by the redefined four-argument version below
368 using InverseOperator<X,Y>::apply;
369#endif
370
376 void apply (X& x, X& b, double reduction, InverseOperatorResult& res) override
377 {
378 scalar_real_type saved_reduction = _reduction;
379 _reduction = reduction;
380 this->apply(x,b,res);
381 _reduction = saved_reduction;
382 }
383
386 {
387 return _category;
388 }
389
390 std::string name() const{
391 std::string name = className(*this);
392 return name.substr(0, name.find("<"));
393 }
394
412 template<class CountType = unsigned int>
413 class Iteration {
414 public:
416 : _i(0)
417 , _res(res)
418 , _parent(parent)
419 , _valid(true)
420 {
421 res.clear();
422 if(_parent._verbose>0){
423 std::cout << "=== " << parent.name() << std::endl;
424 if(_parent._verbose > 1)
425 _parent.printHeader(std::cout);
426 }
427 }
428
429 Iteration(const Iteration&) = delete;
430 Iteration(Iteration&& other)
431 : _def0(other._def0)
432 , _def(other._def)
433 , _i(other._i)
434 , _watch(other._watch)
435 , _res(other._res)
436 , _parent(other._parent)
437 , _valid(other._valid)
438 {
439 other._valid = false;
440 }
441
442 ~Iteration(){
443 if(_valid)
444 finalize();
445 }
446
457 bool step(CountType i, real_type def){
458 if (!Simd::allTrue(isFinite(def))) // check for inf or NaN
459 {
460 if (_parent._verbose>0)
461 std::cout << "=== " << _parent.name() << ": abort due to infinite or NaN defect"
462 << std::endl;
464 _parent.name() << ": defect=" << Simd::io(def)
465 << " is infinite or NaN");
466 }
467 if(i == 0)
468 _def0 = def;
469 if(_parent._verbose > 1){
470 if(i!=0)
471 _parent.printOutput(std::cout,i,def,_def);
472 else
473 _parent.printOutput(std::cout,i,def);
474 }
475 _def = def;
476 _i = i;
477 _res.converged = (Simd::allTrue(def<_def0*_parent._reduction || def<real_type(1E-30))); // convergence check
478 return _res.converged;
479 }
480
481 protected:
482 void finalize(){
483 _res.converged = (Simd::allTrue(_def<_def0*_parent._reduction || _def<real_type(1E-30)));
484 _res.iterations = _i;
485 _res.reduction = static_cast<double>(Simd::max(_def/_def0));
486 _res.conv_rate = pow(_res.reduction,1.0/_i);
487 _res.elapsed = _watch.elapsed();
488 if (_parent._verbose>0) // final print
489 {
490 std::cout << "=== rate=" << _res.conv_rate
491 << ", T=" << _res.elapsed
492 << ", TIT=" << _res.elapsed/_res.iterations
493 << ", IT=" << _res.iterations << std::endl;
494 }
495 }
496
497 real_type _def0 = 0.0, _def = 0.0;
498 CountType _i;
499 Timer _watch;
500 InverseOperatorResult& _res;
501 const IterativeSolver& _parent;
502 bool _valid;
503 };
504
505 protected:
506 std::shared_ptr<const LinearOperator<X,Y>> _op;
507 std::shared_ptr<Preconditioner<X,Y>> _prec;
508 std::shared_ptr<const ScalarProduct<X>> _sp;
509 scalar_real_type _reduction;
510 int _maxit;
511 int _verbose;
512 SolverCategory::Category _category;
513 };
514
522 template <typename ISTLLinearSolver, typename BCRSMatrix>
524 {
525 public:
526 static void setMatrix (ISTLLinearSolver& solver,
527 const BCRSMatrix& matrix)
528 {
529 static const bool is_direct_solver
530 = Dune::IsDirectSolver<ISTLLinearSolver>::value;
533 }
534
535 protected:
540 template <bool is_direct_solver, typename Dummy = void>
542 {
543 static void setMatrix (ISTLLinearSolver&,
544 const BCRSMatrix&)
545 {}
546 };
547
552 template <typename Dummy>
553 struct Implementation<true,Dummy>
554 {
555 static void setMatrix (ISTLLinearSolver& solver,
556 const BCRSMatrix& matrix)
557 {
558 solver.setMatrix(matrix);
559 }
560 };
561 };
562
566}
567
568#endif
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:466
Base class for Dune-Exceptions.
Definition: exceptions.hh:96
Abstract base class for all solvers.
Definition: solver.hh:101
void printHeader(std::ostream &s) const
helper function for printing header of solver output
Definition: solver.hh:165
virtual ~InverseOperator()
Destructor.
Definition: solver.hh:158
void printOutput(std::ostream &s, const CountType &iter, const DataType &norm) const
helper function for printing solver output
Definition: solver.hh:187
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
virtual void apply(X &x, Y &b, double reduction, InverseOperatorResult &res)=0
apply inverse operator, with given convergence criteria.
Simd::Scalar< real_type > scalar_real_type
scalar type underlying the field_type
Definition: solver.hh:116
Y range_type
Type of the range of the operator to be inverted.
Definition: solver.hh:107
X domain_type
Type of the domain of the operator to be inverted.
Definition: solver.hh:104
virtual void apply(X &x, Y &b, InverseOperatorResult &res)=0
Apply inverse operator,.
X::field_type field_type
The field type of the operator.
Definition: solver.hh:110
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:113
virtual SolverCategory::Category category() const =0
Category of the solver (see SolverCategory::Category)
Class for controlling iterative methods.
Definition: solver.hh:413
bool step(CountType i, real_type def)
registers the iteration step, checks for invalid defect norm and convergence.
Definition: solver.hh:457
Base class for all implementations of iterative solvers.
Definition: solver.hh:205
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
SolverCategory::Category category() const override
Category of the solver (see SolverCategory::Category)
Definition: solver.hh:385
IterativeSolver(std::shared_ptr< const LinearOperator< X, Y > > op, std::shared_ptr< Preconditioner< X, X > > prec, const ParameterTree &configuration)
Constructor.
Definition: solver.hh:292
void apply(X &x, X &b, double reduction, InverseOperatorResult &res) override
Apply inverse operator with given reduction factor.
Definition: solver.hh:376
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
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
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
A linear operator.
Definition: operators.hh:68
Hierarchical structure of string parameters.
Definition: parametertree.hh:37
Base class for matrix free definition of preconditioners.
Definition: preconditioner.hh:33
Base class for scalar product and norm computation.
Definition: scalarproducts.hh:52
Default implementation for the scalar case.
Definition: scalarproducts.hh:168
Thrown when a solver aborts due to some problem.
Definition: istlexception.hh:46
Helper class for notifying a DUNE-ISTL linear solver about a change of the iteration matrix object in...
Definition: solver.hh:524
double elapsed() const noexcept
Get elapsed user-time from last reset until now/last stop in seconds.
Definition: timer.hh:77
A few common exception classes.
IO interface of the SIMD abstraction.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
auto io(const V &v)
construct a stream inserter
Definition: io.hh:106
bool allTrue(const Mask &mask)
Whether all entries are true
Definition: interface.hh:439
auto max(const V &v1, const V &v2)
The binary maximum value over two simd objects.
Definition: interface.hh:409
typename Overloads::ScalarType< std::decay_t< V > >::type Scalar
Element type of some SIMD type.
Definition: interface.hh:235
auto isFinite(const FieldVector< K, SIZE > &b, PriorityTag< 2 >, ADLTag)
Returns whether all entries are finite.
Definition: fvector.hh:604
Dune namespace.
Definition: alignedallocator.hh:13
std::shared_ptr< T > stackobject_to_shared_ptr(T &t)
Create a shared_ptr for a stack-allocated object.
Definition: shared_ptr.hh:72
std::string className()
Provide the demangled class name of a type T as a string.
Definition: classname.hh:47
constexpr auto get(std::integer_sequence< T, II... >, std::integral_constant< std::size_t, pos >={})
Return the entry at position pos of the given sequence.
Definition: integersequence.hh:22
STL namespace.
Define general, extensible interface for operators. The available implementation wraps a matrix.
A hierarchical structure of string parameters.
Define base class for scalar product and norm.
This file implements several utilities related to std::shared_ptr.
Include file for users of the SIMD abstraction layer.
Templates characterizing the type of a solver.
Statistics about the application of an inverse operator.
Definition: solver.hh:50
InverseOperatorResult()
Default constructor.
Definition: solver.hh:52
double condition_estimate
Estimate of condition number.
Definition: solver.hh:81
double elapsed
Elapsed time in seconds.
Definition: solver.hh:84
int iterations
Number of iterations.
Definition: solver.hh:69
double reduction
Reduction achieved: .
Definition: solver.hh:72
void clear()
Resets all data.
Definition: solver.hh:58
double conv_rate
Convergence rate (average reduction per step)
Definition: solver.hh:78
bool converged
True if convergence criterion has been met.
Definition: solver.hh:75
Categories for the solvers.
Definition: solvercategory.hh:22
Category
Definition: solvercategory.hh:23
@ sequential
Category for sequential solvers.
Definition: solvercategory.hh:25
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
Implementation that works together with iterative ISTL solvers, e.g. Dune::CGSolver or Dune::BiCGSTAB...
Definition: solver.hh:542
A simple timing class.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 24, 23:30, 2024)