Reference documentation for deal.II version Git 193422c69f 2020-07-08 17:07:46 +0200
\(\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 - 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>
23 
24 #include <deal.II/lac/solver.h>
26 
27 #include <cmath>
28 
30 
33 
91 template <typename VectorType = Vector<double>>
92 class SolverQMRS : public SolverBase<VectorType>
93 {
94 public:
120  {
127  explicit AdditionalData(const bool left_preconditioning = false,
128  const double solver_tolerance = 1.e-9,
129  const bool breakdown_testing = true,
130  const double breakdown_threshold = 1.e-16)
135  {}
136 
141 
146 
151 
157  };
158 
164  const AdditionalData & data = AdditionalData());
165 
171 
175  template <typename MatrixType, typename PreconditionerType>
176  void
177  solve(const MatrixType & A,
178  VectorType & x,
179  const VectorType & b,
180  const PreconditionerType &preconditioner);
181 
187  virtual void
188  print_vectors(const unsigned int step,
189  const VectorType & x,
190  const VectorType & r,
191  const VectorType & d) const;
192 
193 protected:
198 
199 private:
205  {
208 
210  const double last_residual);
211  };
212 
217  template <typename MatrixType, typename PreconditionerType>
219  iterate(const MatrixType & A,
220  VectorType & x,
221  const VectorType & b,
222  const PreconditionerType &preconditioner,
223  VectorType & r,
224  VectorType & u,
225  VectorType & q,
226  VectorType & t,
227  VectorType & d);
228 
232  unsigned int step;
233 };
234 
236 /*------------------------- Implementation ----------------------------*/
237 
238 #ifndef DOXYGEN
239 
240 
241 template <class VectorType>
243  const SolverControl::State state,
244  const double last_residual)
245  : state(state)
246  , last_residual(last_residual)
247 {}
248 
249 
250 
251 template <class VectorType>
254  const AdditionalData & data)
255  : SolverBase<VectorType>(cn, mem)
256  , additional_data(data)
257  , step(0)
258 {}
259 
260 template <class VectorType>
262  const AdditionalData &data)
264  , additional_data(data)
265  , step(0)
266 {}
267 
268 template <class VectorType>
269 void
270 SolverQMRS<VectorType>::print_vectors(const unsigned int,
271  const VectorType &,
272  const VectorType &,
273  const VectorType &) const
274 {}
275 
276 template <class VectorType>
277 template <typename MatrixType, typename PreconditionerType>
278 void
279 SolverQMRS<VectorType>::solve(const MatrixType & A,
280  VectorType & x,
281  const VectorType & b,
282  const PreconditionerType &preconditioner)
283 {
284  LogStream::Prefix prefix("SQMR");
285 
286 
287  // temporary vectors, allocated trough the @p VectorMemory object at the
288  // start of the actual solution process and deallocated at the end.
289  typename VectorMemory<VectorType>::Pointer Vr(this->memory);
290  typename VectorMemory<VectorType>::Pointer Vu(this->memory);
291  typename VectorMemory<VectorType>::Pointer Vq(this->memory);
292  typename VectorMemory<VectorType>::Pointer Vt(this->memory);
293  typename VectorMemory<VectorType>::Pointer Vd(this->memory);
294 
295 
296  // resize the vectors, but do not set
297  // the values since they'd be overwritten
298  // soon anyway.
299  Vr->reinit(x, true);
300  Vu->reinit(x, true);
301  Vq->reinit(x, true);
302  Vt->reinit(x, true);
303  Vd->reinit(x, true);
304 
305  step = 0;
306 
308 
309  do
310  {
311  if (step > 0)
312  deallog << "Restart step " << step << std::endl;
313  state = iterate(A, x, b, preconditioner, *Vr, *Vu, *Vq, *Vt, *Vd);
314  }
315  while (state.state == SolverControl::iterate);
316 
317 
318  // in case of failure: throw exception
321  // otherwise exit as normal
322 }
323 
324 template <class VectorType>
325 template <typename MatrixType, typename PreconditionerType>
327 SolverQMRS<VectorType>::iterate(const MatrixType & A,
328  VectorType & x,
329  const VectorType & b,
330  const PreconditionerType &preconditioner,
331  VectorType & r,
332  VectorType & u,
333  VectorType & q,
334  VectorType & t,
335  VectorType & d)
336 {
338 
339  int it = 0;
340 
341  double tau, rho, theta = 0;
342  double res;
343 
344  // Compute the start residual
345  A.vmult(r, x);
346  r.sadd(-1., 1., b);
347 
348  // Doing the initial preconditioning
350  {
351  // Left preconditioning
352  preconditioner.vmult(t, r);
353  q = t;
354  }
355  else
356  {
357  // Right preconditioning
358  t = r;
359  preconditioner.vmult(q, t);
360  }
361 
362  tau = t.norm_sqr();
363  res = std::sqrt(tau);
364 
365  if (this->iteration_status(step, res, x) == SolverControl::success)
367 
368  rho = q * r;
369 
370  while (state == SolverControl::iterate)
371  {
372  step++;
373  it++;
374  //--------------------------------------------------------------
375  // Step 1: apply the system matrix and compute one inner product
376  //--------------------------------------------------------------
377  A.vmult(t, q);
378  const double sigma = q * t;
379 
380  // Check the breakdown criterion
381  if (additional_data.breakdown_testing == true &&
384  // Update the residual
385  const double alpha = rho / sigma;
386  r.add(-alpha, t);
387 
388  //--------------------------------------------------------------
389  // Step 2: update the solution vector
390  //--------------------------------------------------------------
391  const double theta_old = theta;
392 
393  // Apply the preconditioner
395  {
396  // Left Preconditioning
397  preconditioner.vmult(t, r);
398  }
399  else
400  {
401  // Right Preconditioning
402  t = r;
403  }
404 
405  // Double updates
406  theta = t * t / tau;
407  const double psi = 1. / (1. + theta);
408  tau *= theta * psi;
409 
410  // Actual update of the solution vector
411  d.sadd(psi * theta_old, psi * alpha, q);
412  x += d;
413 
414  print_vectors(step, x, r, d);
415 
416  // Check for convergence
417  // Compute a simple and cheap upper bound of the norm of the residual
418  // vector b-Ax
419  res = std::sqrt((it + 1) * tau);
420  // If res lies close enough, within the desired tolerance, calculate the
421  // exact residual
423  {
424  A.vmult(u, x);
425  u.sadd(-1., 1., b);
426  res = u.l2_norm();
427  }
428  state = this->iteration_status(step, res, x);
429  if ((state == SolverControl::success) ||
430  (state == SolverControl::failure))
431  return IterationResult(state, res);
432 
433  //--------------------------------------------------------------
434  // Step 3: check breakdown criterion and update the vectors
435  //--------------------------------------------------------------
436  if (additional_data.breakdown_testing == true &&
439 
440  const double rho_old = rho;
441 
442  // Applying the preconditioner
444  {
445  // Left preconditioning
446  u = t;
447  }
448  else
449  {
450  // Right preconditioning
451  preconditioner.vmult(u, t);
452  }
453 
454  // Double and vector updates
455  rho = u * r;
456  const double beta = rho / rho_old;
457  q.sadd(beta, 1., u);
458  }
460 }
461 
462 #endif // DOXYGEN
463 
465 
466 #endif
Stop iteration, goal not reached.
LogStream deallog
Definition: logstream.cc:37
Continue iteration.
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
boost::signals2::signal< SolverControl::State(const unsigned int iteration, const double check_value, const VectorType &current_iterate), StateCombiner > iteration_status
Definition: solver.h:471
#define AssertThrow(cond, exc)
Definition: exceptions.h:1513
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:232
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
Expression fabs(const Expression &x)
SolverControl::State state
Definition: solver_qmrs.h:206
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)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
static const char A
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:127
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:362
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:197
IterationResult(const SolverControl::State state, const double last_residual)
VectorMemory< VectorType > & memory
Definition: solver.h:420