From a4917cfe85a4f797a97443fca1e3d16decfd0d57 Mon Sep 17 00:00:00 2001 From: Matthias Maier Date: Wed, 26 Mar 2025 16:32:40 -0500 Subject: [PATCH] SolverMPGMRES: implement a static dispatch strategy --- include/deal.II/lac/solver_gmres.h | 92 ++++++++++++++++++++---------- 1 file changed, 63 insertions(+), 29 deletions(-) diff --git a/include/deal.II/lac/solver_gmres.h b/include/deal.II/lac/solver_gmres.h index 319cc160ab..5efac2498d 100644 --- a/include/deal.II/lac/solver_gmres.h +++ b/include/deal.II/lac/solver_gmres.h @@ -37,6 +37,7 @@ #include #include #include +#include #include DEAL_II_NAMESPACE_OPEN @@ -695,25 +696,23 @@ public: /** * Solve the linear system $Ax=b$ for x. */ - 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); + template + void + solve(const MatrixType &A, + VectorType &x, + const VectorType &b, + const PreconditionerTypes &...preconditioners); protected: /** * Solve the linear system $Ax=b$ for x. */ - template + template void - solve_internal(const MatrixType &A, - VectorType &x, - const VectorType &b, - const PreconditionerType &preconditioner); + solve_internal(const MatrixType &A, + VectorType &x, + const VectorType &b, + const PreconditionerTypes &...preconditioners); private: /** @@ -2053,31 +2052,66 @@ SolverMPGMRES::SolverMPGMRES(SolverControl &cn, 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 SolverMPGMRES::solve(const MatrixType &A, - VectorType &x, - const VectorType &b, - const PreconditionerType &preconditioner) +template +void +SolverMPGMRES::solve(const MatrixType &A, + VectorType &x, + const VectorType &b, + const PreconditionerTypes &...preconditioners) { LogStream::Prefix prefix("MPGMRES"); - SolverMPGMRES::solve_internal(A, x, b, preconditioner); + SolverMPGMRES::solve_internal(A, x, b, preconditioners...); } template -template +template void SolverMPGMRES::solve_internal( - const MatrixType &A, - VectorType &x, - const VectorType &b, - const PreconditionerType &preconditioner) + const MatrixType &A, + VectorType &x, + const VectorType &b, + const PreconditionerTypes &...preconditioners) { + constexpr std::size_t n_preconditioners = sizeof...(PreconditionerTypes); + + // + // With each invocation of preconditioner_vmult we increase the + // variable current_index by one modulo n_preconditioners. + // + std::size_t current_index = 0; + + // + // A lambda that cycles through all preconditioners in sequence applying + // one to the vector src and storing the result in dst: + // + const auto preconditioner_vmult = [&](auto &dst, const auto &src) { + // We have no preconditioner that we could apply + if (n_preconditioners == 0) + { + dst = src; + return; + } + + // + // We cycle through all preconditioners and call the one with a + // matching index: + // + std::size_t n = 0; + + const auto call_current_preconditioner = [&](const auto &preconditioner) { + if (n++ == current_index) + preconditioner.vmult(dst, src); + }; + + // https://en.cppreference.com/w/cpp/language/fold + (call_current_preconditioner(preconditioners), ...); + + current_index = (current_index + 1) % n_preconditioners; + }; + + SolverControl::State iteration_state = SolverControl::iterate; const unsigned int basis_size = additional_data.max_basis_size; @@ -2126,7 +2160,7 @@ SolverMPGMRES::solve_internal( iteration_state == SolverControl::iterate); ++inner_iteration) { - preconditioner.vmult(z(inner_iteration, x), v[inner_iteration]); + preconditioner_vmult(z(inner_iteration, x), v[inner_iteration]); A.vmult(v(inner_iteration + 1, x), z[inner_iteration]); res = -- 2.39.5