Reference documentation for deal.II version GIT relicensing-1062-gc06da148b8 2024-07-15 19:20: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\}}\)
Loading...
Searching...
No Matches
solver_qmrs.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 1999 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_solver_qmrs_h
16#define dealii_solver_qmrs_h
17
18#include <deal.II/base/config.h>
19
23
24#include <deal.II/lac/solver.h>
26
27#include <cmath>
28
30
93template <typename VectorType = Vector<double>>
95class SolverQMRS : public SolverBase<VectorType>
96{
97public:
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)
134 : left_preconditioning(left_preconditioning)
135 , solver_tolerance(solver_tolerance)
136 , breakdown_testing(breakdown_testing)
137 , breakdown_threshold(breakdown_threshold)
138 {}
139
144
149
154
160 };
161
167 const AdditionalData &data = AdditionalData());
168
174
178 template <typename MatrixType, typename PreconditionerType>
182 void solve(const MatrixType &A,
183 VectorType &x,
184 const VectorType &b,
185 const PreconditionerType &preconditioner);
186
192 virtual void
193 print_vectors(const unsigned int step,
194 const VectorType &x,
195 const VectorType &r,
196 const VectorType &d) const;
197
198protected:
202 AdditionalData additional_data;
203
204private:
210 {
213
215 const double last_residual);
216 };
217
222 template <typename MatrixType, typename PreconditionerType>
224 iterate(const MatrixType &A,
225 VectorType &x,
226 const VectorType &b,
227 const PreconditionerType &preconditioner,
228 VectorType &r,
229 VectorType &u,
230 VectorType &q,
231 VectorType &t,
232 VectorType &d);
233
237 unsigned int step;
238};
239
241/*------------------------- Implementation ----------------------------*/
242
243#ifndef DOXYGEN
244
245
246template <typename VectorType>
249 const SolverControl::State state,
250 const double last_residual)
251 : state(state)
252 , last_residual(last_residual)
253{}
254
255
256
257template <typename VectorType>
261 const AdditionalData &data)
262 : SolverBase<VectorType>(cn, mem)
263 , additional_data(data)
264 , step(0)
265{}
266
267
268
269template <typename VectorType>
272 const AdditionalData &data)
274 , additional_data(data)
275 , step(0)
276{}
277
278
279
280template <typename VectorType>
282void SolverQMRS<VectorType>::print_vectors(const unsigned int,
283 const VectorType &,
284 const VectorType &,
285 const VectorType &) const
286{}
287
288
289
290template <typename VectorType>
292template <typename MatrixType, typename PreconditionerType>
296void SolverQMRS<VectorType>::solve(const MatrixType &A,
297 VectorType &x,
298 const VectorType &b,
299 const PreconditionerType &preconditioner)
300{
301 LogStream::Prefix prefix("SQMR");
302
303
304 // temporary vectors, allocated through the @p VectorMemory object at the
305 // start of the actual solution process and deallocated at the end.
306 typename VectorMemory<VectorType>::Pointer Vr(this->memory);
307 typename VectorMemory<VectorType>::Pointer Vu(this->memory);
308 typename VectorMemory<VectorType>::Pointer Vq(this->memory);
309 typename VectorMemory<VectorType>::Pointer Vt(this->memory);
310 typename VectorMemory<VectorType>::Pointer Vd(this->memory);
311
312
313 // resize the vectors, but do not set
314 // the values since they'd be overwritten
315 // soon anyway.
316 Vr->reinit(x, true);
317 Vu->reinit(x, true);
318 Vq->reinit(x, true);
319 Vt->reinit(x, true);
320 Vd->reinit(x, true);
321
322 step = 0;
323
324 IterationResult state(SolverControl::failure, 0);
325
326 do
327 {
328 if (step > 0)
329 deallog << "Restart step " << step << std::endl;
330 state = iterate(A, x, b, preconditioner, *Vr, *Vu, *Vq, *Vt, *Vd);
331 }
332 while (state.state == SolverControl::iterate);
333
334
335 // in case of failure: throw exception
336 AssertThrow(state.state == SolverControl::success,
337 SolverControl::NoConvergence(step, state.last_residual));
338 // otherwise exit as normal
339}
340
341
342
343template <typename VectorType>
345template <typename MatrixType, typename PreconditionerType>
347 SolverQMRS<VectorType>::iterate(const MatrixType &A,
348 VectorType &x,
349 const VectorType &b,
350 const PreconditionerType &preconditioner,
351 VectorType &r,
352 VectorType &u,
353 VectorType &q,
354 VectorType &t,
355 VectorType &d)
356{
358
359 int it = 0;
360
361 double tau, rho, theta = 0;
362 double res;
363
364 // Compute the start residual
365 A.vmult(r, x);
366 r.sadd(-1., 1., b);
367
368 // Doing the initial preconditioning
369 if (additional_data.left_preconditioning)
370 {
371 // Left preconditioning
372 preconditioner.vmult(t, r);
373 q = t;
374 }
375 else
376 {
377 // Right preconditioning
378 t = r;
379 preconditioner.vmult(q, t);
380 }
381
382 tau = t.norm_sqr();
383 res = std::sqrt(tau);
384
385 if (this->iteration_status(step, res, x) == SolverControl::success)
386 return IterationResult(SolverControl::success, res);
387
388 rho = q * r;
389
390 while (state == SolverControl::iterate)
391 {
392 ++step;
393 ++it;
394 //--------------------------------------------------------------
395 // Step 1: apply the system matrix and compute one inner product
396 //--------------------------------------------------------------
397 A.vmult(t, q);
398 const double sigma = q * t;
399
400 // Check the breakdown criterion
401 if (additional_data.breakdown_testing == true &&
402 std::fabs(sigma) < additional_data.breakdown_threshold)
403 return IterationResult(SolverControl::iterate, res);
404 // Update the residual
405 const double alpha = rho / sigma;
406 r.add(-alpha, t);
407
408 //--------------------------------------------------------------
409 // Step 2: update the solution vector
410 //--------------------------------------------------------------
411 const double theta_old = theta;
412
413 // Apply the preconditioner
414 if (additional_data.left_preconditioning)
415 {
416 // Left Preconditioning
417 preconditioner.vmult(t, r);
418 }
419 else
420 {
421 // Right Preconditioning
422 t = r;
423 }
424
425 // Double updates
426 theta = t * t / tau;
427 const double psi = 1. / (1. + theta);
428 tau *= theta * psi;
429
430 // Actual update of the solution vector
431 d.sadd(psi * theta_old, psi * alpha, q);
432 x += d;
433
434 print_vectors(step, x, r, d);
435
436 // Check for convergence
437 // Compute a simple and cheap upper bound of the norm of the residual
438 // vector b-Ax
439 res = std::sqrt((it + 1) * tau);
440 // If res lies close enough, within the desired tolerance, calculate the
441 // exact residual
442 if (res < additional_data.solver_tolerance)
443 {
444 A.vmult(u, x);
445 u.sadd(-1., 1., b);
446 res = u.l2_norm();
447 }
448 state = this->iteration_status(step, res, x);
449 if ((state == SolverControl::success) ||
450 (state == SolverControl::failure))
451 return IterationResult(state, res);
452
453 //--------------------------------------------------------------
454 // Step 3: check breakdown criterion and update the vectors
455 //--------------------------------------------------------------
456 if (additional_data.breakdown_testing == true &&
457 std::fabs(sigma) < additional_data.breakdown_threshold)
458 return IterationResult(SolverControl::iterate, res);
459
460 const double rho_old = rho;
461
462 // Applying the preconditioner
463 if (additional_data.left_preconditioning)
464 {
465 // Left preconditioning
466 u = t;
467 }
468 else
469 {
470 // Right preconditioning
471 preconditioner.vmult(u, t);
472 }
473
474 // Double and vector updates
475 rho = u * r;
476 const double beta = rho / rho_old;
477 q.sadd(beta, 1., u);
478 }
479 return IterationResult(SolverControl::success, res);
480}
481
482#endif // DOXYGEN
483
485
486#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
SolverQMRS(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
SolverQMRS(SolverControl &cn, const AdditionalData &data=AdditionalData())
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_CXX20_REQUIRES(condition)
Definition config.h:177
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
#define AssertThrow(cond, exc)
LogStream deallog
Definition logstream.cc:36
Tpetra::Vector< Number, LO, GO, NodeType< MemorySpace > > VectorType
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
::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)
SolverControl::State state
IterationResult(const SolverControl::State state, const double last_residual)