From: Peter Munch <peterrmuench@gmail.com>
Date: Sat, 18 Mar 2023 14:47:49 +0000 (+0100)
Subject: SolverFMRES: refactor orthogonalization
X-Git-Tag: v9.5.0-rc1~416^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F14913%2Fhead;p=dealii.git

SolverFMRES: refactor orthogonalization
---

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