SolverGMRES: Switch default orthogonalization strategy to delayed CGS2

diff --git a/doc/news/changes/incompatibilities/20240312Kronbichler b/doc/news/changes/incompatibilities/20240312Kronbichler
index bcc4ad89e6bf4a317ce3ef35e702636750c643d3..05c708e1961de2722b3f926d9b4b3fd0b92adf52 100644 (file)
--- a/doc/news/changes/incompatibilities/20240312Kronbichler
+++ b/doc/news/changes/incompatibilities/20240312Kronbichler
@@ -6,6 +6,9 @@ value of the basis size is now 30, compared to 28 used before. The old
 variable SolverGMRES::AdditionalData::max_n_tmp_vectors is still available,
 but whenever SolverGMRES::AdditionalData::max_basis_size is set to a non-zero
 value (including the value set by the default constructor), the latter takes
-precedence.
+precedence. Furthermore, the default algorithm has been changed from
+the modified Gram-Schmidt algorithm to the classical Gram-Schmidt algorithm.
+The latter uses unconditional reorthogonalization delayed by one step,
+following the algorithm described in @cite Bielich2022.
 <br>
-(Martin Kronbichler, 2024/04/12)
+(Martin Kronbichler, 2024/03/12)

diff --git a/doc/news/changes/minor/20240312Kronbichler b/doc/news/changes/minor/20240312Kronbichler
index 77a5da56df3f0635b3c7dbd79db1966e629b4bab..10d90b7d2b25c78a88a428f55c60a84fef392fb5 100644 (file)
--- a/doc/news/changes/minor/20240312Kronbichler
+++ b/doc/news/changes/minor/20240312Kronbichler
@@ -6,4 +6,4 @@ SolverGMRES. Since the Arnoldi process is sensitive to roundoff errors, this
 change might slightly affect iteration counts (often giving slightly better
 results).
 <br>
-(Martin Kronbichler, 2024/04/12)
+(Martin Kronbichler, 2024/03/12)

