LinearSolver::Serial::Gmres< VarsType > Struct Template Reference

A serial GMRES iterative solver for nonsymmetric linear systems of equations. More...

#include <Gmres.hpp>

## Public Types

using options = tmpl::list< ConvergenceCriteria, Verbosity, Restart >

## Public Member Functions

Gmres (Convergence::Criteria convergence_criteria, ::Verbosity verbosity, size_t restart=0) noexcept

Gmres (const Gmres &)=default

Gmresoperator= (const Gmres &)=default

Gmres (Gmres &&)=default

Gmresoperator= (Gmres &&)=default

void initialize () noexcept

const Convergence::Criteriaconvergence_criteria () const noexcept

void pup (PUP::er &p) noexcept

template<typename LinearOperator , typename SourceType , typename Preconditioner = IdentityPreconditioner<VarsType>, typename IterationCallback = NoIterationCallback>
std::pair< Convergence::HasConverged, VarsType > operator() (LinearOperator &&linear_operator, const SourceType &source, const VarsType &initial_guess, Preconditioner &&preconditioner=IdentityPreconditioner< VarsType >{}, IterationCallback &&iteration_callback=NoIterationCallback{}) const noexcept
Iteratively solve the problem $$Ax=b$$ for $$x$$ where $$A$$ is the linear_operator and $$b$$ is the source, starting $$x$$ at initial_guess. More...

## Static Public Attributes

static constexpr OptionString help

## Detailed Description

### template<typename VarsType> struct LinearSolver::Serial::Gmres< VarsType >

A serial GMRES iterative solver for nonsymmetric linear systems of equations.

This is an iterative algorithm to solve general linear equations $$Ax=b$$ where $$A$$ is a linear operator. See [55], chapter 6.5 for a description of the GMRES algorithm and Algorithm 9.6 for this implementation. It is matrix-free, which means the operator $$A$$ needs not be provided explicity as a matrix but only the operator action $$A(x)$$ must be provided for an argument $$x$$.

The GMRES algorithm does not require the operator $$A$$ to be symmetric or positive-definite. Note that other algorithms such as conjugate gradients may be more efficient for symmetric positive-definite operators.

Convergence:
Given a set of $$N_A$$ equations (e.g. through an $$N_A\times N_A$$ matrix) the GMRES algorithm will converge to numerical precision in at most $$N_A$$ iterations. However, depending on the properties of the linear operator, an approximate solution can ideally be obtained in only a few iterations. See [55], section 6.11.4 for details on the convergence of the GMRES algorithm.
Restarting:
This implementation of the GMRES algorithm supports restarting, as detailed in [55], section 6.5.5. Since the GMRES algorithm iteratively builds up an orthogonal basis of the solution space the cost of each iteration increases linearly with the number of iterations. Therefore it is sometimes helpful to restart the algorithm every $$N_\mathrm{restart}$$ iterations, discarding the set of basis vectors and starting again from the current solution estimate. This strategy can improve the performance of the solver, but note that the solver can stagnate for non-positive-definite operators and is not guaranteed to converge within $$N_A$$ iterations anymore. Set the restart argument of the constructor to $$N_\mathrm{restart}$$ to activate restarting, or set it to zero to deactivate restarting (default behaviour).
Preconditioning:
This implementation of the GMRES algorithm also supports preconditioning. You can provide a linear operator $$P$$ that approximates the inverse of the operator $$A$$ to accelerate the convergence of the linear solve. The algorithm is right-preconditioned, which allows the preconditioner to change in every iteration ("flexible" variant). See [55], sections 9.3.2 and 9.4.1 for details. This implementation follows Algorithm 9.6 in [55].
Improvements:
Further improvements can potentially be implemented for this algorithm, see e.g. [3].

### Example

