Reference documentation for deal.II version GIT relicensing-478-g3275795167 2024-04-24 07:10: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
22
23#include <deal.II/lac/solver.h>
25
26#include <cmath>
27
29
92template <typename VectorType = Vector<double>>
93class SolverQMRS : public SolverBase<VectorType>
94{
95public:
159
165 const AdditionalData &data = AdditionalData());
166
172
176 template <typename MatrixType, typename PreconditionerType>
177 void
178 solve(const MatrixType &A,
179 VectorType &x,
180 const VectorType &b,
181 const PreconditionerType &preconditioner);
182
188 virtual void
189 print_vectors(const unsigned int step,
190 const VectorType &x,
191 const VectorType &r,
192 const VectorType &d) const;
193
194protected:
199
200private:
213
218 template <typename MatrixType, typename PreconditionerType>
220 iterate(const MatrixType &A,
221 VectorType &x,
222 const VectorType &b,
223 const PreconditionerType &preconditioner,
224 VectorType &r,
225 VectorType &u,
226 VectorType &q,
227 VectorType &t,
228 VectorType &d);
229
233 unsigned int step;
234};
235
237/*------------------------- Implementation ----------------------------*/
238
239#ifndef DOXYGEN
240
241
242template <typename VectorType>
244 const SolverControl::State state,
245 const double last_residual)
246 : state(state)
247 , last_residual(last_residual)
248{}
249
250
251
252template <typename VectorType>
255 const AdditionalData &data)
256 : SolverBase<VectorType>(cn, mem)
257 , additional_data(data)
258 , step(0)
259{}
260
261template <typename VectorType>
263 const AdditionalData &data)
264 : SolverBase<VectorType>(cn)
265 , additional_data(data)
266 , step(0)
267{}
268
269template <typename VectorType>
270void
272 const VectorType &,
273 const VectorType &,
274 const VectorType &) const
275{}
276
277template <typename VectorType>
278template <typename MatrixType, typename PreconditionerType>
279void
280SolverQMRS<VectorType>::solve(const MatrixType &A,
281 VectorType &x,
282 const VectorType &b,
283 const PreconditionerType &preconditioner)
284{
285 LogStream::Prefix prefix("SQMR");
286
287
288 // temporary vectors, allocated through the @p VectorMemory object at the
289 // start of the actual solution process and deallocated at the end.
290 typename VectorMemory<VectorType>::Pointer Vr(this->memory);
291 typename VectorMemory<VectorType>::Pointer Vu(this->memory);
292 typename VectorMemory<VectorType>::Pointer Vq(this->memory);
293 typename VectorMemory<VectorType>::Pointer Vt(this->memory);
294 typename VectorMemory<VectorType>::Pointer Vd(this->memory);
295
296
297 // resize the vectors, but do not set
298 // the values since they'd be overwritten
299 // soon anyway.
300 Vr->reinit(x, true);
301 Vu->reinit(x, true);
302 Vq->reinit(x, true);
303 Vt->reinit(x, true);
304 Vd->reinit(x, true);
305
306 step = 0;
307
308 IterationResult state(SolverControl::failure, 0);
309
310 do
311 {
312 if (step > 0)
313 deallog << "Restart step " << step << std::endl;
314 state = iterate(A, x, b, preconditioner, *Vr, *Vu, *Vq, *Vt, *Vd);
315 }
316 while (state.state == SolverControl::iterate);
317
318
319 // in case of failure: throw exception
320 AssertThrow(state.state == SolverControl::success,
321 SolverControl::NoConvergence(step, state.last_residual));
322 // otherwise exit as normal
323}
324
325template <typename VectorType>
326template <typename MatrixType, typename PreconditionerType>
328SolverQMRS<VectorType>::iterate(const MatrixType &A,
329 VectorType &x,
330 const VectorType &b,
331 const PreconditionerType &preconditioner,
332 VectorType &r,
333 VectorType &u,
334 VectorType &q,
335 VectorType &t,
336 VectorType &d)
337{
339
340 int it = 0;
341
342 double tau, rho, theta = 0;
343 double res;
344
345 // Compute the start residual
346 A.vmult(r, x);
347 r.sadd(-1., 1., b);
348
349 // Doing the initial preconditioning
350 if (additional_data.left_preconditioning)
351 {
352 // Left preconditioning
353 preconditioner.vmult(t, r);
354 q = t;
355 }
356 else
357 {
358 // Right preconditioning
359 t = r;
360 preconditioner.vmult(q, t);
361 }
362
363 tau = t.norm_sqr();
364 res = std::sqrt(tau);
365
366 if (this->iteration_status(step, res, x) == SolverControl::success)
367 return IterationResult(SolverControl::success, res);
368
369 rho = q * r;
370
371 while (state == SolverControl::iterate)
372 {
373 ++step;
374 ++it;
375 //--------------------------------------------------------------
376 // Step 1: apply the system matrix and compute one inner product
377 //--------------------------------------------------------------
378 A.vmult(t, q);
379 const double sigma = q * t;
380
381 // Check the breakdown criterion
382 if (additional_data.breakdown_testing == true &&
383 std::fabs(sigma) < additional_data.breakdown_threshold)
384 return IterationResult(SolverControl::iterate, res);
385 // Update the residual
386 const double alpha = rho / sigma;
387 r.add(-alpha, t);
388
389 //--------------------------------------------------------------
390 // Step 2: update the solution vector
391 //--------------------------------------------------------------
392 const double theta_old = theta;
393
394 // Apply the preconditioner
395 if (additional_data.left_preconditioning)
396 {
397 // Left Preconditioning
398 preconditioner.vmult(t, r);
399 }
400 else
401 {
402 // Right Preconditioning
403 t = r;
404 }
405
406 // Double updates
407 theta = t * t / tau;
408 const double psi = 1. / (1. + theta);
409 tau *= theta * psi;
410
411 // Actual update of the solution vector
412 d.sadd(psi * theta_old, psi * alpha, q);
413 x += d;
414
415 print_vectors(step, x, r, d);
416
417 // Check for convergence
418 // Compute a simple and cheap upper bound of the norm of the residual
419 // vector b-Ax
420 res = std::sqrt((it + 1) * tau);
421 // If res lies close enough, within the desired tolerance, calculate the
422 // exact residual
423 if (res < additional_data.solver_tolerance)
424 {
425 A.vmult(u, x);
426 u.sadd(-1., 1., b);
427 res = u.l2_norm();
428 }
429 state = this->iteration_status(step, res, x);
430 if ((state == SolverControl::success) ||
431 (state == SolverControl::failure))
432 return IterationResult(state, res);
433
434 //--------------------------------------------------------------
435 // Step 3: check breakdown criterion and update the vectors
436 //--------------------------------------------------------------
437 if (additional_data.breakdown_testing == true &&
438 std::fabs(sigma) < additional_data.breakdown_threshold)
439 return IterationResult(SolverControl::iterate, res);
440
441 const double rho_old = rho;
442
443 // Applying the preconditioner
444 if (additional_data.left_preconditioning)
445 {
446 // Left preconditioning
447 u = t;
448 }
449 else
450 {
451 // Right preconditioning
452 preconditioner.vmult(u, t);
453 }
454
455 // Double and vector updates
456 rho = u * r;
457 const double beta = rho / rho_old;
458 q.sadd(beta, 1., u);
459 }
460 return IterationResult(SolverControl::success, res);
461}
462
463#endif // DOXYGEN
464
466
467#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())
AdditionalData additional_data
SolverQMRS(SolverControl &cn, const AdditionalData &data=AdditionalData())
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
#define AssertThrow(cond, exc)
LogStream deallog
Definition logstream.cc:36
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)