From 6bccdba496c3ea324152d3264f9d76a2ad2b3c12 Mon Sep 17 00:00:00 2001 From: Wyatt Smith <216wsmith@gmail.com> Date: Wed, 26 Mar 2025 12:56:44 -0500 Subject: [PATCH] add SolverMPGMRES skeleton, refactor code --- doc/doxygen/references.bib | 13 ++ include/deal.II/lac/solver_gmres.h | 196 +++++++++++++++++++++++------ 2 files changed, 172 insertions(+), 37 deletions(-) diff --git a/doc/doxygen/references.bib b/doc/doxygen/references.bib index a55bc59367..ce9007f0d9 100644 --- a/doc/doxygen/references.bib +++ b/doc/doxygen/references.bib @@ -2066,6 +2066,19 @@ Url = {citeseer.ist.psu.edu/saad93flexible.html} } +@article{Greif2017, + title = {GMRES with multiple preconditioners}, + author = {Greif, Chen and Rees, Tyrone and Szyld, Daniel B}, + journal = {SeMA Journal}, + volume = {74}, + pages = {213--231}, + year = {2017}, + publisher = {Springer}, + doi = {10.1007/s40324-016-0088-7}, + url = {https://doi.org/10.1007/s40324-016-0088-7} +} + + @article{ainsworth1998hp, author = {M. Ainsworth and B. Senior}, title = {An adaptive refinement strategy for \textit{hp}-finite element computations}, diff --git a/include/deal.II/lac/solver_gmres.h b/include/deal.II/lac/solver_gmres.h index 396014161e..319cc160ab 100644 --- a/include/deal.II/lac/solver_gmres.h +++ b/include/deal.II/lac/solver_gmres.h @@ -127,11 +127,12 @@ namespace internal /** * Class that performs the Arnoldi orthogonalization process within the - * SolverGMRES and SolverFGMRES classes. It uses one of the algorithms in - * LinearAlgebra::OrthogonalizationStrategy for the work on the global - * vectors, transforms the resulting Hessenberg matrix into an upper - * triangular matrix by Givens rotations, and eventually solves the - * minimization problem in the projected Krylov space. + * SolverGMRES, SolverFGMRES, and SolverMPGMRES classes. It uses one of + * the algorithms in LinearAlgebra::OrthogonalizationStrategy for the + * work on the global vectors, transforms the resulting Hessenberg + * matrix into an upper triangular matrix by Givens rotations, and + * eventually solves the minimization problem in the projected Krylov + * space. */ template class ArnoldiProcess @@ -632,28 +633,21 @@ protected: /** - * Implementation of the Generalized minimal residual method with flexible - * preconditioning (flexible GMRES or FGMRES). + * Implementation of the multiply preconditioned generalized minimal + * residual method (MPGMRES). * - * This flexible version of the GMRES method allows for the use of a different - * preconditioner in each iteration step. Therefore, it is also more robust - * with respect to inaccurate evaluation of the preconditioner. An important - * application is the use of a Krylov space method inside the - * preconditioner. As opposed to SolverGMRES which allows one to choose - * between left and right preconditioning, this solver always applies the - * preconditioner from the right. + * This is a variant of the flexible GMRES method that uses multiple + * preconditioners to construct independent Krylov subspaces, one for each + * preconditioner. In contrast, the flexible GMRES method implemented in + * SolverFGMRES constructs only one "Krylov subspace." * - * FGMRES needs two vectors in each iteration steps yielding a total of - * 2*SolverFGMRES::%AdditionalData::%max_basis_size+1 auxiliary - * vectors. Otherwise, FGMRES requires roughly the same number of operations - * per iteration compared to GMRES, except one application of the - * preconditioner less at each restart and at the end of solve(). + * TODO add longer description * - * For more details see @cite Saad1991. + * For more details see @cite Greif2017. */ template > DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector) -class SolverFGMRES : public SolverBase +class SolverMPGMRES : public SolverBase { public: /** @@ -687,16 +681,16 @@ public: /** * Constructor. */ - SolverFGMRES(SolverControl &cn, - VectorMemory &mem, - const AdditionalData &data = AdditionalData()); + SolverMPGMRES(SolverControl &cn, + VectorMemory &mem, + const AdditionalData &data = AdditionalData()); /** * Constructor. Use an object of type GrowingVectorMemory as a default to * allocate memory. */ - SolverFGMRES(SolverControl &cn, - const AdditionalData &data = AdditionalData()); + SolverMPGMRES(SolverControl &cn, + const AdditionalData &data = AdditionalData()); /** * Solve the linear system $Ax=b$ for x. @@ -710,6 +704,17 @@ public: const VectorType &b, const PreconditionerType &preconditioner); +protected: + /** + * Solve the linear system $Ax=b$ for x. + */ + template + void + solve_internal(const MatrixType &A, + VectorType &x, + const VectorType &b, + const PreconditionerType &preconditioner); + private: /** * Additional flags. @@ -725,10 +730,71 @@ private: arnoldi_process; }; + + +/** + * Implementation of the generalized minimal residual method with flexible + * preconditioning (flexible GMRES or FGMRES). + * + * This flexible version of the GMRES method allows for the use of a different + * preconditioner in each iteration step. Therefore, it is also more robust + * with respect to inaccurate evaluation of the preconditioner. An important + * application is the use of a Krylov space method inside the + * preconditioner. As opposed to SolverGMRES which allows one to choose + * between left and right preconditioning, this solver always applies the + * preconditioner from the right. + * + * FGMRES needs two vectors in each iteration steps yielding a total of + * 2*SolverFGMRES::%AdditionalData::%max_basis_size+1 auxiliary + * vectors. Otherwise, FGMRES requires roughly the same number of operations + * per iteration compared to GMRES, except one application of the + * preconditioner less at each restart and at the end of solve(). + * + * For more details see @cite Saad1991. + */ +template > +DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector) +class SolverFGMRES : public SolverMPGMRES +{ +public: + using AdditionalData = typename SolverMPGMRES::AdditionalData; + + /** + * Constructor. + */ + SolverFGMRES(SolverControl &cn, + VectorMemory &mem, + const AdditionalData &data = AdditionalData()); + + /** + * Constructor. Use an object of type GrowingVectorMemory as a default to + * allocate memory. + */ + SolverFGMRES(SolverControl &cn, + const AdditionalData &data = AdditionalData()); + + /** + * Solve the linear system $Ax=b$ for x. + * + * @note If you want to use more than one preconditioner, then you will + * need to supply a @p preconditioner object that switches + * preconditioners for each vmult() invocation. + */ + template + DEAL_II_CXX20_REQUIRES( + (concepts::is_linear_operator_on && + concepts::is_linear_operator_on)) + void solve(const MatrixType &A, + VectorType &x, + const VectorType &b, + const PreconditionerType &preconditioner); +}; + /** @} */ -/* --------------------- Inline and template functions ------------------- */ +/* --------------------- Inline and template functions ------------------- */ + #ifndef DOXYGEN template @@ -1963,11 +2029,13 @@ double SolverGMRES::criterion() //----------------------------------------------------------------------// + + template DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector) -SolverFGMRES::SolverFGMRES(SolverControl &cn, - VectorMemory &mem, - const AdditionalData &data) +SolverMPGMRES::SolverMPGMRES(SolverControl &cn, + VectorMemory &mem, + const AdditionalData &data) : SolverBase(cn, mem) , additional_data(data) {} @@ -1976,8 +2044,8 @@ SolverFGMRES::SolverFGMRES(SolverControl &cn, template DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector) -SolverFGMRES::SolverFGMRES(SolverControl &cn, - const AdditionalData &data) +SolverMPGMRES::SolverMPGMRES(SolverControl &cn, + const AdditionalData &data) : SolverBase(cn) , additional_data(data) {} @@ -1990,13 +2058,26 @@ template DEAL_II_CXX20_REQUIRES( (concepts::is_linear_operator_on && concepts::is_linear_operator_on)) -void SolverFGMRES::solve(const MatrixType &A, - VectorType &x, - const VectorType &b, - const PreconditionerType &preconditioner) +void SolverMPGMRES::solve(const MatrixType &A, + VectorType &x, + const VectorType &b, + const PreconditionerType &preconditioner) { - LogStream::Prefix prefix("FGMRES"); + LogStream::Prefix prefix("MPGMRES"); + SolverMPGMRES::solve_internal(A, x, b, preconditioner); +} + + +template +template +void +SolverMPGMRES::solve_internal( + const MatrixType &A, + VectorType &x, + const VectorType &b, + const PreconditionerType &preconditioner) +{ SolverControl::State iteration_state = SolverControl::iterate; const unsigned int basis_size = additional_data.max_basis_size; @@ -2079,6 +2160,47 @@ void SolverFGMRES::solve(const MatrixType &A, SolverControl::NoConvergence(accumulated_iterations, res)); } + + +//----------------------------------------------------------------------// + + + +template +DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector) +SolverFGMRES::SolverFGMRES(SolverControl &cn, + VectorMemory &mem, + const AdditionalData &data) + : SolverMPGMRES(cn, mem, data) +{} + + + +template +DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector) +SolverFGMRES::SolverFGMRES(SolverControl &cn, + const AdditionalData &data) + : SolverMPGMRES(cn, data) +{} + + + +template +DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector) +template +DEAL_II_CXX20_REQUIRES( + (concepts::is_linear_operator_on && + concepts::is_linear_operator_on)) +void SolverFGMRES::solve(const MatrixType &A, + VectorType &x, + const VectorType &b, + const PreconditionerType &preconditioner) +{ + LogStream::Prefix prefix("FGMRES"); + SolverMPGMRES::solve_internal(A, x, b, preconditioner); +} + + #endif // DOXYGEN DEAL_II_NAMESPACE_CLOSE -- 2.39.5