Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
solver_qmrs.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2019 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_solver_qmrs_h
17 #define dealii_solver_qmrs_h
18 
19 #include <deal.II/base/config.h>
20 
21 #include <deal.II/base/logstream.h>
22 #include <deal.II/base/subscriptor.h>
23 
24 #include <deal.II/lac/solver.h>
25 #include <deal.II/lac/solver_control.h>
26 
27 #include <cmath>
28 
29 DEAL_II_NAMESPACE_OPEN
30 
33 
94 template <typename VectorType = Vector<double>>
95 class SolverQMRS : public SolverBase<VectorType>
96 {
97 public:
123  {
130  explicit AdditionalData(const bool left_preconditioning = false,
131  const double solver_tolerance = 1.e-9,
132  const bool breakdown_testing = true,
133  const double breakdown_threshold = 1.e-16)
138  {}
139 
144 
149 
154 
160  };
161 
167  const AdditionalData & data = AdditionalData());
168 
174 
178  template <typename MatrixType, typename PreconditionerType>
179  void
180  solve(const MatrixType & A,
181  VectorType & x,
182  const VectorType & b,
183  const PreconditionerType &preconditioner);
184 
190  virtual void
191  print_vectors(const unsigned int step,
192  const VectorType & x,
193  const VectorType & r,
194  const VectorType & d) const;
195 
196 protected:
201 
202 private:
208  {
209  SolverControl::State state;
210  double last_residual;
211 
213  const double last_residual);
214  };
215 
220  template <typename MatrixType, typename PreconditionerType>
222  iterate(const MatrixType & A,
223  VectorType & x,
224  const VectorType & b,
225  const PreconditionerType &preconditioner,
226  VectorType & r,
227  VectorType & u,
228  VectorType & q,
229  VectorType & t,
230  VectorType & d);
231 
235  unsigned int step;
236 };
237 
239 /*------------------------- Implementation ----------------------------*/
240 
241 #ifndef DOXYGEN
242 
243 
244 template <class VectorType>
246  const SolverControl::State state,
247  const double last_residual)
248  : state(state)
249  , last_residual(last_residual)
250 {}
251 
252 
253 
254 template <class VectorType>
257  const AdditionalData & data)
258  : SolverBase<VectorType>(cn, mem)
259  , additional_data(data)
260  , step(0)
261 {}
262 
263 template <class VectorType>
265  const AdditionalData &data)
266  : SolverBase<VectorType>(cn)
267  , additional_data(data)
268  , step(0)
269 {}
270 
271 template <class VectorType>
272 void
273 SolverQMRS<VectorType>::print_vectors(const unsigned int,
274  const VectorType &,
275  const VectorType &,
276  const VectorType &) const
277 {}
278 
279 template <class VectorType>
280 template <typename MatrixType, typename PreconditionerType>
281 void
282 SolverQMRS<VectorType>::solve(const MatrixType & A,
283  VectorType & x,
284  const VectorType & b,
285  const PreconditionerType &preconditioner)
286 {
287  LogStream::Prefix prefix("SQMR");
288 
289 
290  // temporary vectors, allocated trough the @p VectorMemory object at the
291  // start of the actual solution process and deallocated at the end.
292  typename VectorMemory<VectorType>::Pointer Vr(this->memory);
293  typename VectorMemory<VectorType>::Pointer Vu(this->memory);
294  typename VectorMemory<VectorType>::Pointer Vq(this->memory);
295  typename VectorMemory<VectorType>::Pointer Vt(this->memory);
296  typename VectorMemory<VectorType>::Pointer Vd(this->memory);
297 
298 
299  // resize the vectors, but do not set
300  // the values since they'd be overwritten
301  // soon anyway.
302  Vr->reinit(x, true);
303  Vu->reinit(x, true);
304  Vq->reinit(x, true);
305  Vt->reinit(x, true);
306  Vd->reinit(x, true);
307 
308  step = 0;
309 
310  IterationResult state(SolverControl::failure, 0);
311 
312  do
313  {
314  if (step > 0)
315  deallog << "Restart step " << step << std::endl;
316  state = iterate(A, x, b, preconditioner, *Vr, *Vu, *Vq, *Vt, *Vd);
317  }
318  while (state.state == SolverControl::iterate);
319 
320 
321  // in case of failure: throw exception
322  AssertThrow(state.state == SolverControl::success,
323  SolverControl::NoConvergence(step, state.last_residual));
324  // otherwise exit as normal
325 }
326 
327 template <class VectorType>
328 template <typename MatrixType, typename PreconditionerType>
330 SolverQMRS<VectorType>::iterate(const MatrixType & A,
331  VectorType & x,
332  const VectorType & b,
333  const PreconditionerType &preconditioner,
334  VectorType & r,
335  VectorType & u,
336  VectorType & q,
337  VectorType & t,
338  VectorType & d)
339 {
341 
342  int it = 0;
343 
344  double tau, rho, theta = 0;
345  double res;
346 
347  // Compute the start residual
348  A.vmult(r, x);
349  r.sadd(-1., 1., b);
350 
351  // Doing the initial preconditioning
352  if (additional_data.left_preconditioning)
353  {
354  // Left preconditioning
355  preconditioner.vmult(t, r);
356  q = t;
357  }
358  else
359  {
360  // Right preconditioning
361  t = r;
362  preconditioner.vmult(q, t);
363  }
364 
365  tau = t.norm_sqr();
366  res = std::sqrt(tau);
367 
368  if (this->iteration_status(step, res, x) == SolverControl::success)
369  return IterationResult(SolverControl::success, res);
370 
371  rho = q * r;
372 
373  while (state == SolverControl::iterate)
374  {
375  step++;
376  it++;
377  //--------------------------------------------------------------
378  // Step 1: apply the system matrix and compute one inner product
379  //--------------------------------------------------------------
380  A.vmult(t, q);
381  const double sigma = q * t;
382 
383  // Check the breakdown criterion
384  if (additional_data.breakdown_testing == true &&
385  std::fabs(sigma) < additional_data.breakdown_threshold)
386  return IterationResult(SolverControl::iterate, res);
387  // Update the residual
388  const double alpha = rho / sigma;
389  r.add(-alpha, t);
390 
391  //--------------------------------------------------------------
392  // Step 2: update the solution vector
393  //--------------------------------------------------------------
394  const double theta_old = theta;
395 
396  // Apply the preconditioner
397  if (additional_data.left_preconditioning)
398  {
399  // Left Preconditioning
400  preconditioner.vmult(t, r);
401  }
402  else
403  {
404  // Right Preconditioning
405  t = r;
406  }
407 
408  // Double updates
409  theta = t * t / tau;
410  const double psi = 1. / (1. + theta);
411  tau *= theta * psi;
412 
413  // Actual update of the solution vector
414  d.sadd(psi * theta_old, psi * alpha, q);
415  x += d;
416 
417  print_vectors(step, x, r, d);
418 
419  // Check for convergence
420  // Compute a simple and cheap upper bound of the norm of the residual
421  // vector b-Ax
422  res = std::sqrt((it + 1) * tau);
423  // If res lies close enough, within the desired tolerance, calculate the
424  // exact residual
425  if (res < additional_data.solver_tolerance)
426  {
427  A.vmult(u, x);
428  u.sadd(-1., 1., b);
429  res = u.l2_norm();
430  }
431  state = this->iteration_status(step, res, x);
432  if ((state == SolverControl::success) ||
433  (state == SolverControl::failure))
434  return IterationResult(state, res);
435 
436  //--------------------------------------------------------------
437  // Step 3: check breakdown criterion and update the vectors
438  //--------------------------------------------------------------
439  if (additional_data.breakdown_testing == true &&
440  std::fabs(sigma) < additional_data.breakdown_threshold)
441  return IterationResult(SolverControl::iterate, res);
442 
443  const double rho_old = rho;
444 
445  // Applying the preconditioner
446  if (additional_data.left_preconditioning)
447  {
448  // Left preconditioning
449  u = t;
450  }
451  else
452  {
453  // Right preconditioning
454  preconditioner.vmult(u, t);
455  }
456 
457  // Double and vector updates
458  rho = u * r;
459  const double beta = rho / rho_old;
460  q.sadd(beta, 1., u);
461  }
462  return IterationResult(SolverControl::success, res);
463 }
464 
465 #endif // DOXYGEN
466 
467 DEAL_II_NAMESPACE_CLOSE
468 
469 #endif
Stop iteration, goal not reached.
Continue iteration.
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
virtual void print_vectors(const unsigned int step, const VectorType &x, const VectorType &r, const VectorType &d) const
SolverQMRS(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
Stop iteration, goal reached.
unsigned int step
Definition: solver_qmrs.h:235
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
AdditionalData(const bool left_preconditioning=false, const double solver_tolerance=1.e-9, const bool breakdown_testing=true, const double breakdown_threshold=1.e-16)
Definition: solver_qmrs.h:130
IterationResult iterate(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner, VectorType &r, VectorType &u, VectorType &q, VectorType &t, VectorType &d)
AdditionalData additional_data
Definition: solver_qmrs.h:200