SpECTRE Documentation Coverage Report
Current view: top level - NumericalAlgorithms/LinearSolver - Gmres.hpp Hit Total Coverage
Commit: ebec864322c50bab8dca0a90baf8d01875114261 Lines: 2 39 5.1 %
Date: 2020-11-25 20:28:50
Legend: Lines: hit not hit

          Line data    Source code
       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <boost/none.hpp>
       7             : #include <boost/optional.hpp>
       8             : #include <boost/optional/optional_io.hpp>
       9             : #include <cstddef>
      10             : 
      11             : #include "DataStructures/DenseMatrix.hpp"
      12             : #include "DataStructures/DenseVector.hpp"
      13             : #include "Informer/Verbosity.hpp"
      14             : #include "NumericalAlgorithms/Convergence/Criteria.hpp"
      15             : #include "NumericalAlgorithms/Convergence/HasConverged.hpp"
      16             : #include "NumericalAlgorithms/Convergence/Reason.hpp"
      17             : #include "NumericalAlgorithms/LinearSolver/InnerProduct.hpp"
      18             : #include "Options/Options.hpp"
      19             : #include "Utilities/Gsl.hpp"
      20             : 
      21             : namespace LinearSolver {
      22             : namespace gmres::detail {
      23             : 
      24             : // Perform an Arnoldi orthogonalization to find a new `operand` that is
      25             : // orthogonal to all vectors in `basis_history`. Appends a new column to the
      26             : // `orthogonalization_history` that holds the inner product of the intermediate
      27             : // `operand` with each vector in the `basis_history` and itself.
      28             : template <typename VarsType>
      29             : void arnoldi_orthogonalize(
      30             :     const gsl::not_null<VarsType*> operand,
      31             :     const gsl::not_null<DenseMatrix<double>*> orthogonalization_history,
      32             :     const std::vector<VarsType>& basis_history,
      33             :     const size_t iteration) noexcept {
      34             :   // Resize matrix and make sure the new entries that are not being filled below
      35             :   // are zero.
      36             :   orthogonalization_history->resize(iteration + 2, iteration + 1);
      37             :   for (size_t j = 0; j < iteration; ++j) {
      38             :     (*orthogonalization_history)(iteration + 1, j) = 0.;
      39             :   }
      40             :   // Arnoldi orthogonalization
      41             :   for (size_t j = 0; j < iteration + 1; ++j) {
      42             :     const double orthogonalization = inner_product(basis_history[j], *operand);
      43             :     (*orthogonalization_history)(j, iteration) = orthogonalization;
      44             :     *operand -= orthogonalization * basis_history[j];
      45             :   }
      46             :   (*orthogonalization_history)(iteration + 1, iteration) =
      47             :       sqrt(inner_product(*operand, *operand));
      48             :   // Avoid an FPE if the new operand norm is exactly zero. In that case the
      49             :   // problem is solved and the algorithm will terminate (see Proposition 9.3 in
      50             :   // \cite Saad2003). Since there will be no next iteration we don't need to
      51             :   // normalize the operand.
      52             :   if (UNLIKELY((*orthogonalization_history)(iteration + 1, iteration) == 0.)) {
      53             :     return;
      54             :   }
      55             :   *operand /= (*orthogonalization_history)(iteration + 1, iteration);
      56             : }
      57             : 
      58             : // Solve the linear least-squares problem `||beta - H * y||` for `y`, where `H`
      59             : // is the Hessenberg matrix given by `orthogonalization_history` and `beta` is
      60             : // the vector `(initial_residual, 0, 0, ...)` by updating the QR decomposition
      61             : // of `H` from the previous iteration with a Givens rotation.
      62             : void solve_minimal_residual(
      63             :     gsl::not_null<DenseMatrix<double>*> orthogonalization_history,
      64             :     gsl::not_null<DenseVector<double>*> residual_history,
      65             :     gsl::not_null<DenseVector<double>*> givens_sine_history,
      66             :     gsl::not_null<DenseVector<double>*> givens_cosine_history,
      67             :     size_t iteration) noexcept;
      68             : 
      69             : // Find the vector that minimizes the residual by inverting the upper
      70             : // triangular matrix obtained above.
      71             : DenseVector<double> minimal_residual_vector(
      72             :     const DenseMatrix<double>& orthogonalization_history,
      73             :     const DenseVector<double>& residual_history) noexcept;
      74             : 
      75             : }  // namespace gmres::detail
      76             : 
      77           0 : namespace Serial {
      78             : 
      79             : template <typename VarsType>
      80           0 : struct IdentityPreconditioner {
      81           0 :   VarsType operator()(const VarsType& arg) const noexcept { return arg; }
      82             : };
      83             : 
      84           0 : struct NoIterationCallback {
      85           0 :   void operator()(const Convergence::HasConverged& /*has_converged*/) const
      86             :       noexcept {}
      87             : };
      88             : 
      89             : /*!
      90             :  * \brief A serial GMRES iterative solver for nonsymmetric linear systems of
      91             :  * equations.
      92             :  *
      93             :  * This is an iterative algorithm to solve general linear equations \f$Ax=b\f$
      94             :  * where \f$A\f$ is a linear operator. See \cite Saad2003, chapter 6.5 for a
      95             :  * description of the GMRES algorithm and Algorithm 9.6 for this implementation.
      96             :  * It is matrix-free, which means the operator \f$A\f$ needs not be provided
      97             :  * explicity as a matrix but only the operator action \f$A(x)\f$ must be
      98             :  * provided for an argument \f$x\f$.
      99             :  *
     100             :  * The GMRES algorithm does not require the operator \f$A\f$ to be symmetric or
     101             :  * positive-definite. Note that other algorithms such as conjugate gradients may
     102             :  * be more efficient for symmetric positive-definite operators.
     103             :  *
     104             :  * \par Convergence:
     105             :  * Given a set of \f$N_A\f$ equations (e.g. through an \f$N_A\times N_A\f$
     106             :  * matrix) the GMRES algorithm will converge to numerical precision in at most
     107             :  * \f$N_A\f$ iterations. However, depending on the properties of the linear
     108             :  * operator, an approximate solution can ideally be obtained in only a few
     109             :  * iterations. See \cite Saad2003, section 6.11.4 for details on the convergence
     110             :  * of the GMRES algorithm.
     111             :  *
     112             :  * \par Restarting:
     113             :  * This implementation of the GMRES algorithm supports restarting, as detailed
     114             :  * in \cite Saad2003, section 6.5.5. Since the GMRES algorithm iteratively
     115             :  * builds up an orthogonal basis of the solution space the cost of each
     116             :  * iteration increases linearly with the number of iterations. Therefore it is
     117             :  * sometimes helpful to restart the algorithm every \f$N_\mathrm{restart}\f$
     118             :  * iterations, discarding the set of basis vectors and starting again from the
     119             :  * current solution estimate. This strategy can improve the performance of the
     120             :  * solver, but note that the solver can stagnate for non-positive-definite
     121             :  * operators and is not guaranteed to converge within \f$N_A\f$ iterations
     122             :  * anymore. Set the `restart` argument of the constructor to
     123             :  * \f$N_\mathrm{restart}\f$ to activate restarting, or set it to zero to
     124             :  * deactivate restarting (default behaviour).
     125             :  *
     126             :  * \par Preconditioning:
     127             :  * This implementation of the GMRES algorithm also supports preconditioning.
     128             :  * You can provide a linear operator \f$P\f$ that approximates the inverse of
     129             :  * the operator \f$A\f$ to accelerate the convergence of the linear solve.
     130             :  * The algorithm is right-preconditioned, which allows the preconditioner to
     131             :  * change in every iteration ("flexible" variant). See \cite Saad2003, sections
     132             :  * 9.3.2 and 9.4.1 for details. This implementation follows Algorithm 9.6 in
     133             :  * \cite Saad2003.
     134             :  *
     135             :  * \par Improvements:
     136             :  * Further improvements can potentially be implemented for this algorithm, see
     137             :  * e.g. \cite Ayachour2003.
     138             :  *
     139             :  * \example
     140             :  * \snippet NumericalAlgorithms/LinearSolver/Test_Gmres.cpp gmres_example
     141             :  */
     142             : template <typename VarsType>
     143           1 : struct Gmres {
     144             :  private:
     145           0 :   struct ConvergenceCriteria {
     146           0 :     using type = Convergence::Criteria;
     147           0 :     static constexpr Options::String help =
     148             :         "Determine convergence of the algorithm";
     149             :   };
     150           0 :   struct Restart {
     151           0 :     using type = size_t;
     152           0 :     static constexpr Options::String help =
     153             :         "Iterations to run before restarting";
     154           0 :     static size_t default_value() noexcept { return 0; }
     155             :   };
     156           0 :   struct Verbosity {
     157           0 :     using type = ::Verbosity;
     158           0 :     static constexpr Options::String help = "Logging verbosity";
     159             :   };
     160             : 
     161             :  public:
     162           0 :   static constexpr Options::String help =
     163             :       "A serial GMRES iterative solver for nonsymmetric linear systems of\n"
     164             :       "equations Ax=b. It will converge to numerical precision in at most N_A\n"
     165             :       "iterations, where N_A is the number of equations represented by the\n"
     166             :       "linear operator A, but will ideally converge to a reasonable\n"
     167             :       "approximation of the solution x in only a few iterations.\n"
     168             :       "\n"
     169             :       "Restarting: It is sometimes helpful to restart the algorithm every\n"
     170             :       "N_restart iterations to speed it up. Note that it can stagnate for\n"
     171             :       "non-positive-definite matrices and is not guaranteed to converge\n"
     172             :       "within N_A iterations anymore when restarting is activated.\n"
     173             :       "Activate restarting by setting the 'Restart' option to N_restart, or\n"
     174             :       "deactivate restarting by setting it to zero (default).";
     175           0 :   using options = tmpl::list<ConvergenceCriteria, Verbosity, Restart>;
     176             : 
     177           0 :   Gmres(Convergence::Criteria convergence_criteria, ::Verbosity verbosity,
     178             :         size_t restart = 0) noexcept
     179             :       // clang-tidy: trivially copyable
     180             :       : convergence_criteria_(std::move(convergence_criteria)),  // NOLINT
     181             :         verbosity_(std::move(verbosity)),                        // NOLINT
     182             :         restart_(restart > 0 ? restart : convergence_criteria_.max_iterations) {
     183             :     initialize();
     184             :   }
     185             : 
     186           0 :   Gmres() = default;
     187           0 :   Gmres(const Gmres& /*rhs*/) = default;
     188           0 :   Gmres& operator=(const Gmres& /*rhs*/) = default;
     189           0 :   Gmres(Gmres&& /*rhs*/) = default;
     190           0 :   Gmres& operator=(Gmres&& /*rhs*/) = default;
     191           0 :   ~Gmres() = default;
     192             : 
     193           0 :   void initialize() noexcept {
     194             :     orthogonalization_history_.reserve(restart_ + 1);
     195             :     residual_history_.reserve(restart_ + 1);
     196             :     givens_sine_history_.reserve(restart_);
     197             :     givens_cosine_history_.reserve(restart_);
     198             :     basis_history_.resize(restart_ + 1);
     199             :     preconditioned_basis_history_.resize(restart_);
     200             :   }
     201             : 
     202           0 :   const Convergence::Criteria& convergence_criteria() const noexcept {
     203             :     return convergence_criteria_;
     204             :   }
     205             : 
     206           0 :   void pup(PUP::er& p) noexcept {  // NOLINT
     207             :     p | convergence_criteria_;
     208             :     p | verbosity_;
     209             :     p | restart_;
     210             :     if (p.isUnpacking()) {
     211             :       initialize();
     212             :     }
     213             :   }
     214             : 
     215             :   /*!
     216             :    * \brief Iteratively solve the problem \f$Ax=b\f$ for \f$x\f$ where \f$A\f$
     217             :    * is the `linear_operator` and \f$b\f$ is the `source`, starting \f$x\f$ at
     218             :    * `initial_guess`.
     219             :    *
     220             :    * Optionally provide a `preconditioner` (see class documentation).
     221             :    *
     222             :    * \return An instance of `Convergence::HasConverged` that provides
     223             :    * information on the convergence status of the completed solve, and the
     224             :    * approximate solution \f$x\f$.
     225             :    */
     226             :   template <typename LinearOperator, typename SourceType,
     227             :             typename Preconditioner = IdentityPreconditioner<VarsType>,
     228             :             typename IterationCallback = NoIterationCallback>
     229           1 :   std::pair<Convergence::HasConverged, VarsType> operator()(
     230             :       LinearOperator&& linear_operator, const SourceType& source,
     231             :       const VarsType& initial_guess,
     232             :       Preconditioner&& preconditioner = IdentityPreconditioner<VarsType>{},
     233             :       IterationCallback&& iteration_callback = NoIterationCallback{}) const
     234             :       noexcept;
     235             : 
     236             :  private:
     237           0 :   Convergence::Criteria convergence_criteria_{};
     238           0 :   ::Verbosity verbosity_{::Verbosity::Verbose};
     239           0 :   size_t restart_{};
     240             : 
     241             :   // Memory buffers to avoid re-allocating memory for successive solves:
     242             :   // The `orthogonalization_history_` is built iteratively from inner products
     243             :   // between existing and potential basis vectors and then Givens-rotated to
     244             :   // become upper-triangular.
     245           0 :   mutable DenseMatrix<double> orthogonalization_history_{};
     246             :   // The `residual_history_` holds the remaining residual in its last entry, and
     247             :   // the other entries `g` "source" the minimum residual vector `y` in
     248             :   // `R * y = g` where `R` is the upper-triangular `orthogonalization_history_`.
     249           0 :   mutable DenseVector<double> residual_history_{};
     250             :   // These represent the accumulated Givens rotations up to the current
     251             :   // iteration.
     252           0 :   mutable DenseVector<double> givens_sine_history_{};
     253           0 :   mutable DenseVector<double> givens_cosine_history_{};
     254             :   // These represent the orthogonal Krylov-subspace basis that is constructed
     255             :   // iteratively by Arnoldi-orthogonalizing a new vector in each iteration and
     256             :   // appending it to the `basis_history_`.
     257           0 :   mutable std::vector<VarsType> basis_history_{};
     258             :   // When a preconditioner is used it is applied to each new basis vector. The
     259             :   // preconditioned basis is used to construct the solution when the algorithm
     260             :   // has converged.
     261           0 :   mutable std::vector<VarsType> preconditioned_basis_history_{};
     262             : };
     263             : 
     264             : template <typename VarsType>
     265             : template <typename LinearOperator, typename SourceType, typename Preconditioner,
     266             :           typename IterationCallback>
     267             : std::pair<Convergence::HasConverged, VarsType> Gmres<VarsType>::operator()(
     268             :     LinearOperator&& linear_operator, const SourceType& source,
     269             :     const VarsType& initial_guess, Preconditioner&& preconditioner,
     270             :     IterationCallback&& iteration_callback) const noexcept {
     271             :   constexpr bool use_preconditioner =
     272             :       not std::is_same_v<Preconditioner, IdentityPreconditioner<VarsType>>;
     273             :   constexpr bool use_iteration_callback =
     274             :       not std::is_same_v<IterationCallback, NoIterationCallback>;
     275             : 
     276             :   auto result = initial_guess;
     277             :   Convergence::HasConverged has_converged{};
     278             :   size_t iteration = 0;
     279             : 
     280             :   while (not has_converged) {
     281             :     auto& initial_operand = basis_history_[0] =
     282             :         source - linear_operator(result);
     283             :     const double initial_residual_magnitude =
     284             :         sqrt(inner_product(initial_operand, initial_operand));
     285             :     has_converged = Convergence::HasConverged{convergence_criteria_, iteration,
     286             :                                               initial_residual_magnitude,
     287             :                                               initial_residual_magnitude};
     288             :     if (use_iteration_callback) {
     289             :       iteration_callback(has_converged);
     290             :     }
     291             :     if (UNLIKELY(has_converged)) {
     292             :       break;
     293             :     }
     294             :     initial_operand /= initial_residual_magnitude;
     295             :     residual_history_.resize(1);
     296             :     residual_history_[0] = initial_residual_magnitude;
     297             :     for (size_t k = 0; k < restart_; ++k) {
     298             :       auto& operand = basis_history_[k + 1];
     299             :       if (use_preconditioner) {
     300             :         preconditioned_basis_history_[k] = preconditioner(basis_history_[k]);
     301             :       }
     302             :       operand =
     303             :           linear_operator(use_preconditioner ? preconditioned_basis_history_[k]
     304             :                                              : basis_history_[k]);
     305             :       // Find a new orthogonal basis vector of the Krylov subspace
     306             :       gmres::detail::arnoldi_orthogonalize(
     307             :           make_not_null(&operand), make_not_null(&orthogonalization_history_),
     308             :           basis_history_, k);
     309             :       // Least-squares solve for the minimal residual
     310             :       gmres::detail::solve_minimal_residual(
     311             :           make_not_null(&orthogonalization_history_),
     312             :           make_not_null(&residual_history_),
     313             :           make_not_null(&givens_sine_history_),
     314             :           make_not_null(&givens_cosine_history_), k);
     315             :       ++iteration;
     316             :       has_converged = Convergence::HasConverged{
     317             :           convergence_criteria_, iteration, abs(residual_history_[k + 1]),
     318             :           initial_residual_magnitude};
     319             :       if (use_iteration_callback) {
     320             :         iteration_callback(has_converged);
     321             :       }
     322             :       if (UNLIKELY(has_converged)) {
     323             :         break;
     324             :       }
     325             :     }
     326             :     // Find the vector w.r.t. the constructed orthogonal basis of the Krylov
     327             :     // subspace that minimizes the residual
     328             :     const auto minres = gmres::detail::minimal_residual_vector(
     329             :         orthogonalization_history_, residual_history_);
     330             :     // Construct the solution from the orthogonal basis and the minimal residual
     331             :     // vector
     332             :     for (size_t i = 0; i < minres.size(); ++i) {
     333             :       result +=
     334             :           minres[i] * gsl::at(use_preconditioner ? preconditioned_basis_history_
     335             :                                                  : basis_history_,
     336             :                               i);
     337             :     }
     338             :   }
     339             :   // NOLINTNEXTLINE(hicpp-move-const-arg,performance-move-const-arg)
     340             :   return {std::move(has_converged), std::move(result)};
     341             : }
     342             : 
     343             : }  // namespace Serial
     344             : }  // namespace LinearSolver

Generated by: LCOV version 1.14