dune-fem 2.8.0
Loading...
Searching...
No Matches
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
9namespace Dune
10{
11namespace Fem
12{
13namespace 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
double real(const complex< Dune::Fem::Double > &x)
Definition: double.hh:983
Definition: bindguard.hh:11
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