From c3cf784368e9088ed23f41d0a5cfe4e3fcec5edc Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Fri, 14 Oct 2022 16:01:31 +0200 Subject: [PATCH] Refactor SolverGMRES::modified_gram_schmidt --- include/deal.II/lac/solver_gmres.h | 91 ++++++++++++++++-------------- 1 file changed, 48 insertions(+), 43 deletions(-) diff --git a/include/deal.II/lac/solver_gmres.h b/include/deal.II/lac/solver_gmres.h index 9b21890a24..74891ad5ca 100644 --- a/include/deal.II/lac/solver_gmres.h +++ b/include/deal.II/lac/solver_gmres.h @@ -418,7 +418,7 @@ protected: * Calls the signal re_orthogonalize_signal if it is connected. */ static double - modified_gram_schmidt( + iterated_modified_gram_schmidt( const internal::SolverGMRESImplementation::TmpVectors & orthogonal_vectors, const unsigned int dim, @@ -712,7 +712,7 @@ SolverGMRES::givens_rotation(Vector &h, template inline double -SolverGMRES::modified_gram_schmidt( +SolverGMRES::iterated_modified_gram_schmidt( const internal::SolverGMRESImplementation::TmpVectors & orthogonal_vectors, const unsigned int dim, @@ -732,40 +732,12 @@ SolverGMRES::modified_gram_schmidt( if (consider_reorthogonalize) norm_vv_start = vv.l2_norm(); - // Orthogonalization - h(0) = vv * orthogonal_vectors[0]; - for (unsigned int i = 1; i < dim; ++i) - h(i) = vv.add_and_dot(-h(i - 1), - orthogonal_vectors[i - 1], - orthogonal_vectors[i]); - double norm_vv = - std::sqrt(vv.add_and_dot(-h(dim - 1), orthogonal_vectors[dim - 1], vv)); - - // Re-orthogonalization if loss of orthogonality detected. For the test, use - // a strategy discussed in C. T. Kelley, Iterative Methods for Linear and - // Nonlinear Equations, SIAM, Philadelphia, 1995: Compare the norm of vv - // after orthogonalization with its norm when starting the - // orthogonalization. If vv became very small (here: less than the square - // root of the machine precision times 10), it is almost in the span of the - // previous vectors, which indicates loss of precision. - if (consider_reorthogonalize) - { - if (norm_vv > - 10. * norm_vv_start * - std::sqrt( - std::numeric_limits::epsilon())) - return norm_vv; + for (unsigned int i = 0; i < dim; ++i) + h[i] = 0; - else - { - reorthogonalize = true; - if (!reorthogonalize_signal.empty()) - reorthogonalize_signal(accumulated_iterations); - } - } - - if (reorthogonalize == true) + for (unsigned int c = 0; c < 2; ++c) // 0: orthogonalize, 1: reorthogonalize { + // Orthogonalization double htmp = vv * orthogonal_vectors[0]; h(0) += htmp; for (unsigned int i = 1; i < dim; ++i) @@ -775,11 +747,43 @@ SolverGMRES::modified_gram_schmidt( orthogonal_vectors[i]); h(i) += htmp; } - norm_vv = + + double norm_vv = std::sqrt(vv.add_and_dot(-htmp, orthogonal_vectors[dim - 1], vv)); + + if (c == 1) + return norm_vv; // reorthogonalization already performed -> finished + + // Re-orthogonalization if loss of orthogonality detected. For the test, + // use a strategy discussed in C. T. Kelley, Iterative Methods for Linear + // and Nonlinear Equations, SIAM, Philadelphia, 1995: Compare the norm of + // vv after orthogonalization with its norm when starting the + // orthogonalization. If vv became very small (here: less than the square + // root of the machine precision times 10), it is almost in the span of + // the previous vectors, which indicates loss of precision. + if (consider_reorthogonalize) + { + if (norm_vv > + 10. * norm_vv_start * + std::sqrt(std::numeric_limits< + typename VectorType::value_type>::epsilon())) + return norm_vv; + + else + { + reorthogonalize = true; + if (!reorthogonalize_signal.empty()) + reorthogonalize_signal(accumulated_iterations); + } + } + + if (reorthogonalize == false) + return norm_vv; // no reorthogonalization needed -> finished } - return norm_vv; + AssertThrow(false, ExcInternalError()); + + return 0.0; } @@ -1013,13 +1017,14 @@ SolverGMRES::solve(const MatrixType & A, dim = inner_iteration + 1; - const double s = modified_gram_schmidt(tmp_vectors, - dim, - accumulated_iterations, - vv, - h, - re_orthogonalize, - re_orthogonalize_signal); + const double s = + iterated_modified_gram_schmidt(tmp_vectors, + dim, + accumulated_iterations, + vv, + h, + re_orthogonalize, + re_orthogonalize_signal); h(inner_iteration + 1) = s; // s=0 is a lucky breakdown, the solver will reach convergence, -- 2.39.5