dune-fem  2.8-git
bicgstab.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_BICGSTAB_HH
2 #define DUNE_FEM_BICGSTAB_HH
3 
4 #include <cmath>
5 #include <cassert>
6 #include <iostream>
7 
8 #include <utility>
9 
11 
12 namespace Dune
13 {
14 namespace Fem
15 {
16 namespace LinearSolver
17 {
18 
19  template <class DiscreteFunction, class FieldType>
20  void scalarProductVecs( const DiscreteFunction& r_df,
21  const DiscreteFunction& s_df,
22  const DiscreteFunction& r_star_df,
23  FieldType* global_dot )
24  {
25  const auto& comm = r_df.space().gridPart().comm();
26 
27  const auto& r = r_df.dofVector();
28  const auto& s = s_df.dofVector();
29  const auto& r_star = r_star_df.dofVector();
30 
31  const auto& auxiliaryDofs = r_df.space().auxiliaryDofs();
32 
33  global_dot[ 4 ] = global_dot[ 3 ] = global_dot[ 2 ] = global_dot[ 1 ] = global_dot[ 0 ] = 0.0;
34 
35  const size_t numAuxiliarys = auxiliaryDofs.size();
36  for( size_t auxiliary = 0, i = 0 ; auxiliary < numAuxiliarys; ++auxiliary )
37  {
38  const size_t nextAuxiliary = auxiliaryDofs[ auxiliary ];
39  for(; i < nextAuxiliary; ++i )
40  {
41  global_dot[ 0 ] += r[ i ] * s[ i ];
42  global_dot[ 1 ] += r[ i ] * r[ i ];
43  global_dot[ 2 ] += s[ i ] * s[ i ];
44  global_dot[ 3 ] += s[ i ] * r_star[ i ];
45  global_dot[ 4 ] += r[ i ] * r_star[ i ];
46  }
47  }
48 
49  // make sum global
50  comm.sum( global_dot, 5 );
51  }
52 
53  /* General implementation of a BiCG-stab algorithm based on Dune::Fem::DiscreteFunction
54  * \param op linear operator A
55  * \param x solution to be sought
56  * \param b right hand side
57  * \param tempMem temporary memory
58  * \param tolerance for the solver
59  * \param maxIter maximal iterations to be performed
60  * \param toleranceCriteria tolerance citeria (see iterativesolvers.hh)
61  * \param os output if enabled
62  */
63  template <class Operator, class Precoditioner, class DiscreteFunction>
64  inline int bicgstab( Operator &op,
65  Precoditioner* preconditioner,
66  std::vector< DiscreteFunction >& tempMem,
67  DiscreteFunction& x,
68  const DiscreteFunction& b,
69  const double tolerance,
70  const int maxIterations,
71  const int toleranceCriteria,
72  std::ostream* os = nullptr )
73  {
74  assert( ( preconditioner ) ? tempMem.size() == 6 : tempMem.size() == 5 );
75 
76  DiscreteFunction& r = tempMem[ 0 ];
77  DiscreteFunction& r_star = tempMem[ 1 ];
78  DiscreteFunction& p = tempMem[ 2 ];
79  DiscreteFunction& s = tempMem[ 3 ];
80  DiscreteFunction& tmp = tempMem[ 4 ];
81 
82  DiscreteFunction& z = ( preconditioner ) ? tempMem[ 5 ] : r;
83 
84  typedef typename DiscreteFunction :: RangeFieldType FieldType;
85 
86  // relative or absolute tolerance
87  FieldType global_dot[5];
88  double _tolerance = tolerance;
89 
90  if (toleranceCriteria == ToleranceCriteria::relative)
91  {
92  global_dot[ 0 ] = b.scalarProductDofs( b );
93  _tolerance *= std::sqrt(global_dot[0]);
94  }
95 
96  // init
97  if( preconditioner ) // right preconditioning
98  {
99  (*preconditioner)(x, z); // z = M^{-1} x
100  op( z, r); // r = A z
101  }
102  else
103  {
104  op(x, r); // r = A x
105  }
106 
107  // residual: r = b - r
108  r *= -1.0;
109  r += b ;
110 
111  // p = r
112  p.assign( r );
113  // r_star = r
114  r_star.assign( r );
115 
116  global_dot[0] = r.scalarProductDofs( r_star );
117 
118  FieldType nu = global_dot[0];
119  if (toleranceCriteria == ToleranceCriteria::residualReduction)
120  {
121  _tolerance *= std::sqrt(nu);
122  }
123 
124  // iterate
125  int iterations = 0;
126  while (true)
127  {
128  // 2x linear operator 1x dot
129  if( preconditioner ) // right preconditioning
130  {
131  (*preconditioner)(p, z); // z = M^{-1} p
132  op( z, tmp); // tmp = A z
133  }
134  else
135  {
136  op(p, tmp); // tmp = A p
137  }
138 
139  global_dot[0] = tmp.scalarProductDofs( r_star );
140 
141  const FieldType alpha = nu / global_dot[0];
142  s.assign( r );
143  s.axpy( -alpha, tmp );
144 
145  if( preconditioner )
146  {
147  // right preconditioning
148  (*preconditioner)(s, z); // z = M^{-1} s
149  op(z, r); // r = A z
150  }
151  else
152  {
153  op(s, r); // r = A s
154  }
155 
156  // 5x dot: r * s | r * r | s * s | s * r_star | r * r_star
157  scalarProductVecs( r, s, r_star, global_dot );
158 
159  // scalars
160  const FieldType omega = global_dot[0] / global_dot[1];
161  const FieldType res = std::sqrt(global_dot[2]
162  -omega*(2.0*global_dot[0] - omega*global_dot[1]) );
163 
164  const FieldType beta = (global_dot[3] - omega*global_dot[4])
165  *alpha / (omega*nu);
166 
167  nu = (global_dot[3] - omega*global_dot[4]);
168 
169  // update
170  // x += alpha * p + omega s
171  x.axpy( alpha, p );
172  x.axpy( omega, s );
173 
174  ++iterations;
175  if (res < _tolerance || iterations >= maxIterations ) break;
176 
177  // r = s - omega r
178  r *= -omega;
179  r += s ;
180 
181  // p = r + beta * ( p - omega * tmp )
182  p *= beta;
183  p.axpy( -omega*beta, tmp );
184  p += r;
185 
186  if ( os )
187  {
188  (*os) << "Fem::BiCGstab it: " << iterations << " : " << res << std::endl;
189  }
190  }
191 
192  // output
193  if ( os )
194  {
195  (*os) << "Fem::BiCGstab: number of iterations: "
196  << iterations
197  << std::endl;
198  }
199 
200  // setup approx solution for right preconditioner use
201  if( preconditioner )
202  {
203  // right preconditioner
204  (*preconditioner)(x, z);
205 
206  // x = z
207  x.assign( z );
208  }
209 
210  if( iterations >= maxIterations )
211  return -iterations;
212 
213  return iterations;
214  }
215 
216 } // end namespace Solver
217 
218 } // end namespace Fem
219 
220 } // end namespace Dune
221 
222 #endif
Definition: bindguard.hh:11
static double sqrt(const Double &v)
Definition: double.hh:886
void scalarProductVecs(const DiscreteFunction &r_df, const DiscreteFunction &s_df, const DiscreteFunction &r_star_df, FieldType *global_dot)
Definition: bicgstab.hh:20
int bicgstab(Operator &op, Precoditioner *preconditioner, std::vector< DiscreteFunction > &tempMem, DiscreteFunction &x, const DiscreteFunction &b, const double tolerance, const int maxIterations, const int toleranceCriteria, std::ostream *os=nullptr)
Definition: bicgstab.hh:64
abstract operator
Definition: operator.hh:34
static const int relative
Definition: cg.hh:18
static const int residualReduction
Definition: cg.hh:19