From 1ef9febfd42d7a9d67eacaea72a6f5e332c92a3e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= <oberhuber.tomas@gmail.com>
Date: Wed, 22 Dec 2021 16:25:06 +0100
Subject: [PATCH] Writtin documentation for GMRES method.
---
src/TNL/Solvers/Linear/GMRES.h | 341 ++++++++++++++++++++-------------
1 file changed, 209 insertions(+), 132 deletions(-)
diff --git a/src/TNL/Solvers/Linear/GMRES.h b/src/TNL/Solvers/Linear/GMRES.h
index 7314a3a05b..7ff51cbb98 100644
--- a/src/TNL/Solvers/Linear/GMRES.h
+++ b/src/TNL/Solvers/Linear/GMRES.h
@@ -12,12 +12,33 @@
#pragma once
-#include "LinearSolver.h"
+#include <TNL/Solvers/Linear/LinearSolver.h>
namespace TNL {
-namespace Solvers {
-namespace Linear {
-
+ namespace Solvers {
+ namespace Linear {
+
+/**
+ * \brief Iterative solver of linear systems based on the Generalized minimal residual (GMRES) method.
+ *
+ * This method can be used for solving of non-semmetric linear systems. This implementation
+ * offers various methods for the orthogonalization of vectors generated by the Arnoldi algorithm.
+ * The user may choose one of the following:
+ *
+ * \e CGS - Classical Gramm-Schmidt,
+ *
+ * \e CGSR - Classical Gramm-Schmidt with reorthogonalization,
+ *
+ * \e MGS - Modified Gramm-Schmidt,
+ *
+ * \e MGSR - Modified Gramm-Schmidt with reorthogonalization (the default one),
+ *
+ * \e CWY - Compact WY form of the Householder reflections).
+ *
+ * See [Wikipedia](https://en.wikipedia.org/wiki/Generalized_minimal_residual_method) for more details.
+ *
+ * \tparam Matrix is type of matrix describing the linear system.
+ */
template< typename Matrix >
class GMRES
: public LinearSolver< Matrix >
@@ -25,139 +46,195 @@ class GMRES
using Base = LinearSolver< Matrix >;
using Traits = Linear::Traits< Matrix >;
-public:
- using RealType = typename Base::RealType;
- using DeviceType = typename Base::DeviceType;
- using IndexType = typename Base::IndexType;
- // distributed vectors/views
- using VectorViewType = typename Base::VectorViewType;
- using ConstVectorViewType = typename Base::ConstVectorViewType;
- using VectorType = typename Traits::VectorType;
-
- static void configSetup( Config::ConfigDescription& config,
- const String& prefix = "" );
-
- bool setup( const Config::ParameterContainer& parameters,
- const String& prefix = "" ) override;
-
- bool solve( ConstVectorViewType b, VectorViewType x ) override;
-
-protected:
- // local vectors/views
- using ConstDeviceView = typename Traits::ConstLocalViewType;
- using DeviceView = typename Traits::LocalViewType;
- using DeviceVector = typename Traits::LocalVectorType;
- using HostView = typename DeviceView::template Self< RealType, Devices::Host >;
- using HostVector = typename DeviceVector::template Self< RealType, Devices::Host >;;
-
- enum class Variant { CGS, CGSR, MGS, MGSR, CWY };
-
-// nvcc allows __cuda_callable__ lambdas only in public methods
-#ifdef __NVCC__
-public:
-#endif
- int orthogonalize_CGS( const int m, const RealType normb, const RealType beta );
-#ifdef __NVCC__
-protected:
-#endif
-
- int orthogonalize_MGS( const int m, const RealType normb, const RealType beta );
-
- int orthogonalize_CWY( const int m, const RealType normb, const RealType beta );
-
- void compute_residue( VectorViewType r, ConstVectorViewType x, ConstVectorViewType b );
-
- void preconditioned_matvec( VectorViewType w, ConstVectorViewType v );
-
-// nvcc allows __cuda_callable__ lambdas only in public methods
-#ifdef __NVCC__
-public:
-#endif
- void hauseholder_generate( const int i,
- VectorViewType y_i,
- ConstVectorViewType z );
-#ifdef __NVCC__
-protected:
-#endif
-
- void hauseholder_apply_trunc( HostView out,
- const int i,
+ public:
+
+ /**
+ * \brief Floating point type used for computations.
+ */
+ using RealType = typename Base::RealType;
+
+ /**
+ * \brief Device where the solver will run on and auxillary data will alloacted on.
+ */
+ using DeviceType = typename Base::DeviceType;
+
+ /**
+ * \brief Type for indexing.
+ */
+ using IndexType = typename Base::IndexType;
+
+ /**
+ * \brief Type for vector view.
+ */
+ using VectorViewType = typename Base::VectorViewType;
+
+ /**
+ * \brief Type for constant vector view.
+ */
+ using ConstVectorViewType = typename Base::ConstVectorViewType;
+
+ /**
+ * \brief This is method defines configuration entries for setup of the linear iterative solver.
+ *
+ * In addition to config entries defined by \ref IterativeSolver::configSetup, this method
+ * defines the following:
+ *
+ * \e gmres-variant - Algorithm used for the orthogonalization - CGS (Classical Gramm-Schmidt),
+ * CGSR(Classical Gramm-Schmidt with reorthogonalization), MGS(Modified Gramm-Schmidt),
+ * MGSR(Modified Gramm-Schmidt with reorthogonalization), CWY(Compact WY form of the Householder reflections).
+ *
+ * \e gmres-restarting-min - minimal number of iterations after which the GMRES restarts.
+ *
+ * \e gmres-restarting-max - maximal number of iterations after which the GMRES restarts.
+ *
+ * \e gmres-restarting-step-min - minimal adjusting step for the adaptivity of the GMRES restarting parameter.
+ *
+ * \e gmres-restarting-step-max - maximal adjusting step for the adaptivity of the GMRES restarting parameter.
+ *
+ * \param config contains description of configuration parameters.
+ * \param prefix is a prefix of particular configuration entries.
+ */
+ static void configSetup( Config::ConfigDescription& config,
+ const String& prefix = "" );
+
+ /**
+ * \brief Method for setup of the linear iterative solver based on configuration parameters.
+ *
+ * \param parameters contains values of the define configuration entries.
+ * \param prefix is a prefix of particular configuration entries.
+ */
+ bool setup( const Config::ParameterContainer& parameters,
+ const String& prefix = "" ) override;
+
+ /**
+ * \brief Method for solving of a linear system.
+ *
+ * See \ref LinearSolver::solve for more details.
+ *
+ * \param b vector with the right-hand side of the linear system.
+ * \param x vector for the solution of the linear system.
+ * \return true if the solver converged.
+ * \return false if the solver did not converge.
+ */
+ bool solve( ConstVectorViewType b, VectorViewType x ) override;
+
+ protected:
+ using VectorType = typename Traits::VectorType;
+ // local vectors/views
+ using ConstDeviceView = typename Traits::ConstLocalViewType;
+ using DeviceView = typename Traits::LocalViewType;
+ using DeviceVector = typename Traits::LocalVectorType;
+ using HostView = typename DeviceView::template Self< RealType, Devices::Host >;
+ using HostVector = typename DeviceVector::template Self< RealType, Devices::Host >;;
+
+ enum class Variant { CGS, CGSR, MGS, MGSR, CWY };
+
+ // nvcc allows __cuda_callable__ lambdas only in public methods
+ #ifdef __NVCC__
+ public:
+ #endif
+ int orthogonalize_CGS( const int m, const RealType normb, const RealType beta );
+ #ifdef __NVCC__
+ protected:
+ #endif
+
+ int orthogonalize_MGS( const int m, const RealType normb, const RealType beta );
+
+ int orthogonalize_CWY( const int m, const RealType normb, const RealType beta );
+
+ void compute_residue( VectorViewType r, ConstVectorViewType x, ConstVectorViewType b );
+
+ void preconditioned_matvec( VectorViewType w, ConstVectorViewType v );
+
+ // nvcc allows __cuda_callable__ lambdas only in public methods
+ #ifdef __NVCC__
+ public:
+ #endif
+ void hauseholder_generate( const int i,
VectorViewType y_i,
ConstVectorViewType z );
+ #ifdef __NVCC__
+ protected:
+ #endif
- void hauseholder_cwy( VectorViewType v,
- const int i );
-
-// nvcc allows __cuda_callable__ lambdas only in public methods
-#ifdef __NVCC__
-public:
-#endif
- void hauseholder_cwy_transposed( VectorViewType z,
+ void hauseholder_apply_trunc( HostView out,
const int i,
- ConstVectorViewType w );
-#ifdef __NVCC__
-protected:
-#endif
-
- template< typename Vector >
- void update( const int k,
- const int m,
- const HostVector& H,
- const HostVector& s,
- DeviceVector& V,
- Vector& x );
-
- void generatePlaneRotation( RealType& dx,
- RealType& dy,
- RealType& cs,
- RealType& sn );
-
- void applyPlaneRotation( RealType& dx,
- RealType& dy,
- RealType& cs,
- RealType& sn );
-
- void apply_givens_rotations( const int i, const int m );
-
- void setSize( const VectorViewType& x );
-
- // Specialized methods to distinguish between normal and distributed matrices
- // in the implementation.
- template< typename M >
- static IndexType getLocalOffset( const M& m )
- {
- return 0;
- }
-
- template< typename M >
- static IndexType getLocalOffset( const Matrices::DistributedMatrix< M >& m )
- {
- return m.getLocalRowRange().getBegin();
- }
-
- // selected GMRES variant
- Variant variant = Variant::CWY;
-
- // single vectors (distributed)
- VectorType r, w, z, _M_tmp;
- // matrices (in column-major format) (local)
- DeviceVector V, Y;
- // (CWY only) duplicate of the upper (m+1)x(m+1) submatrix of Y (it is lower triangular) for fast access
- HostVector YL, T;
- // host-only storage for Givens rotations and the least squares problem
- HostVector cs, sn, H, s;
-
- IndexType size = 0;
- IndexType ldSize = 0;
- IndexType localOffset = 0;
- int restarting_min = 10;
- int restarting_max = 10;
- int restarting_step_min = 3;
- int restarting_step_max = 3;
+ VectorViewType y_i,
+ ConstVectorViewType z );
+
+ void hauseholder_cwy( VectorViewType v,
+ const int i );
+
+ // nvcc allows __cuda_callable__ lambdas only in public methods
+ #ifdef __NVCC__
+ public:
+ #endif
+ void hauseholder_cwy_transposed( VectorViewType z,
+ const int i,
+ ConstVectorViewType w );
+ #ifdef __NVCC__
+ protected:
+ #endif
+
+ template< typename Vector >
+ void update( const int k,
+ const int m,
+ const HostVector& H,
+ const HostVector& s,
+ DeviceVector& V,
+ Vector& x );
+
+ void generatePlaneRotation( RealType& dx,
+ RealType& dy,
+ RealType& cs,
+ RealType& sn );
+
+ void applyPlaneRotation( RealType& dx,
+ RealType& dy,
+ RealType& cs,
+ RealType& sn );
+
+ void apply_givens_rotations( const int i, const int m );
+
+ void setSize( const VectorViewType& x );
+
+ // Specialized methods to distinguish between normal and distributed matrices
+ // in the implementation.
+ template< typename M >
+ static IndexType getLocalOffset( const M& m )
+ {
+ return 0;
+ }
+
+ template< typename M >
+ static IndexType getLocalOffset( const Matrices::DistributedMatrix< M >& m )
+ {
+ return m.getLocalRowRange().getBegin();
+ }
+
+ // selected GMRES variant
+ Variant variant = Variant::CWY;
+
+ // single vectors (distributed)
+ VectorType r, w, z, _M_tmp;
+ // matrices (in column-major format) (local)
+ DeviceVector V, Y;
+ // (CWY only) duplicate of the upper (m+1)x(m+1) submatrix of Y (it is lower triangular) for fast access
+ HostVector YL, T;
+ // host-only storage for Givens rotations and the least squares problem
+ HostVector cs, sn, H, s;
+
+ IndexType size = 0;
+ IndexType ldSize = 0;
+ IndexType localOffset = 0;
+ int restarting_min = 10;
+ int restarting_max = 10;
+ int restarting_step_min = 3;
+ int restarting_step_max = 3;
};
-} // namespace Linear
-} // namespace Solvers
+ } // namespace Linear
+ } // namespace Solvers
} // namespace TNL
-#include "GMRES.hpp"
+#include <TNL/Solvers/Linear/GMRES.hpp>
--
GitLab