Move OrthogonalizationStrategy into new file

diff --git a/doc/news/changes/minor/20221208Munch b/doc/news/changes/minor/20221208Munch
index e90bef2f41e6f6ca5c73b897ee0c98c2e4614d58..4c7959e946f242aa1d62d62cf5e20ac58860b84b 100644 (file)
--- a/doc/news/changes/minor/20221208Munch
+++ b/doc/news/changes/minor/20221208Munch
@@ -1,6 +1,7 @@
-New: SolverGMRES now also supports classical Gram-Schmidt orthonormalization
+New: SolverGMRES and SolverFGMRES now also support classical 
+Gram-Schmidt orthonormalization
 alongside to the existing modified Gram-Schmidt algorithm. This
-allows to reduce the cost of vector operations in terms of
+algorithm allows to reduce the cost of vector operations in terms of
 communication latency and memory transfer significantly.
 <br>
 (Peter Munch, Martin Kronbichler, 2022/12/08)

diff --git a/include/deal.II/lac/orthogonalization.h b/include/deal.II/lac/orthogonalization.h
new file mode 100644 (file)
index 0000000..69ff82c
--- /dev/null
+++ b/include/deal.II/lac/orthogonalization.h
@@ -0,0 +1,50 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2023 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii_lac_orthogonalization_h
+#define dealii_lac_orthogonalization_h
+
+
+#include <deal.II/base/config.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace LinearAlgebra
+{
+  /**
+   * Supported orthogonalization strategies within SolverGMRES and
+   * SolverFGMRES.
+   */
+  enum class OrthogonalizationStrategy
+  {
+    /**
+     * Use modified Gram-Schmidt algorithm.
+     */
+    modified_gram_schmidt,
+    /**
+     * Use classical Gram-Schmidt algorithm. Since this approach works on
+     * multi-vectors with a single global reduction (of multiple elements), it
+     * is more efficient than the modified Gram-Schmidt algorithm. However, it
+     * is less stable in terms of roundoff error propagation, requiring
+     * additional re-orthogonalization steps more frequently.
+     */
+    classical_gram_schmidt
+  };
+} // namespace LinearAlgebra
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif

diff --git a/include/deal.II/lac/solver_gmres.h b/include/deal.II/lac/solver_gmres.h
index 54d2393af1e12825ab679838b7f201a40d7c3a56..bdeb3884c9af1630fb2cf7e2e3358a6341b8f962 100644 (file)
--- a/include/deal.II/lac/solver_gmres.h
+++ b/include/deal.II/lac/solver_gmres.h
@@ -28,6 +28,7 @@
 #include <deal.II/lac/full_matrix.h>
 #include <deal.II/lac/householder.h>
 #include <deal.II/lac/lapack_full_matrix.h>
+#include <deal.II/lac/orthogonalization.h>
 #include <deal.II/lac/solver.h>
 #include <deal.II/lac/solver_control.h>
 #include <deal.II/lac/vector.h>
@@ -198,21 +199,6 @@ public:
    */
   struct AdditionalData
   {
-    enum class OrthogonalizationStrategy
-    {
-      /**
-       * Use modified Gram-Schmidt algorithm.
-       */
-      modified_gram_schmidt,
-      /**
-       * Use classical Gram-Schmidt algorithm. Since this approach works on
-       * multi-vectors and performs a global reduction only once, it is
-       * more efficient than the modified Gram-Schmidt algorithm.
-       * However, it might be numerically unstable.
-       */
-      classical_gram_schmidt
-    };
-
     /**
      * Constructor. By default, set the number of temporary vectors to 30,
      * i.e. do a restart every 28 iterations. Also set preconditioning from
@@ -221,13 +207,14 @@ public:
      * reduced functionality to track information is disabled by default.
      */
     explicit AdditionalData(
-      const unsigned int              max_n_tmp_vectors          = 30,
-      const bool                      right_preconditioning      = false,
-      const bool                      use_default_residual       = true,
-      const bool                      force_re_orthogonalization = false,
-      const bool                      batched_mode               = false,
-      const OrthogonalizationStrategy orthogonalization_strategy =
-        OrthogonalizationStrategy::modified_gram_schmidt);
+      const unsigned int max_n_tmp_vectors          = 30,
+      const bool         right_preconditioning      = false,
+      const bool         use_default_residual       = true,
+      const bool         force_re_orthogonalization = false,
+      const bool         batched_mode               = false,
+      const LinearAlgebra::OrthogonalizationStrategy
+        orthogonalization_strategy =
+          LinearAlgebra::OrthogonalizationStrategy::modified_gram_schmidt);
 
     /**
      * Maximum number of temporary vectors. This parameter controls the size
@@ -270,7 +257,7 @@ public:
     /**
      * Strategy to orthogonalize vectors.
      */
-    OrthogonalizationStrategy orthogonalization_strategy;
+    LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy;
   };
 
   /**
@@ -520,14 +507,24 @@ public:
     /**
      * Constructor. By default, set the maximum basis size to 30.
      */
-    explicit AdditionalData(const unsigned int max_basis_size = 30)
+    explicit AdditionalData(
+      const unsigned int max_basis_size = 30,
+      const LinearAlgebra::OrthogonalizationStrategy
+        orthogonalization_strategy =
+          LinearAlgebra::OrthogonalizationStrategy::modified_gram_schmidt)
       : max_basis_size(max_basis_size)
+      , orthogonalization_strategy(orthogonalization_strategy)
     {}
 
     /**
      * Maximum basis size.
      */
     unsigned int max_basis_size;
