LinearSolver::Serial::Gmres< VarsType, Preconditioner, LinearSolverRegistrars > Class Template Reference

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

#include <Gmres.hpp>

Public Types
using	options = implementation defined

Public Member Functions
	Gmres (Convergence::Criteria convergence_criteria, ::Verbosity verbosity, std::optional< size_t > restart=std::nullopt, std::optional< typename Base::PreconditionerType > local_preconditioner=std::nullopt, const Options::Context &context={})
	Gmres (Gmres &&)=default
Gmres &	operator= (Gmres &&)=default
	Gmres (const Gmres &rhs)
Gmres &	operator= (const Gmres &rhs)
std::unique_ptr< LinearSolver< LinearSolverRegistrars > >	get_clone () const override
void	initialize ()
const Convergence::Criteria &	convergence_criteria () const
::Verbosity	verbosity () const
size_t	restart () const
void	pup (PUP::er &p) override
template<typename LinearOperator , typename SourceType , typename... OperatorArgs, typename IterationCallback = NoIterationCallback>
Convergence::HasConverged	solve (gsl::not_null< VarsType * > initial_guess_in_solution_out, const LinearOperator &linear_operator, const SourceType &source, const std::tuple< OperatorArgs... > &operator_args=std::tuple{}, const IterationCallback &iteration_callback=NoIterationCallback{}) const
void	reset () override
template<typename LinearOperator , typename SourceType , typename... OperatorArgs, typename IterationCallback >
Convergence::HasConverged	solve (const gsl::not_null< VarsType * > initial_guess_in_solution_out, const LinearOperator &linear_operator, const SourceType &source, const std::tuple< OperatorArgs... > &operator_args, const IterationCallback &iteration_callback) const

Static Public Attributes
static constexpr Options::String	help

Detailed Description

template<typename VarsType, typename Preconditioner = NoPreconditioner, typename LinearSolverRegistrars = tmpl::list<Registrars::Gmres<VarsType>>>
class LinearSolver::Serial::Gmres< VarsType, Preconditioner, LinearSolverRegistrars >

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 [174], chapter 6.5 for a description of the GMRES algorithm, Algorithm 9.6 for this implementation, and Sec. 6.5.9 for the extension to complex arithmetic. The GMRES method is matrix-free, which means the operator $A$ needs not be provided explicitly 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 [174], 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 [174], 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{r}\mathrm{e}\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{t}}$ 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{r}\mathrm{e}\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{t}}$ to activate restarting, or set it to 'None' to deactivate restarting.

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 [174], sections 9.3.2 and 9.4.1 for details. This implementation follows Algorithm 9.6 in [174].

Improvements:

Further improvements can potentially be implemented for this algorithm, see e.g. [7].

Example

    INFO("Solve a symmetric 2x2 matrix");
    blaze::DynamicMatrix<double> matrix{{4., 1.}, {1., 3.}};
    const helpers::ApplyMatrix<double> linear_operator{std::move(matrix)};
    const blaze::DynamicVector<double> source{1., 2.};
    blaze::DynamicVector<double> initial_guess_in_solution_out{2., 1.};
    const blaze::DynamicVector<double> expected_solution{0.0909090909090909,
                                                         0.6363636363636364};
    const Convergence::Criteria convergence_criteria{2, 1.e-14, 0.};
    const Gmres<blaze::DynamicVector<double>> gmres{convergence_criteria,
                                                    ::Verbosity::Verbose};
    CHECK_FALSE(gmres.has_preconditioner());
    std::vector<double> recorded_residuals;
    const auto has_converged =
        gmres.solve(make_not_null(&initial_guess_in_solution_out),
                    linear_operator, source, std::tuple{},
                    [&recorded_residuals](
                        const Convergence::HasConverged& local_has_converged) {
                      recorded_residuals.push_back(
                          local_has_converged.residual_magnitude());
                    });
    REQUIRE(has_converged);
    CHECK(linear_operator.invocations == 3);
    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(initial_guess_in_solution_out, 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]);
    }

Template Parameters

VarsType

The type of the vectors x and b in the linear equation $Ax=b$. This type must support addition, subtraction, scalar multiplication, and an inner product. It must expose an ElementType alias that is the type of the components of the vector, e.g. double or std::complex<double>.

Member Data Documentation

help

template<typename VarsType , typename Preconditioner = NoPreconditioner, typename LinearSolverRegistrars = tmpl::list<Registrars::Gmres<VarsType>>>

constexpr Options::String LinearSolver::Serial::Gmres< VarsType, Preconditioner, LinearSolverRegistrars >::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"
      "Preconditioning: Specify a preconditioner to run in every GMRES "
      "iteration to accelerate the solve, or 'None' to disable "
      "preconditioning. The choice of preconditioner can be crucial to obtain "
      "good convergence.\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 'None'."

---

The documentation for this class was generated from the following file:

• src/NumericalAlgorithms/LinearSolver/Gmres.hpp