Reference documentation for deal.II version 8.1.0
solver_qmrs.h
1 // ---------------------------------------------------------------------
2 // @f$Id: solver_qmrs.h 31349 2013-10-20 19:07:06Z maier @f$
3 //
4 // Copyright (C) 1999 - 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_qmrs_h
18 #define __deal2__solver_qmrs_h
19 
20 #include <deal.II/base/config.h>
21 #include <deal.II/lac/solver.h>
22 #include <deal.II/lac/solver_control.h>
23 #include <deal.II/base/logstream.h>
24 #include <cmath>
25 #include <deal.II/base/subscriptor.h>
26 
27 #include <cmath>
28 
30 
33 
66 template <class VECTOR = Vector<double> >
67 class SolverQMRS : public Solver<VECTOR>
68 {
69 public:
91  {
100  double breakdown=1.e-16) :
103  {}
104 
109 
113  double breakdown;
114  };
115 
121  const AdditionalData &data=AdditionalData());
122 
129  const AdditionalData &data=AdditionalData());
130 
135  template<class MATRIX, class PRECONDITIONER>
136  void
137  solve (const MATRIX &A,
138  VECTOR &x,
139  const VECTOR &b,
140  const PRECONDITIONER &precondition);
141 
151  virtual void print_vectors(const unsigned int step,
152  const VECTOR &x,
153  const VECTOR &r,
154  const VECTOR &d) const;
155 protected:
160  virtual double criterion();
161 
168  VECTOR *Vv;
169  VECTOR *Vp;
170  VECTOR *Vq;
171  VECTOR *Vt;
172  VECTOR *Vd;
176  VECTOR *Vx;
180  const VECTOR *Vb;
181 
192  double res2;
197 private:
201  template<class MATRIX, class PRECONDITIONER>
202  bool
203  iterate(const MATRIX &A, const PRECONDITIONER &precondition);
207  unsigned int step;
208 };
209 
211 /*------------------------- Implementation ----------------------------*/
212 
213 #ifndef DOXYGEN
214 
215 template<class VECTOR>
218  const AdditionalData &data)
219  :
220  Solver<VECTOR>(cn,mem),
221  additional_data(data)
222 {}
223 
224 
225 
226 template<class VECTOR>
228  const AdditionalData &data)
229  :
230  Solver<VECTOR>(cn),
231  additional_data(data)
232 {}
233 
234 
235 
236 template<class VECTOR>
237 double
239 {
240  return std::sqrt(res2);
241 }
242 
243 
244 
245 template<class VECTOR>
246 void
247 SolverQMRS<VECTOR>::print_vectors(const unsigned int,
248  const VECTOR &,
249  const VECTOR &,
250  const VECTOR &) const
251 {}
252 
253 
254 
255 template<class VECTOR>
256 template<class MATRIX, class PRECONDITIONER>
257 void
259  VECTOR &x,
260  const VECTOR &b,
261  const PRECONDITIONER &precondition)
262 {
263  deallog.push("QMRS");
264 
265  // Memory allocation
266  Vv = this->memory.alloc();
267  Vp = this->memory.alloc();
268  Vq = this->memory.alloc();
269  Vt = this->memory.alloc();
270  Vd = this->memory.alloc();
271 
272  Vx = &x;
273  Vb = &b;
274  // resize the vectors, but do not set
275  // the values since they'd be overwritten
276  // soon anyway.
277  Vv->reinit(x, true);
278  Vp->reinit(x, true);
279  Vq->reinit(x, true);
280  Vt->reinit(x, true);
281 
282  step = 0;
283 
284  bool state;
285 
286  do
287  {
288  if (step)
289  deallog << "Restart step " << step << std::endl;
290  state = iterate(A, precondition);
291  }
292  while (state);
293 
294  // Deallocate Memory
295  this->memory.free(Vv);
296  this->memory.free(Vp);
297  this->memory.free(Vq);
298  this->memory.free(Vt);
299  this->memory.free(Vd);
300 
301  // Output
302  deallog.pop();
303 
304  // in case of failure: throw exception
305  if (this->control().last_check() != SolverControl::success)
306  AssertThrow(false, SolverControl::NoConvergence (this->control().last_step(),
307  this->control().last_value()));
308  // otherwise exit as normal
309 }
310 
311 
312 
313 template<class VECTOR>
314 template<class MATRIX, class PRECONDITIONER>
315 bool
317  const PRECONDITIONER &precondition)
318 {
319  /* Remark: the matrix A in the article is the preconditioned matrix.
320  * Therefore, we have to precondition x before we compute the first residual.
321  * In step 1 we replace p by q to avoid one preconditioning step.
322  * There are still two steps left, making this algorithm expensive.
323  */
324 
326 
327  // define some aliases for simpler access
328  VECTOR &v = *Vv;
329  VECTOR &p = *Vp;
330  VECTOR &q = *Vq;
331  VECTOR &t = *Vt;
332  VECTOR &d = *Vd;
333  VECTOR &x = *Vx;
334  const VECTOR &b = *Vb;
335 
336  int it=0;
337 
338  double tau, rho, theta=0, sigma, alpha, psi, theta_old, rho_old, beta;
339  double res;
340 
341  d.reinit(x);
342 
343  // Apply right preconditioning to x
344  precondition.vmult(q,x);
345  // Preconditioned residual
346  A.vmult(v,q);
347  v.sadd(-1.,1.,b);
348  res = v.l2_norm();
349 
350  if (this->control().check(step, res) == SolverControl::success)
351  return false;
352 
353  p = v;
354 
355  precondition.vmult(q,p);
356 
357  tau = v.norm_sqr();
358  rho = q*v;
359 
360  while (state == SolverControl::iterate)
361  {
362  step++;
363  it++;
364  // Step 1
365  A.vmult(t,q);
366  // Step 2
367  sigma = q*t;
368 
369 //TODO:[?] Find a really good breakdown criterion. The absolute one detects breakdown instead of convergence
370  if (std::fabs(sigma/rho) < additional_data.breakdown)
371  return true;
372  // Step 3
373  alpha = rho/sigma;
374 
375  v.add(-alpha,t);
376  // Step 4
377  theta_old = theta;
378  theta = v*v/tau;
379  psi = 1./(1.+theta);
380  tau *= theta*psi;
381 
382  d.sadd(psi*theta_old, psi*alpha, p);
383  x.add(d);
384 
385  print_vectors(step,x,v,d);
386  // Step 5
387  if (additional_data.exact_residual)
388  {
389  A.vmult(q,x);
390  q.sadd(-1.,1.,b);
391  res = q.l2_norm();
392  }
393  else
394  res = std::sqrt((it+1)*tau);
395  state = this->control().check(step,res);
396  if ((state == SolverControl::success)
397  || (state == SolverControl::failure))
398  return false;
399  // Step 6
400  if (std::fabs(rho) < additional_data.breakdown)
401  return true;
402  // Step 7
403  rho_old = rho;
404  precondition.vmult(q,v);
405  rho = q*v;
406 
407  beta = rho/rho_old;
408  p.sadd(beta,v);
409  precondition.vmult(q,p);
410  }
411  return false;
412 }
413 
414 #endif // DOXYGEN
415 
416 DEAL_II_NAMESPACE_CLOSE
417 
418 #endif
void pop()
const VECTOR * Vb
Definition: solver_qmrs.h:180
void vmult(VECTOR &u, const VECTOR &v) const
Continue iteration.
virtual void print_vectors(const unsigned int step, const VECTOR &x, const VECTOR &r, const VECTOR &d) const
AdditionalData additional_data
Definition: solver_qmrs.h:196
#define AssertThrow(cond, exc)
Definition: exceptions.h:362
VECTOR * Vx
Definition: solver_qmrs.h:176
virtual double criterion()
VECTOR * Vv
Definition: solver_qmrs.h:168
bool iterate(const MATRIX &A, const PRECONDITIONER &precondition)
Definition: solver.h:147
Stop iteration, goal not reached.
void push(const std::string &text)
unsigned int step
Definition: solver_qmrs.h:207
double res2
Definition: solver_qmrs.h:192
void solve(const MATRIX &A, VECTOR &x, const VECTOR &b, const PRECONDITIONER &precondition)
AdditionalData(bool exact_residual=false, double breakdown=1.e-16)
Definition: solver_qmrs.h:99
SolverQMRS(SolverControl &cn, VectorMemory< VECTOR > &mem, const AdditionalData &data=AdditionalData())
Stop iteration, goal reached.