From: Martin Kronbichler <martin.kronbichler@rub.de>
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<VectorType>::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<VectorType>::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<VectorType>::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
     }