Reference documentation for deal.II version 9.5.0
\(\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// Copyright (C) 1999 - 2022 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
23
24#include <deal.II/lac/solver.h>
26
27#include <cmath>
28
30
93template <typename VectorType = Vector<double>>
94class SolverQMRS : public SolverBase<VectorType>
95{
96public:
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
195protected:
200
201private:
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
243template <class VectorType>
245 const SolverControl::State state,
246 const double last_residual)
247 : state(state)
248 , last_residual(last_residual)
249{}
250
251
252
253template <class VectorType>
256 const AdditionalData & data)
257 : SolverBase<VectorType>(cn, mem)
258 , additional_data(data)
259 , step(0)
260{}
261
262template <class VectorType>
264 const AdditionalData &data)
265 : SolverBase<VectorType>(cn)
266 , additional_data(data)
267 , step(0)
268{}
269
270template <class VectorType>
271void
273 const VectorType &,
274 const VectorType &,
275 const VectorType &) const
276{}
277
278template <class VectorType>
279template <typename MatrixType, typename PreconditionerType>
280void
281SolverQMRS<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
326template <class VectorType>
327template <typename MatrixType, typename PreconditionerType>
329SolverQMRS<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
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:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
#define AssertThrow(cond, exc)
LogStream deallog
Definition logstream.cc:37
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)