+
+    /**
+     * Strategy to orthogonalize vectors.
+     */
+    LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy;
   };
 
   /**
@@ -659,12 +656,12 @@ namespace internal
 
 template <class VectorType>
 inline SolverGMRES<VectorType>::AdditionalData::AdditionalData(
-  const unsigned int              max_n_tmp_vectors,
-  const bool                      right_preconditioning,
-  const bool                      use_default_residual,
-  const bool                      force_re_orthogonalization,
-  const bool                      batched_mode,
-  const OrthogonalizationStrategy orthogonalization_strategy)
+  const unsigned int                             max_n_tmp_vectors,
+  const bool                                     right_preconditioning,
+  const bool                                     use_default_residual,
+  const bool                                     force_re_orthogonalization,
+  const bool                                     batched_mode,
+  const LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy)
   : max_n_tmp_vectors(max_n_tmp_vectors)
   , right_preconditioning(right_preconditioning)
   , use_default_residual(use_default_residual)
@@ -1133,8 +1130,7 @@ namespace internal
     template <class VectorType>
     inline double
     iterated_gram_schmidt(
-      const typename SolverGMRES<VectorType>::AdditionalData::
-        OrthogonalizationStrategy orthogonalization_strategy,
+      const LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy,
       const internal::SolverGMRESImplementation::TmpVectors<VectorType>
         &                                       orthogonal_vectors,
       const unsigned int                        dim,
@@ -1165,8 +1161,7 @@ namespace internal
           double norm_vv = 0.0;
 
           if (orthogonalization_strategy ==
-              SolverGMRES<VectorType>::AdditionalData::
-                OrthogonalizationStrategy::modified_gram_schmidt)
+              LinearAlgebra::OrthogonalizationStrategy::modified_gram_schmidt)
             {
               double htmp = vv * orthogonal_vectors[0];
               h(0) += htmp;
@@ -1182,8 +1177,8 @@ namespace internal
                 vv.add_and_dot(-htmp, orthogonal_vectors[dim - 1], vv));
             }
           else if (orthogonalization_strategy ==
-                   SolverGMRES<VectorType>::AdditionalData::
-                     OrthogonalizationStrategy::classical_gram_schmidt)
+                   LinearAlgebra::OrthogonalizationStrategy::
+                     classical_gram_schmidt)
             {
               Tvmult_add(dim, vv, orthogonal_vectors, h);
               norm_vv = substract_and_norm(dim, orthogonal_vectors, h, vv);
@@ -1778,8 +1773,7 @@ SolverFGMRES<VectorType>::solve(const MatrixType &        A,
           bool         re_orthogonalize = false;
           const double s =
             internal::SolverGMRESImplementation::iterated_gram_schmidt<
-              VectorType>(SolverGMRES<VectorType>::AdditionalData::
-                            OrthogonalizationStrategy::modified_gram_schmidt,
+              VectorType>(additional_data.orthogonalization_strategy,
                           v,
                           j + 1,
                           0,

diff --git a/tests/lac/solver.cc b/tests/lac/solver.cc
index 0ddfc6c3d344a3a64e29f54f3330685058c93194..ff81c85d8dd6373a3be13ed26e1011a6b3f5483b 100644 (file)
--- a/tests/lac/solver.cc
+++ b/tests/lac/solver.cc
@@ -101,8 +101,8 @@ main()
   SolverFIRE<>                  fire(control, mem);
 
   SolverGMRES<>::AdditionalData data3(8);
-  data3.orthogonalization_strategy = SolverGMRES<>::AdditionalData::
-    OrthogonalizationStrategy::classical_gram_schmidt;
+  data3.orthogonalization_strategy =
+    LinearAlgebra::OrthogonalizationStrategy::classical_gram_schmidt;
   SolverGMRES<> gmresclassical(control, mem, data3);
 
   for (unsigned int size = 4; size <= 30; size *= 3)

diff --git a/tests/lac/solver_gmres_01.cc b/tests/lac/solver_gmres_01.cc
index fb7ccad9662a2a7c9a9c131d57492a2e6926a886..a79f12e0f2088beb20debe20bc5e2dc8ac871b0e 100644 (file)
--- a/tests/lac/solver_gmres_01.cc
+++ b/tests/lac/solver_gmres_01.cc
@@ -14,7 +14,7 @@
 // ---------------------------------------------------------------------
 
 // Check the path of SolverGMRES with distributed vectors and
-// OrthogonalizationStrategy::classical_gram_schmidt.
+// LinearAlgebra::OrthogonalizationStrategy::classical_gram_schmidt.
 
 
 #include <deal.II/lac/diagonal_matrix.h>
@@ -108,8 +108,7 @@ main()
     SolverGMRES<LinearAlgebra::distributed::Vector<double>>::AdditionalData
       data3(8);
     data3.orthogonalization_strategy =
-      SolverGMRES<LinearAlgebra::distributed::Vector<double>>::AdditionalData::
-        OrthogonalizationStrategy::classical_gram_schmidt;
+      LinearAlgebra::OrthogonalizationStrategy::classical_gram_schmidt;
     SolverGMRES<LinearAlgebra::distributed::Vector<double>> solver(control,
                                                                    data3);
     solver.connect(&monitor_norm);
@@ -131,8 +130,7 @@ main()
     SolverGMRES<LinearAlgebra::distributed::BlockVector<double>>::AdditionalData
       data3(8);
     data3.orthogonalization_strategy =
-      SolverGMRES<LinearAlgebra::distributed::BlockVector<double>>::
-        AdditionalData::OrthogonalizationStrategy::classical_gram_schmidt;
+      LinearAlgebra::OrthogonalizationStrategy::classical_gram_schmidt;
     SolverGMRES<LinearAlgebra::distributed::BlockVector<double>> solver(control,
                                                                         data3);
     solver.connect(&monitor_norm_block);