From: Martin Kronbichler Date: Mon, 11 Mar 2024 12:29:52 +0000 (+0100) Subject: SolverGMRES: Fix loop bound for orthogonalization X-Git-Tag: v9.6.0-rc1~486^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=597da73acd9265589db9858a3747d0c904e99a9d;p=dealii.git SolverGMRES: Fix loop bound for orthogonalization --- diff --git a/include/deal.II/lac/solver_gmres.h b/include/deal.II/lac/solver_gmres.h index 7ac8083e12..1941456675 100644 --- a/include/deal.II/lac/solver_gmres.h +++ b/include/deal.II/lac/solver_gmres.h @@ -930,10 +930,11 @@ namespace internal { Assert(dim > 0, ExcInternalError()); - for (unsigned int i = 0; i < dim; ++i) + for (unsigned int i = 0; i < dim - 1; ++i) vv.add(-h(i), orthogonal_vectors[i]); - return std::sqrt(vv.add_and_dot(-h(dim), orthogonal_vectors[dim], vv)); + return std::sqrt( + vv.add_and_dot(-h(dim - 1), orthogonal_vectors[dim - 1], vv)); } @@ -1366,6 +1367,8 @@ SolverGMRES::solve(const MatrixType &A, } bool re_orthogonalize = additional_data.force_re_orthogonalization; + // reset this vector to the right size + h.reinit(n_tmp_vectors - 1); /////////////////////////////////////////////////////////////////////////// // outer iteration: loop until we either reach convergence or the maximum @@ -1373,9 +1376,6 @@ SolverGMRES::solve(const MatrixType &A, // restart do { - // reset this vector to the right size - h.reinit(n_tmp_vectors - 1); - double rho = 0.0; if (left_precondition) @@ -1577,7 +1577,7 @@ SolverGMRES::solve(const MatrixType &A, p, dim, h, tmp_vectors, true); preconditioner.vmult(v, p); x.add(1., v); - }; + } // end of outer iteration. restart if no convergence and the number of // iterations is not exceeded }