dune-fem  2.6-git
cginverseoperator.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_CGINVERSEOPERATOR_HH
2 #define DUNE_FEM_CGINVERSEOPERATOR_HH
3 
4 #include <limits>
5 #include <memory>
6 #include <type_traits>
7 
12 
13 namespace Dune
14 {
15 
16  namespace Fem
17  {
18 
19  // ConjugateGradientSolver
20  // -----------------------
21 
28  template< class Operator>
30  {
32 
33  public:
36 
41  typedef typename Dune::FieldTraits< RangeFieldType >::real_type RealType;
42 
47 
50 
51 
52  private:
53  static_assert( (std::is_same< DomainFunctionType, RangeFunctionType >::value),
54  "DomainFunctionType must equal RangeFunctionType." );
55 
56  public:
64  unsigned int maxIterations,
65  bool verbose,
66  const ParameterReader &parameter = Parameter::container() )
67  : epsilon_( epsilon ),
68  maxIterations_( maxIterations ),
69  verbose_( verbose ),
70  averageCommTime_( 0.0 ),
71  realCount_( 0 )
72 
73  {}
74 
81  unsigned int maxIterations,
82  const ParameterReader &parameter = Parameter::container() )
83  : epsilon_( epsilon ),
84  maxIterations_( maxIterations ),
85  verbose_( parameter.getValue< bool >( "fem.solver.verbose", false ) ),
86  averageCommTime_( 0.0 ),
87  realCount_( 0 )
88  {}
89 
90  // ConjugateGradientSolver ( const ThisType & )=delete;
91 
103  void solve ( const OperatorType &op, const RangeFunctionType &b, DomainFunctionType &x ) const;
104 
117  void solve ( const OperatorType &op, const PreconditionerType &p,
118  const RangeFunctionType &b, DomainFunctionType &x ) const;
119 
121  unsigned int iterations () const
122  {
123  return realCount_;
124  }
125 
126  void setMaxIterations ( unsigned int maxIterations ) { maxIterations_ = maxIterations; }
127 
128 
130  double averageCommTime() const
131  {
132  return averageCommTime_;
133  }
134 
135  protected:
137  unsigned int maxIterations_;
138  const bool verbose_;
139  mutable double averageCommTime_;
140  mutable unsigned int realCount_;
141  };
142 
143 
144  namespace Solver
145  {
146  // CGInverseOperator
147  // -----------------
148 
154  template< class DiscreteFunction >
156  : public Fem::Operator< DiscreteFunction, DiscreteFunction >
157  {
159 
160  public:
163 
166 
168  typedef typename Dune::FieldTraits< RangeFieldType >::real_type RealType;
169 
170  private:
172 
173  public:
181  CGInverseOperator ( RealType redEps, RealType absLimit,
182  unsigned int maxIter, bool verbose,
183  const ParameterReader &parameter = Parameter::container() )
184  : solver_( absLimit, maxIter, verbose, parameter ),
185  parameter_( parameter )
186  {}
187 
188 
195  CGInverseOperator ( RealType redEps, RealType absLimit,
196  unsigned int maxIter,
197  const ParameterReader &parameter = Parameter::container() )
198  : solver_( absLimit, maxIter, parameter ),
199  parameter_( parameter )
200  {}
201 
202  CGInverseOperator ( RealType redEps, RealType absLimit,
203  const ParameterReader &parameter = Parameter::container() )
204  : CGInverseOperator( redEps, absLimit, std::numeric_limits< unsigned int >::max(), parameter )
205  {}
206 
216  RealType redEps, RealType absLimit,
217  unsigned int maxIter, bool verbose,
218  const ParameterReader &parameter = Parameter::container() )
219  : CGInverseOperator( redEps, absLimit, maxIter, verbose, parameter )
220  {
221  bind( op );
222  }
223 
224 
233  RealType redEps, RealType absLimit,
234  unsigned int maxIter,
235  const ParameterReader &parameter = Parameter::container() )
236  : CGInverseOperator( redEps, absLimit, maxIter, parameter )
237  {
238  bind( op );
239  }
240 
242  RealType redEps, RealType absLimit,
243  const ParameterReader &parameter = Parameter::container() )
244  : CGInverseOperator( redEps, absLimit, parameter )
245  {
246  bind( op );
247  }
248 
259  const PreconditionerType &precond,
260  RealType redEps, RealType absLimit,
261  unsigned int maxIter, bool verbose,
262  const ParameterReader &parameter = Parameter::container() )
263  : CGInverseOperator( redEps, absLimit, maxIter, verbose )
264  {
265  bind( op, precond );
266  }
267 
277  const PreconditionerType &precond,
278  RealType redEps, RealType absLimit,
279  const ParameterReader &parameter = Parameter::container() )
280  : CGInverseOperator( redEps, absLimit, parameter )
281  {
282  bind( op, precond );
283  }
284 
286  const PreconditionerType &precond,
287  RealType redEps, RealType absLimit,
288  unsigned int maxIter,
289  const ParameterReader &parameter = Parameter::container() )
290  : CGInverseOperator( redEps, absLimit, maxIter, parameter )
291  {
292  bind( op, precond );
293  }
294 
295  void bind ( const OperatorType &op ) { operator_ = &op; preconditioner_ = nullptr; }
296  void bind ( const OperatorType &op, const PreconditionerType &precond )
297  {
298  operator_ = &op;
299  preconditioner_ = &precond;
300  }
301  void unbind () { operator_ = nullptr; preconditioner_ = nullptr; }
302 
311  virtual void operator()( const DomainFunctionType &arg, RangeFunctionType &dest ) const
312  {
313  prepare();
314  apply(arg,dest);
315  finalize();
316  }
317 
318  template<typename... A>
319  inline void prepare(A... ) const
320  {}
321 
322  inline void finalize() const
323  {}
324 
333  virtual void apply( const DomainFunctionType &arg, RangeFunctionType &dest ) const
334  {
335  assert(operator_);
336  if(preconditioner_)
337  solver_.solve( *operator_, *preconditioner_, arg, dest );
338  else
339  solver_.solve( *operator_, arg, dest );
340  }
341 
343  unsigned int iterations () const
344  {
345  return solver_.iterations();
346  }
347 
348  void setMaxIterations ( unsigned int maxIter ) { solver_.setMaxIterations( maxIter ); }
349 
351  double averageCommTime() const
352  {
353  return solver_.averageCommTime();
354  }
355 
356  protected:
357  const OperatorType *operator_ = nullptr;
361  };
362  }
363 
364  // CGInverseOperator
365  // -----------------
366 
374  template< class DiscreteFunction,
377  : public Fem::Solver::CGInverseOperator< DiscreteFunction >
378  {
380 
381  public:
384 
386 
388  typedef Op OperatorType;
389 
391  typedef typename Dune::FieldTraits< RangeFieldType >::real_type RealType;
392 
393  // Preconditioner is to approximate op^-1 !
395 
403  CGInverseOperator ( RealType redEps, RealType absLimit,
404  unsigned int maxIter, bool verbose,
405  const ParameterReader &parameter = Parameter::container() )
406  : BaseType( redEps, absLimit, maxIter, verbose, parameter )
407  {}
408 
415  CGInverseOperator ( RealType redEps, RealType absLimit,
416  unsigned int maxIter,
417  const ParameterReader &parameter = Parameter::container() )
418  : BaseType( redEps, absLimit, maxIter, parameter )
419  {}
420 
421  CGInverseOperator ( RealType redEps, RealType absLimit,
422  const ParameterReader &parameter = Parameter::container() )
423  : BaseType( redEps, absLimit, parameter )
424  {}
425 
434  template< class LinearOperator, std::enable_if_t< std::is_base_of< OperatorType, LinearOperator >::value, int > = 0 >
436  RealType redEps, RealType absLimit,
437  unsigned int maxIter, bool verbose,
438  const ParameterReader &parameter = Parameter::container() )
439  : BaseType( redEps, absLimit, maxIter, verbose, parameter )
440  {
441  bind( op );
442  }
443 
451  template< class LinearOperator, std::enable_if_t< std::is_base_of< OperatorType, LinearOperator >::value, int > = 0 >
453  RealType redEps, RealType absLimit,
454  unsigned int maxIter,
455  const ParameterReader &parameter = Parameter::container() )
456  : BaseType( redEps, absLimit, maxIter, parameter )
457  {
458  bind( op );
459  }
460 
461  template< class LinearOperator, std::enable_if_t< std::is_base_of< OperatorType, LinearOperator >::value, int > = 0 >
463  RealType redEps, RealType absLimit,
464  const ParameterReader &parameter = Parameter::container() )
465  : BaseType( redEps, absLimit, parameter )
466  {
467  bind( op );
468  }
469 
479  const PreconditioningType &precond,
480  RealType redEps, RealType absLimit,
481  unsigned int maxIter, bool verbose,
482  const ParameterReader &parameter = Parameter::container() )
483  : BaseType( op, precond, redEps, absLimit, maxIter, verbose, parameter )
484  {}
485 
495  const PreconditioningType &precond,
496  RealType redEps, RealType absLimit,
497  unsigned int maxIter,
498  const ParameterReader &parameter = Parameter::container() )
499  : BaseType( op, precond, redEps, absLimit, maxIter, parameter )
500  {}
501 
503  const PreconditioningType &precond,
504  RealType redEps, RealType absLimit,
505  const ParameterReader &parameter = Parameter::container() )
506  : BaseType( op, precond, redEps, absLimit, parameter )
507  {}
508 
509 
510 
511  template< class LinearOperator, std::enable_if_t< std::is_base_of< OperatorType, LinearOperator >::value, int > = 0 >
512  void bind ( const LinearOperator &op )
513  {
514  BaseType::bind( op );
515  checkPreconditioning( op );
516  }
517  void unbind ()
518  {
520  precondObj_.reset();
521  }
522 
523  protected:
524  template< class LinearOperator >
525  void checkPreconditioning( const LinearOperator &linearOp )
526  {
527  const bool preconditioning = parameter_.template getValue< bool >( "fem.preconditioning", false );
528  if( preconditioning && std::is_base_of< AssembledOperator< DomainFunctionType, DomainFunctionType > ,LinearOperator > :: value )
529  {
530  // create diagonal preconditioner
533  }
534  }
535 
537  using BaseType::parameter_;
538  std::unique_ptr< PreconditioningType > precondObj_;
539  };
540 
541  // Implementation of ConjugateGradientSolver
542  // -----------------------------------------
543 
544  template< class Operator >
547  {
548  const bool verbose = (verbose_ && (b.space().gridPart().comm().rank() == 0));
549 
550  const RealType tolerance = (epsilon_ * epsilon_) * b.normSquaredDofs( );
551 
552  averageCommTime_ = 0.0;
553 
554  RangeFunctionType h( b );
555  op( x, h );
556 
557  RangeFunctionType r( h );
558  r -= b;
559 
560  RangeFunctionType p( b );
561  p -= h;
562 
563  RealType prevResiduum = 0;
564  RealType residuum = r.normSquaredDofs( );
565 
566  for( realCount_ = 0; (residuum > tolerance) && (realCount_ < maxIterations_); ++realCount_ )
567  {
568  if( realCount_ > 0 )
569  {
570  assert( residuum/prevResiduum == residuum/prevResiduum );
571  p *= (residuum / prevResiduum);
572  p -= r;
573  }
574 
575  op( p, h );
576 
577  RangeFieldType pdoth = p.scalarProductDofs( h );
578  const RangeFieldType alpha = residuum / pdoth;
579  assert( alpha == alpha );
580  x.axpy( alpha, p );
581  r.axpy( alpha, h );
582 
583  prevResiduum = residuum;
584  residuum = r.normSquaredDofs( );
585 
586  double exchangeTime = h.space().communicator().exchangeTime();
587  if( verbose )
588  {
589  std::cerr << "CG-Iteration: " << realCount_ << ", sqr(Residuum): " << residuum << std::endl;
590  // only for parallel apps
591  if( b.space().gridPart().comm().size() > 1 )
592  std::cerr << "Communication needed: " << exchangeTime << " s" << std::endl;
593  }
594 
595  averageCommTime_ += exchangeTime;
596  }
597  }
598 
599 
600  template< class Operator >
602  ::solve ( const OperatorType &op, const PreconditionerType &precond, const RangeFunctionType &b, DomainFunctionType &x ) const
603  {
604  const bool verbose = (verbose_ && (b.space().gridPart().comm().rank() == 0));
605 
606  const RealType tolerance = (epsilon_ * epsilon_) * b.normSquaredDofs( );
607 
608  averageCommTime_ = 0.0;
609 
610  RangeFunctionType h( b );
611  //h=Ax
612  op( x, h );
613 
614  //r=Ax-b
615  RangeFunctionType r( h );
616  r -= b;
617 
618  //p=b-A*x <= r_0 Deufelhard
619  RangeFunctionType p( b );
620  p -= h;
621 
622  //q=B*p <=q Deuf
623  RangeFunctionType q ( b );
624  precond(p,q);
625 
626  RangeFunctionType s (q);
627 
628  RangeFieldType prevResiduum = 0; // note that these will be real_type but require scalar product evaluation
629  RangeFieldType residuum = p.scalarProductDofs( q );//<p,Bp>
630 
631  for( realCount_ = 0; (std::real(residuum) > tolerance) && (realCount_ < maxIterations_); ++realCount_ )
632  {
633  if( realCount_ > 0 )
634  {
635  assert( residuum/prevResiduum == residuum/prevResiduum );
636  const RangeFieldType beta=residuum/prevResiduum;
637  q*=beta;
638  q+=(s);
639  }
640 
641  op( q, h );
642 
643  RangeFieldType qdoth = q.scalarProductDofs( h );
644  const RangeFieldType alpha = residuum / qdoth;//<p,Bp>/<q,Aq>
645  assert( alpha == alpha );
646  x.axpy( alpha, q );
647 
648  p.axpy( -alpha, h );//r_k
649 
650  precond(p,s); //B*r_k
651 
652  prevResiduum = residuum;//<rk-1,B*rk-1>
653 
654  residuum = p.scalarProductDofs( s );//<rk,B*rk>
655 
656  double exchangeTime = h.space().communicator().exchangeTime();
657  if( verbose )
658  {
659  std::cerr << "CG-Iteration: " << realCount_ << ", Residuum: " << residuum << std::endl;
660  // only for parallel apps
661  if( b.space().gridPart().comm().size() > 1 )
662  std::cerr << "Communication needed: " << exchangeTime << " s" << std::endl;
663  }
664 
665  averageCommTime_ += exchangeTime;
666  }
667  }
668 
669  } // namespace Fem
670 
671 } // namespace Dune
672 
673 #endif // #ifndef DUNE_FEM_CGINVERSEOPERATOR_HH
Definition: bindguard.hh:11
double real(const std::complex< Double > &x)
Definition: double.hh:907
static double max(const Double &v, const double p)
Definition: double.hh:399
static DUNE_EXPORT ParameterContainer & container()
Definition: io/parameter.hh:192
abstract operator
Definition: operator.hh:26
DomainFunction DomainFunctionType
type of discrete function in the operator's domain
Definition: operator.hh:28
RangeFunction::RangeFieldType RangeFieldType
field type of the operator's range
Definition: operator.hh:35
DomainFunction::RangeFieldType DomainFieldType
field type of the operator's domain
Definition: operator.hh:33
RangeFunction RangeFunctionType
type of discrete function in the operator's range
Definition: operator.hh:30
abstract affine-linear operator
Definition: operator.hh:73
abstract matrix operator
Definition: operator.hh:110
linear solver using the CG algorithm
Definition: cginverseoperator.hh:30
ConjugateGradientSolver(const RealType &epsilon, unsigned int maxIterations, bool verbose, const ParameterReader &parameter=Parameter::container())
constructor
Definition: cginverseoperator.hh:63
OperatorType::RangeFunctionType RangeFunctionType
type of the operator's range vectors
Definition: cginverseoperator.hh:46
ConjugateGradientSolver(RealType epsilon, unsigned int maxIterations, const ParameterReader &parameter=Parameter::container())
constructor
Definition: cginverseoperator.hh:80
Dune::FieldTraits< RangeFieldType >::real_type RealType
Definition: cginverseoperator.hh:41
unsigned int iterations() const
number of iterations needed for last solve
Definition: cginverseoperator.hh:121
void solve(const OperatorType &op, const RangeFunctionType &b, DomainFunctionType &x) const
solve
Definition: cginverseoperator.hh:546
double averageCommTime() const
return average communication time during last solve
Definition: cginverseoperator.hh:130
void setMaxIterations(unsigned int maxIterations)
Definition: cginverseoperator.hh:126
void solve(const OperatorType &op, const PreconditionerType &p, const RangeFunctionType &b, DomainFunctionType &x) const
solve
Definition: cginverseoperator.hh:602
Fem::Operator< RangeFunctionType, DomainFunctionType > PreconditionerType
type of the preconditioner, maps from the range of the operator (the dual space) in it's domain
Definition: cginverseoperator.hh:49
const bool verbose_
Definition: cginverseoperator.hh:138
OperatorType::RangeFieldType RangeFieldType
field type of the operator's range vectors
Definition: cginverseoperator.hh:40
OperatorType::DomainFieldType DomainFieldType
field type of the operator's domain vectors
Definition: cginverseoperator.hh:38
double averageCommTime_
Definition: cginverseoperator.hh:139
OperatorType::DomainFunctionType DomainFunctionType
type of the operator's domain vectors
Definition: cginverseoperator.hh:44
const RealType epsilon_
Definition: cginverseoperator.hh:136
Operator OperatorType
type of the operators to invert
Definition: cginverseoperator.hh:35
unsigned int realCount_
Definition: cginverseoperator.hh:140
unsigned int maxIterations_
Definition: cginverseoperator.hh:137
Inverse operator base on CG method. This is the base class for the cg solver and does not imvolve any...
Definition: cginverseoperator.hh:157
CGInverseOperator(const OperatorType &op, RealType redEps, RealType absLimit, unsigned int maxIter, const ParameterReader &parameter=Parameter::container())
constructor of CGInverseOperator
Definition: cginverseoperator.hh:232
const PreconditionerType * preconditioner_
Definition: cginverseoperator.hh:358
void prepare(A...) const
Definition: cginverseoperator.hh:319
void unbind()
Definition: cginverseoperator.hh:301
CGInverseOperator(RealType redEps, RealType absLimit, unsigned int maxIter, bool verbose, const ParameterReader &parameter=Parameter::container())
constructor of CGInverseOperator
Definition: cginverseoperator.hh:181
virtual void operator()(const DomainFunctionType &arg, RangeFunctionType &dest) const
application operator
Definition: cginverseoperator.hh:311
const OperatorType * operator_
Definition: cginverseoperator.hh:357
CGInverseOperator(const OperatorType &op, const PreconditionerType &precond, RealType redEps, RealType absLimit, const ParameterReader &parameter=Parameter::container())
constructor of CGInverseOperator
Definition: cginverseoperator.hh:276
CGInverseOperator(const OperatorType &op, RealType redEps, RealType absLimit, const ParameterReader &parameter=Parameter::container())
Definition: cginverseoperator.hh:241
CGInverseOperator(RealType redEps, RealType absLimit, unsigned int maxIter, const ParameterReader &parameter=Parameter::container())
constructor of CGInverseOperator
Definition: cginverseoperator.hh:195
Fem::Operator< DomainFunctionType, RangeFunctionType > OperatorType
Definition: cginverseoperator.hh:164
unsigned int iterations() const
number of iterations needed for last solve
Definition: cginverseoperator.hh:343
virtual void apply(const DomainFunctionType &arg, RangeFunctionType &dest) const
application operator
Definition: cginverseoperator.hh:333
void bind(const OperatorType &op, const PreconditionerType &precond)
Definition: cginverseoperator.hh:296
CGInverseOperator(const OperatorType &op, const PreconditionerType &precond, RealType redEps, RealType absLimit, unsigned int maxIter, const ParameterReader &parameter=Parameter::container())
Definition: cginverseoperator.hh:285
CGInverseOperator(RealType redEps, RealType absLimit, const ParameterReader &parameter=Parameter::container())
Definition: cginverseoperator.hh:202
OperatorType::RangeFieldType RangeFieldType
Definition: cginverseoperator.hh:167
ParameterReader parameter_
Definition: cginverseoperator.hh:360
BaseType::DomainFunctionType DomainFunctionType
Definition: cginverseoperator.hh:161
CGInverseOperator(const OperatorType &op, const PreconditionerType &precond, RealType redEps, RealType absLimit, unsigned int maxIter, bool verbose, const ParameterReader &parameter=Parameter::container())
constructor of CGInverseOperator
Definition: cginverseoperator.hh:258
double averageCommTime() const
return average communication time during last solve
Definition: cginverseoperator.hh:351
Fem::Operator< RangeFunctionType, DomainFunctionType > PreconditionerType
Definition: cginverseoperator.hh:165
Dune::FieldTraits< RangeFieldType >::real_type RealType
Definition: cginverseoperator.hh:168
BaseType::RangeFunctionType RangeFunctionType
Definition: cginverseoperator.hh:162
SolverType solver_
Definition: cginverseoperator.hh:359
void setMaxIterations(unsigned int maxIter)
Definition: cginverseoperator.hh:348
void bind(const OperatorType &op)
Definition: cginverseoperator.hh:295
CGInverseOperator(const OperatorType &op, RealType redEps, RealType absLimit, unsigned int maxIter, bool verbose, const ParameterReader &parameter=Parameter::container())
constructor of CGInverseOperator
Definition: cginverseoperator.hh:215
void finalize() const
Definition: cginverseoperator.hh:322
Inverse operator base on CG method. Uses a runtime parameter fem.preconditioning which enables diagon...
Definition: cginverseoperator.hh:378
CGInverseOperator(RealType redEps, RealType absLimit, unsigned int maxIter, bool verbose, const ParameterReader &parameter=Parameter::container())
constructor of CGInverseOperator
Definition: cginverseoperator.hh:403
CGInverseOperator(const OperatorType &op, const PreconditioningType &precond, RealType redEps, RealType absLimit, const ParameterReader &parameter=Parameter::container())
Definition: cginverseoperator.hh:502
Op OperatorType
type of operator
Definition: cginverseoperator.hh:388
CGInverseOperator(RealType redEps, RealType absLimit, const ParameterReader &parameter=Parameter::container())
Definition: cginverseoperator.hh:421
CGInverseOperator(const LinearOperator &op, RealType redEps, RealType absLimit, unsigned int maxIter, const ParameterReader &parameter=Parameter::container())
constructor of CGInverseOperator
Definition: cginverseoperator.hh:452
Dune::FieldTraits< RangeFieldType >::real_type RealType
Definition: cginverseoperator.hh:391
void bind(const LinearOperator &op)
Definition: cginverseoperator.hh:512
CGInverseOperator(RealType redEps, RealType absLimit, unsigned int maxIter, const ParameterReader &parameter=Parameter::container())
constructor of CGInverseOperator
Definition: cginverseoperator.hh:415
CGInverseOperator(const LinearOperator &op, RealType redEps, RealType absLimit, const ParameterReader &parameter=Parameter::container())
Definition: cginverseoperator.hh:462
OperatorType::RangeFieldType RangeFieldType
Definition: cginverseoperator.hh:390
CGInverseOperator(const OperatorType &op, const PreconditioningType &precond, RealType redEps, RealType absLimit, unsigned int maxIter, bool verbose, const ParameterReader &parameter=Parameter::container())
constructor of CGInverseOperator
Definition: cginverseoperator.hh:478
Fem::Operator< RangeFunctionType, DomainFunctionType > PreconditioningType
Definition: cginverseoperator.hh:394
CGInverseOperator(const LinearOperator &op, RealType redEps, RealType absLimit, unsigned int maxIter, bool verbose, const ParameterReader &parameter=Parameter::container())
constructor of CGInverseOperator
Definition: cginverseoperator.hh:435
BaseType::RangeFunctionType RangeFunctionType
Definition: cginverseoperator.hh:383
CGInverseOperator(const OperatorType &op, const PreconditioningType &precond, RealType redEps, RealType absLimit, unsigned int maxIter, const ParameterReader &parameter=Parameter::container())
constructor of CGInverseOperator
Definition: cginverseoperator.hh:494
DomainFunctionType DestinationType
Definition: cginverseoperator.hh:385
void checkPreconditioning(const LinearOperator &linearOp)
Definition: cginverseoperator.hh:525
BaseType::DomainFunctionType DomainFunctionType
Definition: cginverseoperator.hh:382
std::unique_ptr< PreconditioningType > precondObj_
Definition: cginverseoperator.hh:538
void unbind()
Definition: cginverseoperator.hh:517
Precondtioner, multiplies with inverse of the diagonal works with.
Definition: diagonalpreconditioner.hh:152