Reference documentation for deal.II version 8.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
solver_minres.h
1 // ---------------------------------------------------------------------
2 // @f$Id: solver_minres.h 31349 2013-10-20 19:07:06Z maier @f$
3 //
4 // Copyright (C) 2000 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 #ifndef __deal2__solver_minres_h
18 #define __deal2__solver_minres_h
19 
20 
21 #include <deal.II/base/config.h>
22 #include <deal.II/lac/solver.h>
23 #include <deal.II/lac/solver_control.h>
24 #include <deal.II/base/logstream.h>
25 #include <cmath>
26 #include <deal.II/base/subscriptor.h>
27 
29 
32 
62 template <class VECTOR = Vector<double> >
63 class SolverMinRes : public Solver<VECTOR>
64 {
65 public:
73  {
74  };
75 
81  const AdditionalData &data=AdditionalData());
82 
89  const AdditionalData &data=AdditionalData());
90 
94  virtual ~SolverMinRes ();
95 
100  template<class MATRIX, class PRECONDITIONER>
101  void
102  solve (const MATRIX &A,
103  VECTOR &x,
104  const VECTOR &b,
105  const PRECONDITIONER &precondition);
106 
113  DeclException0 (ExcPreconditionerNotDefinite);
115 
116 protected:
121  virtual double criterion();
131  virtual void print_vectors(const unsigned int step,
132  const VECTOR &x,
133  const VECTOR &r,
134  const VECTOR &d) const;
135 
142  VECTOR *Vu0, *Vu1, *Vu2;
143  VECTOR *Vm0, *Vm1, *Vm2;
144  VECTOR *Vv;
145 
156  double res2;
157 };
158 
160 /*------------------------- Implementation ----------------------------*/
161 
162 #ifndef DOXYGEN
163 
164 template<class VECTOR>
167  const AdditionalData &)
168  :
169  Solver<VECTOR>(cn,mem)
170 {}
171 
172 
173 
174 template<class VECTOR>
176  const AdditionalData &)
177  :
178  Solver<VECTOR>(cn)
179 {}
180 
181 
182 template<class VECTOR>
184 {}
185 
186 
187 
188 template<class VECTOR>
189 double
191 {
192  return res2;
193 }
194 
195 
196 template<class VECTOR>
197 void
198 SolverMinRes<VECTOR>::print_vectors(const unsigned int,
199  const VECTOR &,
200  const VECTOR &,
201  const VECTOR &) const
202 {}
203 
204 
205 
206 template<class VECTOR>
207 template<class MATRIX, class PRECONDITIONER>
208 void
210  VECTOR &x,
211  const VECTOR &b,
212  const PRECONDITIONER &precondition)
213 {
215 
216  deallog.push("minres");
217 
218  // Memory allocation
219  Vu0 = this->memory.alloc();
220  Vu1 = this->memory.alloc();
221  Vu2 = this->memory.alloc();
222  Vv = this->memory.alloc();
223  Vm0 = this->memory.alloc();
224  Vm1 = this->memory.alloc();
225  Vm2 = this->memory.alloc();
226  // define some aliases for simpler access
227  typedef VECTOR *vecptr;
228  vecptr u[3] = {Vu0, Vu1, Vu2};
229  vecptr m[3] = {Vm0, Vm1, Vm2};
230  VECTOR &v = *Vv;
231  // resize the vectors, but do not set
232  // the values since they'd be overwritten
233  // soon anyway.
234  u[0]->reinit(b,true);
235  u[1]->reinit(b,true);
236  u[2]->reinit(b,true);
237  m[0]->reinit(b,true);
238  m[1]->reinit(b,true);
239  m[2]->reinit(b,true);
240  v.reinit(b,true);
241 
242  // some values needed
243  double delta[3] = { 0, 0, 0 };
244  double f[2] = { 0, 0 };
245  double e[2] = { 0, 0 };
246 
247  double r_l2 = 0;
248  double r0 = 0;
249  double tau = 0;
250  double c = 0;
251  double gamma = 0;
252  double s = 0;
253  double d_ = 0;
254  double d = 0;
255 
256  // The iteration step.
257  unsigned int j = 1;
258 
259 
260  // Start of the solution process
261  A.vmult(*m[0],x);
262  *u[1] = b;
263  *u[1] -= *m[0];
264  // Precondition is applied.
265  // The preconditioner has to be
266  // positiv definite and symmetric
267 
268  // M v = u[1]
269  precondition.vmult (v,*u[1]);
270 
271  delta[1] = v * (*u[1]);
272  // Preconditioner positive
273  Assert (delta[1]>=0, ExcPreconditionerNotDefinite());
274 
275  r0 = std::sqrt(delta[1]);
276  r_l2 = r0;
277 
278 
279  u[0]->reinit(b);
280  delta[0] = 1.;
281  m[0]->reinit(b);
282  m[1]->reinit(b);
283  m[2]->reinit(b);
284 
285  conv = this->control().check(0,r_l2);
286 
287  while (conv==SolverControl::iterate)
288  {
289  if (delta[1]!=0)
290  v *= 1./std::sqrt(delta[1]);
291  else
292  v.reinit(b);
293 
294  A.vmult(*u[2],v);
295  u[2]->add (-std::sqrt(delta[1]/delta[0]), *u[0]);
296 
297  gamma = *u[2] * v;
298  u[2]->add (-gamma / std::sqrt(delta[1]), *u[1]);
299  *m[0] = v;
300 
301  // precondition: solve M v = u[2]
302  // Preconditioner has to be positiv
303  // definite and symmetric.
304  precondition.vmult(v,*u[2]);
305 
306  delta[2] = v * (*u[2]);
307 
308  Assert (delta[2]>=0, ExcPreconditionerNotDefinite());
309 
310  if (j==1)
311  {
312  d_ = gamma;
313  e[1] = std::sqrt(delta[2]);
314  }
315  if (j>1)
316  {
317  d_ = s * e[0] - c * gamma;
318  e[0] = c * e[0] + s * gamma;
319  f[1] = s * std::sqrt(delta[2]);
320  e[1] = -c * std::sqrt(delta[2]);
321  }
322 
323  d = std::sqrt (d_*d_ + delta[2]);
324 
325  if (j>1) tau *= s / c;
326  c = d_ / d;
327  tau *= c;
328 
329  s = std::sqrt(delta[2]) / d;
330 
331  if (j==1)
332  tau = r0 * c;
333 
334  m[0]->add (-e[0], *m[1]);
335  if (j>1)
336  m[0]->add (-f[0], *m[2]);
337  *m[0] *= 1./d;
338  x.add (tau, *m[0]);
339  r_l2 *= std::fabs(s);
340 
341  conv = this->control().check(j,r_l2);
342 
343  // next iteration step
344  ++j;
345  // All vectors have to be shifted
346  // one iteration step.
347  // This should be changed one time.
348  //
349  // the previous code was like this:
350  // m[2] = m[1];
351  // m[1] = m[0];
352  // but it can be made more efficient,
353  // since the value of m[0] is no more
354  // needed in the next iteration
355  swap (*m[2], *m[1]);
356  swap (*m[1], *m[0]);
357 
358  // likewise, but reverse direction:
359  // u[0] = u[1];
360  // u[1] = u[2];
361  swap (*u[0], *u[1]);
362  swap (*u[1], *u[2]);
363 
364  // these are scalars, so need
365  // to bother
366  f[0] = f[1];
367  e[0] = e[1];
368  delta[0] = delta[1];
369  delta[1] = delta[2];
370  }
371 
372  // Deallocation of Memory
373  this->memory.free(Vu0);
374  this->memory.free(Vu1);
375  this->memory.free(Vu2);
376  this->memory.free(Vv);
377  this->memory.free(Vm0);
378  this->memory.free(Vm1);
379  this->memory.free(Vm2);
380  // Output
381  deallog.pop ();
382 
383  // in case of failure: throw exception
384  if (this->control().last_check() != SolverControl::success)
385  AssertThrow(false, SolverControl::NoConvergence (this->control().last_step(),
386  this->control().last_value()));
387  // otherwise exit as normal
388 }
389 
390 #endif // DOXYGEN
391 
392 DEAL_II_NAMESPACE_CLOSE
393 
394 #endif
void pop()
void vmult(VECTOR &u, const VECTOR &v) const
Continue iteration.
SolverMinRes(SolverControl &cn, VectorMemory< VECTOR > &mem, const AdditionalData &data=AdditionalData())
virtual void print_vectors(const unsigned int step, const VECTOR &x, const VECTOR &r, const VECTOR &d) const
#define AssertThrow(cond, exc)
Definition: exceptions.h:362
DeclException0(ExcPreconditionerNotDefinite)
VECTOR * Vu0
#define Assert(cond, exc)
Definition: exceptions.h:299
virtual ~SolverMinRes()
void swap(BlockVector &u, BlockVector &v)
Definition: solver.h:147
void push(const std::string &text)
void solve(const MATRIX &A, VECTOR &x, const VECTOR &b, const PRECONDITIONER &precondition)
virtual double criterion()
Stop iteration, goal reached.