dune-fem  2.8-git
cg.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_CG_HH
2 #define DUNE_FEM_CG_HH
3 
4 #include <type_traits>
5 #include <vector>
6 
7 #include <dune/common/ftraits.hh>
8 
9 namespace Dune
10 {
11 namespace Fem
12 {
13 namespace LinearSolver
14 {
16  {
17  static const int absolute = 0;
18  static const int relative = 1;
19  static const int residualReduction = 2;
20  };
21 
22 
23  template <class Operator, class Precoditioner, class DiscreteFunction>
24  inline int cg( Operator &op,
25  Precoditioner* preconditioner,
26  std::vector< DiscreteFunction >& tempMem,
27  DiscreteFunction& x,
28  const DiscreteFunction& b,
29  const double epsilon,
30  const int maxIterations,
31  const int toleranceCriteria,
32  std::ostream* os = nullptr )
33  {
34  typedef typename DiscreteFunction::RangeFieldType RangeFieldType;
35  typedef typename Dune::FieldTraits< RangeFieldType >::real_type RealType;
36 
37  const RealType tolerance = (epsilon * epsilon) * b.normSquaredDofs( );
38 
39  assert( preconditioner ? tempMem.size() == 5 : tempMem.size() == 3 );
40 
41  DiscreteFunction& h = tempMem[ 0 ];
42 
43  //h=Ax
44  op( x, h );
45 
46  //r=Ax-b
47  DiscreteFunction& r = tempMem[ 1 ];
48  r.assign( h );
49  r -= b;
50 
51  //p=b-A*x <= r_0 Deufelhard
52  DiscreteFunction& p = tempMem[ 2 ];
53  p.assign( b );
54  p -= h;
55 
56  DiscreteFunction& s = ( preconditioner ) ? tempMem[ 3 ] : p;
57 
58  //q=B*p <=q Deufelhard
59  DiscreteFunction& q = ( preconditioner ) ? tempMem[ 4 ] : p;
60  if( preconditioner )
61  {
62  (*preconditioner)( p, q );
63  s.assign( q );
64  }
65 
66  RangeFieldType prevResiduum = 0; // note that these will be real_type but require scalar product evaluation
67  RangeFieldType residuum = p.scalarProductDofs( q );//<p,Bp>
68 
69  int iterations = 0;
70  for( iterations = 0; (std::real(residuum) > tolerance) && (iterations < maxIterations); ++iterations )
71  {
72  if( iterations > 0 )
73  {
74  assert( residuum/prevResiduum == residuum/prevResiduum );
75  const RangeFieldType beta= residuum / prevResiduum;
76  q *= beta;
77  if( preconditioner )
78  {
79  q += s;
80  }
81  else
82  {
83  p -= r;
84  }
85  }
86 
87  op( q, h );
88 
89  RangeFieldType qdoth = q.scalarProductDofs( h );
90  const RangeFieldType alpha = residuum / qdoth;//<p,Bp>/<q,Aq>
91  assert( !std::isnan( alpha ) );
92  x.axpy( alpha, q );
93 
94  if( preconditioner )
95  {
96  p.axpy( -alpha, h );//r_k
97  (*preconditioner)( p, s); //B*r_k
98 
99  prevResiduum = residuum;//<rk-1,B*rk-1>
100  residuum = p.scalarProductDofs( s );//<rk,B*rk>
101  }
102  else
103  {
104  r.axpy( alpha, h );
105 
106  prevResiduum = residuum;
107  residuum = r.normSquaredDofs( );
108  }
109 
110  if( os )
111  {
112  (*os) << "Fem::CG it: " << iterations << " : sqr(residuum) " << residuum << std::endl;
113  }
114  }
115 
116  return (iterations < maxIterations) ? iterations : -iterations;
117  }
118 
119 } // naemspace LinearSolver
120 } // namespace Fem
121 
122 } // namespace Dune
123 
124 #endif // #ifndef DUNE_FEM_CGINVERSEOPERATOR_HH
Definition: bindguard.hh:11
double real(const std::complex< Double > &x)
Definition: double.hh:906
int cg(Operator &op, Precoditioner *preconditioner, std::vector< DiscreteFunction > &tempMem, DiscreteFunction &x, const DiscreteFunction &b, const double epsilon, const int maxIterations, const int toleranceCriteria, std::ostream *os=nullptr)
Definition: cg.hh:24
abstract operator
Definition: operator.hh:34
static const int absolute
Definition: cg.hh:17
static const int relative
Definition: cg.hh:18
static const int residualReduction
Definition: cg.hh:19