diff --git a/doc/news/changes/minor/20240319Kronbichler b/doc/news/changes/minor/20240319Kronbichler
new file mode 100644 (file)
index 0000000..e858a30
--- /dev/null
+++ b/doc/news/changes/minor/20240319Kronbichler
@@ -0,0 +1,8 @@
+New: SolverGMRES and SolverFGMRES can now use an additional orthogonalization
+strategy, controlled by
+LinearAlgebra::OrthogonalizationStrategy::delayed_classical_gram_schmidt. This
+implements the classical Gram-Schmidt method with delayed reorthogonalization,
+a low-synchronization algorithm (performing one global reduction per GMRES
+iteration for deal.II's own vectors) that has excellent stability properties.
+<br>
+(Martin Kronbichler, 2024/03/19)

diff --git a/include/deal.II/lac/solver_gmres.h b/include/deal.II/lac/solver_gmres.h
index 667bb99a033227410eb1bc77485f131985fa9c0c..7229f3ce9cc12dcf0e3480e9fadec720f448b87e 100644 (file)
--- a/include/deal.II/lac/solver_gmres.h
+++ b/include/deal.II/lac/solver_gmres.h
@@ -358,17 +358,19 @@ public:
      * to the default residual, and re-orthogonalization only if
      * necessary. Also, the batched mode with reduced functionality to track
      * information is disabled by default. Finally, the default
-     * orthogonalization algorithm is the modified Gram-Schmidt method.
+     * orthogonalization algorithm is the classical Gram-Schmidt method with
+     * delayed reorthogonalization, which combines stability with fast
+     * execution, especially in parallel.
      */
-    explicit AdditionalData(
-      const unsigned int max_basis_size             = 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);
+    explicit AdditionalData(const unsigned int max_basis_size        = 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::
+                                  delayed_classical_gram_schmidt);
 
     /**
      * Maximum number of temporary vectors. Together with max_basis_size, this

diff --git a/tests/lac/gmres_reorthogonalize_01.cc b/tests/lac/gmres_reorthogonalize_01.cc
index 78a5b995a0579bddc71e280f6932437c407190f9..204feebe090819c94c031523f7df7214e25b64b5 100644 (file)
--- a/tests/lac/gmres_reorthogonalize_01.cc
+++ b/tests/lac/gmres_reorthogonalize_01.cc
@@ -15,7 +15,7 @@
 
 // tests that GMRES builds an orthonormal basis properly for a few difficult
 // test matrices. In particular, this test monitors when re-orthogonalization
-// kicks in.
+// kicks in for the modified Gram-Schmidt algorithm.
 
 #include <deal.II/lac/full_matrix.h>
 #include <deal.II/lac/precondition.h>
@@ -70,6 +70,8 @@ test(unsigned int variant, unsigned int min_convergence_steps)
   SolverControl control(1000, 1e2 * std::numeric_limits<number>::epsilon());
   typename SolverGMRES<Vector<number>>::AdditionalData data;
   data.max_basis_size = 80;
+  data.orthogonalization_strategy =
+    LinearAlgebra::OrthogonalizationStrategy::modified_gram_schmidt;
 
   SolverGMRES<Vector<number>> solver(control, data);
   auto print_re_orthogonalization = [](int accumulated_iterations) {

diff --git a/tests/lac/gmres_reorthogonalize_02.cc b/tests/lac/gmres_reorthogonalize_02.cc
index 572485a7e67f9bc259488f252e0391ad5d1b842c..576b9765d5ce7184361d986c54eeca658b4e69c2 100644 (file)
--- a/tests/lac/gmres_reorthogonalize_02.cc
+++ b/tests/lac/gmres_reorthogonalize_02.cc
@@ -44,6 +44,8 @@ test()
   SolverControl control(1000, 1e3 * std::numeric_limits<number>::epsilon());
   typename SolverGMRES<Vector<number>>::AdditionalData data;
   data.max_basis_size = 200;
+  data.orthogonalization_strategy =
+    LinearAlgebra::OrthogonalizationStrategy::modified_gram_schmidt;
 
   SolverGMRES<Vector<number>> solver(control, data);
   auto print_re_orthogonalization = [](int accumulated_iterations) {

diff --git a/tests/lac/gmres_reorthogonalize_03.cc b/tests/lac/gmres_reorthogonalize_03.cc
index bccc1a35244d671bb9c02281498942e95924ced7..9f441b638506135b40b6f03e559c0a4a23f82201 100644 (file)
--- a/tests/lac/gmres_reorthogonalize_03.cc
+++ b/tests/lac/gmres_reorthogonalize_03.cc
@@ -69,6 +69,8 @@ test(unsigned int variant)
   typename SolverGMRES<Vector<number>>::AdditionalData data;
   data.max_basis_size             = 80;
   data.force_re_orthogonalization = true;
+  data.orthogonalization_strategy =
+    LinearAlgebra::OrthogonalizationStrategy::modified_gram_schmidt;
 
   SolverGMRES<Vector<number>> solver(control, data);
   auto print_re_orthogonalization = [](int accumulated_iterations) {

diff --git a/tests/lac/gmres_reorthogonalize_04.cc b/tests/lac/gmres_reorthogonalize_04.cc
index fd41939e9244ce7115882c5611ba01831610149a..883f84124ecd3140c7e047f04b36a5df25cd90d2 100644 (file)
--- a/tests/lac/gmres_reorthogonalize_04.cc
+++ b/tests/lac/gmres_reorthogonalize_04.cc
@@ -74,6 +74,8 @@ test(const unsigned int n_expected_steps)
   SolverControl control(1000, 1e3 * std::numeric_limits<number>::epsilon());
   typename SolverGMRES<Vector<number>>::AdditionalData data;
   data.max_basis_size = 200;
+  data.orthogonalization_strategy =
+    LinearAlgebra::OrthogonalizationStrategy::modified_gram_schmidt;
 
   SolverGMRES<Vector<number>> solver(control, data);
   auto print_re_orthogonalization = [](int accumulated_iterations) {

diff --git a/tests/lac/gmres_reorthogonalize_05.cc b/tests/lac/gmres_reorthogonalize_05.cc
index 37774657d89112e14e20c6b308c8aa565df6d52a..a251c0d7da0c1b2dd00d1b9ac9f2c12f0f0177fc 100644 (file)
--- a/tests/lac/gmres_reorthogonalize_05.cc
+++ b/tests/lac/gmres_reorthogonalize_05.cc
@@ -72,6 +72,8 @@ test()
   typename SolverGMRES<Vector<number>>::AdditionalData data;
   data.max_basis_size             = 200;
   data.force_re_orthogonalization = true;
+  data.orthogonalization_strategy =
+    LinearAlgebra::OrthogonalizationStrategy::modified_gram_schmidt;
 
   SolverGMRES<Vector<number>> solver(control, data);
   auto print_re_orthogonalization = [](int accumulated_iterations) {