From d9554eb4a8ba5e31a01c9744dd5119b43e3c4629 Mon Sep 17 00:00:00 2001
From: Peter Munch <peterrmuench@gmail.com>
Date: Sat, 18 Mar 2023 15:47:49 +0100
Subject: [PATCH] SolverFMRES: refactor orthogonalization

---
 include/deal.II/lac/solver_gmres.h | 20 ++++++++++++++++----
 1 file changed, 16 insertions(+), 4 deletions(-)

diff --git a/include/deal.II/lac/solver_gmres.h b/include/deal.II/lac/solver_gmres.h
index c386a07553..54d2393af1 100644
--- a/include/deal.II/lac/solver_gmres.h
+++ b/include/deal.II/lac/solver_gmres.h
@@ -1738,6 +1738,8 @@ SolverFGMRES<VectorType>::solve(const MatrixType &        A,
   // matrix used for the orthogonalization process later
   H.reinit(basis_size + 1, basis_size);
 
+  Vector<double> h(basis_size + 1);
+
   // Vectors for projected system
   Vector<double> projected_rhs;
   Vector<double> y;
@@ -1773,10 +1775,20 @@ SolverFGMRES<VectorType>::solve(const MatrixType &        A,
           A.vmult(*aux, z[j]);
 
           // Gram-Schmidt
-          H(0, j) = *aux * v[0];
-          for (unsigned int i = 1; i <= j; ++i)
-            H(i, j) = aux->add_and_dot(-H(i - 1, j), v[i - 1], v[i]);
-          H(j + 1, j) = a = std::sqrt(aux->add_and_dot(-H(j, j), v[j], *aux));
+          bool         re_orthogonalize = false;
+          const double s =
+            internal::SolverGMRESImplementation::iterated_gram_schmidt<
+              VectorType>(SolverGMRES<VectorType>::AdditionalData::
+                            OrthogonalizationStrategy::modified_gram_schmidt,
+                          v,
+                          j + 1,
+                          0,
+                          *aux,
+                          h,
+                          re_orthogonalize);
+          for (unsigned int i = 0; i <= j; ++i)
+            H(i, j) = h(i);
+          H(j + 1, j) = a = s;
 
           // Compute projected solution
 
-- 
2.39.5