INFO("Solve a symmetric 2x2 matrix");
DenseMatrix<double> matrix(2, 2);
matrix(0, 0) = 4.;
matrix(0, 1) = 1.;
matrix(1, 0) = 1.;
matrix(1, 1) = 3.;
const DenseVector<double> source{1., 2.};
const DenseVector<double> initial_guess{2., 1.};
const DenseVector<double> expected_solution{0.0909090909090909,
0.6363636363636364};
const Convergence::Criteria convergence_criteria{2, 1.e-14, 0.};
const Gmres<DenseVector<double>> gmres{convergence_criteria,
::Verbosity::Verbose};
const auto linear_operator =
[&matrix](const DenseVector<double>& arg) noexcept {
return matrix * arg;
};
std::vector<double> recorded_residuals;
const auto result = gmres(
linear_operator, source, initial_guess,
IdentityPreconditioner<DenseVector<double>>{},
[&recorded_residuals](const Convergence::HasConverged& has_converged) {
recorded_residuals.push_back(has_converged.residual_magnitude());
});
const auto& has_converged = result.first;
const auto& solution = result.second;
REQUIRE(has_converged);
CHECK(has_converged.reason() == Convergence::Reason::AbsoluteResidual);
CHECK(has_converged.num_iterations() == 2);
CHECK(has_converged.residual_magnitude() <= 1.e-14);
CHECK(has_converged.initial_residual_magnitude() ==
approx(8.54400374531753));
CHECK_ITERABLE_APPROX(solution, expected_solution);
// The residuals should decrease monotonically
CHECK(recorded_residuals[0] == has_converged.initial_residual_magnitude());
for (size_t i = 1; i < has_converged.num_iterations(); ++i) {
CHECK(recorded_residuals[i] <= recorded_residuals[i - 1]);
}

## ◆ operator()()

template<typename VarsType >
template<typename LinearOperator , typename SourceType , typename Preconditioner , typename IterationCallback >
 std::pair< Convergence::HasConverged, VarsType > LinearSolver::Serial::Gmres< VarsType >::operator() ( LinearOperator && linear_operator, const SourceType & source, const VarsType & initial_guess, Preconditioner && preconditioner = IdentityPreconditioner{}, IterationCallback && iteration_callback = NoIterationCallback{} ) const
noexcept

Iteratively solve the problem $$Ax=b$$ for $$x$$ where $$A$$ is the linear_operator and $$b$$ is the source, starting $$x$$ at initial_guess.

Optionally provide a preconditioner (see class documentation).

Returns: An instance of Convergence::HasConverged that provides information on the convergence status of the completed solve, and the approximate solution $$x$$.

## ◆ help

template<typename VarsType >
 constexpr OptionString LinearSolver::Serial::Gmres< VarsType >::help
staticconstexpr
Initial value:
=
"A serial GMRES iterative solver for nonsymmetric linear systems of\n"
"equations Ax=b. It will converge to numerical precision in at most N_A\n"
"iterations, where N_A is the number of equations represented by the\n"
"linear operator A, but will ideally converge to a reasonable\n"
"approximation of the solution x in only a few iterations.\n"
"\n"
"Restarting: It is sometimes helpful to restart the algorithm every\n"
"N_restart iterations to speed it up. Note that it can stagnate for\n"
"non-positive-definite matrices and is not guaranteed to converge\n"
"within N_A iterations anymore when restarting is activated.\n"
"Activate restarting by setting the 'Restart' option to N_restart, or\n"
"deactivate restarting by setting it to zero (default)."

The documentation for this struct was generated from the following file:
• src/NumericalAlgorithms/LinearSolver/Gmres.hpp
std::vector< double >
Convergence::Criteria
Criteria that determine an iterative algorithm has converged.
Definition: Criteria.hpp:35
DenseVector< double >
Convergence::HasConverged
Signals convergence of the algorithm.
Definition: HasConverged.hpp:59
CHECK_ITERABLE_APPROX
#define CHECK_ITERABLE_APPROX(a, b)
A wrapper around Catch's CHECK macro that checks approximate equality of entries in iterable containe...
Definition: TestingFramework.hpp:139
DenseMatrix< double >