Reference documentation for deal.II version GIT 7deb6c54a6 2023-06-09 18:50:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
solver_qmrs.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2020 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>
23 
24 #include <deal.II/lac/solver.h>
26 
27 #include <cmath>
28 
30 
93 template <typename VectorType = Vector<double>>
94 class SolverQMRS : public SolverBase<VectorType>
95 {
96 public:
122  {
129  explicit AdditionalData(const bool left_preconditioning = false,
130  const double solver_tolerance = 1.e-9,
131  const bool breakdown_testing = true,
132  const double breakdown_threshold = 1.e-16)
137  {}
138 
143 
148 
153 
159  };
160 
166  const AdditionalData & data = AdditionalData());
167 
173 
177  template <typename MatrixType, typename PreconditionerType>
178  void
179  solve(const MatrixType & A,
180  VectorType & x,
181  const VectorType & b,
182  const PreconditionerType &preconditioner);
183 
189  virtual void
190  print_vectors(const unsigned int step,
191  const VectorType & x,
192  const VectorType & r,
193  const VectorType & d) const;
194 
195 protected:
200 
201 private:
207  {
210 
212  const double last_residual);
213  };
214 
219  template <typename MatrixType, typename PreconditionerType>
221  iterate(const MatrixType & A,
222  VectorType & x,
223  const VectorType & b,
224  const PreconditionerType &preconditioner,
225  VectorType & r,
226  VectorType & u,
227  VectorType & q,
228  VectorType & t,
229  VectorType & d);
230 
234  unsigned int step;
235 };
236 
238 /*------------------------- Implementation ----------------------------*/
239 
240 #ifndef DOXYGEN
241 
242 
243 template <class VectorType>
245  const SolverControl::State state,
246  const double last_residual)
247  : state(state)
248  , last_residual(last_residual)
249 {}
250 
251 
252 
253 template <class VectorType>
256  const AdditionalData & data)
257  : SolverBase<VectorType>(cn, mem)
258  , additional_data(data)
259  , step(0)
260 {}
261 
262 template <class VectorType>
264  const AdditionalData &data)
265  : SolverBase<VectorType>(cn)
266  , additional_data(data)
267  , step(0)
268 {}
269 
270 template <class VectorType>
271 void
272 SolverQMRS<VectorType>::print_vectors(const unsigned int,
273  const VectorType &,
274  const VectorType &,
275  const VectorType &) const
276 {}
277 
278 template <class VectorType>
279 template <typename MatrixType, typename PreconditionerType>
280 void
281 SolverQMRS<VectorType>::solve(const MatrixType & A,
282  VectorType & x,
283  const VectorType & b,
284  const PreconditionerType &preconditioner)
285 {
286  LogStream::Prefix prefix("SQMR");
287 
288 
289  // temporary vectors, allocated through the @p VectorMemory object at the
290  // start of the actual solution process and deallocated at the end.
291  typename VectorMemory<VectorType>::Pointer Vr(this->memory);
292  typename VectorMemory<VectorType>::Pointer Vu(this->memory);
293  typename VectorMemory<VectorType>::Pointer Vq(this->memory);
294  typename VectorMemory<VectorType>::Pointer Vt(this->memory);
295  typename VectorMemory<VectorType>::Pointer Vd(this->memory);
296 
297 
298  // resize the vectors, but do not set
299  // the values since they'd be overwritten
300  // soon anyway.
301  Vr->reinit(x, true);
302  Vu->reinit(x, true);
303  Vq->reinit(x, true);
304  Vt->reinit(x, true);
305  Vd->reinit(x, true);
306 
307  step = 0;
308 
309  IterationResult state(SolverControl::failure, 0);
310 
311  do
312  {
313  if (step > 0)
314  deallog << "Restart step " << step << std::endl;
315  state = iterate(A, x, b, preconditioner, *Vr, *Vu, *Vq, *Vt, *Vd);
316  }
317  while (state.state == SolverControl::iterate);
318 
319 
320  // in case of failure: throw exception
321  AssertThrow(state.state == SolverControl::success,
322  SolverControl::NoConvergence(step, state.last_residual));
323  // otherwise exit as normal
324 }
325 
326 template <class VectorType>
327 template <typename MatrixType, typename PreconditionerType>
329 SolverQMRS<VectorType>::iterate(const MatrixType & A,
330  VectorType & x,
331  const VectorType & b,
332  const PreconditionerType &preconditioner,
333  VectorType & r,
334  VectorType & u,
335  VectorType & q,
336  VectorType & t,
337  VectorType & d)
338 {
340 
341  int it = 0;
342 
343  double tau, rho, theta = 0;
344  double res;
345 
346  // Compute the start residual
347  A.vmult(r, x);
348  r.sadd(-1., 1., b);
349 
350  // Doing the initial preconditioning
351  if (additional_data.left_preconditioning)
352  {
353  // Left preconditioning
354  preconditioner.vmult(t, r);
355  q = t;
356  }
357  else
358  {
359  // Right preconditioning
360  t = r;
361  preconditioner.vmult(q, t);
362  }
363 
364  tau = t.norm_sqr();
365  res = std::sqrt(tau);
366 
367  if (this->iteration_status(step, res, x) == SolverControl::success)
368  return IterationResult(SolverControl::success, res);
369 
370  rho = q * r;
371 
372  while (state == SolverControl::iterate)
373  {
374  step++;
375  it++;
376  //--------------------------------------------------------------
377  // Step 1: apply the system matrix and compute one inner product
378  //--------------------------------------------------------------
379  A.vmult(t, q);
380  const double sigma = q * t;
381 
382  // Check the breakdown criterion
383  if (additional_data.breakdown_testing == true &&
384  std::fabs(sigma) < additional_data.breakdown_threshold)
385  return IterationResult(SolverControl::iterate, res);
386  // Update the residual
387  const double alpha = rho / sigma;
388  r.add(-alpha, t);
389 
390  //--------------------------------------------------------------
391  // Step 2: update the solution vector
392  //--------------------------------------------------------------
393  const double theta_old = theta;
394 
395  // Apply the preconditioner
396  if (additional_data.left_preconditioning)
397  {
398  // Left Preconditioning
399  preconditioner.vmult(t, r);
400  }
401  else
402  {
403  // Right Preconditioning
404  t = r;
405  }
406 
407  // Double updates
408  theta = t * t / tau;
409  const double psi = 1. / (1. + theta);
410  tau *= theta * psi;
411 
412  // Actual update of the solution vector
413  d.sadd(psi * theta_old, psi * alpha, q);
414  x += d;
415 
416  print_vectors(step, x, r, d);
417 
418  // Check for convergence
419  // Compute a simple and cheap upper bound of the norm of the residual
420  // vector b-Ax
421  res = std::sqrt((it + 1) * tau);
422  // If res lies close enough, within the desired tolerance, calculate the
423  // exact residual
424  if (res < additional_data.solver_tolerance)
425  {
426  A.vmult(u, x);
427  u.sadd(-1., 1., b);
428  res = u.l2_norm();
429  }
430  state = this->iteration_status(step, res, x);
431  if ((state == SolverControl::success) ||
432  (state == SolverControl::failure))
433  return IterationResult(state, res);
434 
435  //--------------------------------------------------------------
436  // Step 3: check breakdown criterion and update the vectors
437  //--------------------------------------------------------------
438  if (additional_data.breakdown_testing == true &&
439  std::fabs(sigma) < additional_data.breakdown_threshold)
440  return IterationResult(SolverControl::iterate, res);
441 
442  const double rho_old = rho;
443 
444  // Applying the preconditioner
445  if (additional_data.left_preconditioning)
446  {
447  // Left preconditioning
448  u = t;
449  }
450  else
451  {
452  // Right preconditioning
453  preconditioner.vmult(u, t);
454  }
455 
456  // Double and vector updates
457  rho = u * r;
458  const double beta = rho / rho_old;
459  q.sadd(beta, 1., u);
460  }
461  return IterationResult(SolverControl::success, res);
462 }
463 
464 #endif // DOXYGEN
465 
467 
468 #endif
@ iterate
Continue iteration.
@ success
Stop iteration, goal reached.
@ failure
Stop iteration, goal not reached.
virtual void print_vectors(const unsigned int step, const VectorType &x, const VectorType &r, const VectorType &d) const
IterationResult iterate(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner, VectorType &r, VectorType &u, VectorType &q, VectorType &t, VectorType &d)
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
unsigned int step
Definition: solver_qmrs.h:234
SolverQMRS(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
AdditionalData additional_data
Definition: solver_qmrs.h:199
SolverQMRS(SolverControl &cn, const AdditionalData &data=AdditionalData())
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:475
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:476
#define AssertThrow(cond, exc)
Definition: exceptions.h:1703
LogStream deallog
Definition: logstream.cc:37
Expression fabs(const Expression &x)
static const char A
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
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:129
SolverControl::State state
Definition: solver_qmrs.h:208
IterationResult(const SolverControl::State state, const double last_residual)