Fix bug in result of SolverCG when run in interleaved mode

diff --git a/include/deal.II/lac/solver_cg.h b/include/deal.II/lac/solver_cg.h
index 63902fd915fa89f61984eeaf2b4fa5d33fe90b26..d4d901860714c8b3cdff642437f997dee594d8dc 100644 (file)
--- a/include/deal.II/lac/solver_cg.h
+++ b/include/deal.II/lac/solver_cg.h
@@ -952,7 +952,7 @@ namespace internal
       std::enable_if_t<has_apply<PreconditionerType>, U>
       finalize_after_convergence(const unsigned int iteration_index)
       {
-        if (iteration_index % 2 == 1 || iteration_index == 2)
+        if (iteration_index % 2 == 1)
           this->x.add(this->alpha, this->p);
         else
           {
@@ -968,9 +968,9 @@ namespace internal
             // previous_beta has not been applied at this stage, allowing us
             // to recover the previous search direction
             const Number alpha_plus_previous_alpha_over_beta =
-              this->alpha + this->previous_alpha / this->previous_beta;
+              this->alpha + this->previous_alpha / this->beta;
             const Number previous_alpha_over_beta =
-              this->previous_alpha / this->previous_beta;
+              this->previous_alpha / this->beta;
 
             DEAL_II_OPENMP_SIMD_PRAGMA
             for (unsigned int j = 0; j < end_range; ++j)

diff --git a/tests/lac/solver_cg_interleave.cc b/tests/lac/solver_cg_interleave_01.cc
similarity index 99%
rename from tests/lac/solver_cg_interleave.cc
rename to tests/lac/solver_cg_interleave_01.cc
index c48efd677d6a5d7bec595dcf6e0c47bc20f13cba..bdd4d5c3501760f4320dec5b1c9f0523c55f2479 100644 (file)
--- a/tests/lac/solver_cg_interleave.cc
+++ b/tests/lac/solver_cg_interleave_01.cc
@@ -69,6 +69,7 @@ monitor_norm(const unsigned int iteration,
 }
 
 
+
 int
 main()
 {

diff --git a/tests/lac/solver_cg_interleave.output b/tests/lac/solver_cg_interleave_01.output
similarity index 100%
rename from tests/lac/solver_cg_interleave.output
rename to tests/lac/solver_cg_interleave_01.output

diff --git a/tests/lac/solver_cg_interleave_02.cc b/tests/lac/solver_cg_interleave_02.cc
new file mode 100644 (file)
index 0000000..04c2b54
--- /dev/null
+++ b/tests/lac/solver_cg_interleave_02.cc
@@ -0,0 +1,97 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+// Check the path of SolverCG that can apply loops to sub-ranges of the matrix
+// in terms of the final solution vector, terminating at arbitrary numbers of
+// iterations also when the solution has not converged yet, and make sure we
+// obtain the same result as without interleaving.
+
+
+#include <deal.II/lac/diagonal_matrix.h>
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+
+#include "../tests.h"
+
+
+struct MyDiagonalMatrix
+{
+  MyDiagonalMatrix(const LinearAlgebra::distributed::Vector<double> &diagonal)
+    : diagonal(diagonal)
+  {}
+
+  void
+  vmult(LinearAlgebra::distributed::Vector<double> &      dst,
+        const LinearAlgebra::distributed::Vector<double> &src) const
+  {
+    dst = src;
+    dst.scale(diagonal);
+  }
+
+  void
+  vmult(
+    LinearAlgebra::distributed::Vector<double> &                       dst,
+    const LinearAlgebra::distributed::Vector<double> &                 src,
+    const std::function<void(const unsigned int, const unsigned int)> &before,
+    const std::function<void(const unsigned int, const unsigned int)> &after)
+    const
+  {
+    before(0, dst.size());
+    vmult(dst, src);
+    after(0, dst.size());
+  }
+
+  const LinearAlgebra::distributed::Vector<double> &diagonal;
+};
+
+
+
+int
+main()
+{
+  initlog();
+
+  // Create diagonal matrix with entries between 1 and 15
+  DiagonalMatrix<LinearAlgebra::distributed::Vector<double>> unit_matrix;
+  unit_matrix.get_vector().reinit(15);
+  unit_matrix.get_vector() = 1.0;
+
+  LinearAlgebra::distributed::Vector<double> matrix_entries(unit_matrix.m());
+  for (unsigned int i = 0; i < unit_matrix.m(); ++i)
+    matrix_entries(i) = i + 1;
+  MyDiagonalMatrix matrix(matrix_entries);
+
+  LinearAlgebra::distributed::Vector<double> rhs(unit_matrix.m()),
+    sol(unit_matrix.m());
+  rhs = 1.;
+
+  for (unsigned int n_iter = 1; n_iter < 12; ++n_iter)
+    {
+      deallog << "Test solution accuracy with " << n_iter << " iterations"
+              << std::endl;
+      deallog << "Solve with PreconditionIdentity: " << std::endl;
+      IterationNumberControl control(n_iter, 1e-12);
+      SolverCG<LinearAlgebra::distributed::Vector<double>> solver(control);
+      sol = 0;
+      solver.solve(matrix, sol, rhs, PreconditionIdentity());
+      sol.print(deallog.get_file_stream());
+
+      deallog << "Solve with diagonal preconditioner: " << std::endl;
+      sol = 0;
+      solver.solve(matrix, sol, rhs, unit_matrix);
+      sol.print(deallog.get_file_stream());
+    }
+}

diff --git a/tests/lac/solver_cg_interleave_02.output b/tests/lac/solver_cg_interleave_02.output
new file mode 100644 (file)
index 0000000..ef5e212
--- /dev/null
+++ b/tests/lac/solver_cg_interleave_02.output
@@ -0,0 +1,166 @@
+
+DEAL::Test solution accuracy with 1 iterations
+DEAL::Solve with PreconditionIdentity: 
+DEAL:cg::Starting value 3.87298
+DEAL:cg::Convergence step 1 value 2.09165
+Process #0
+Local range: [0, 15), global size: 15
+Vector data:
+1.250e-01 1.250e-01 1.250e-01 1.250e-01 1.250e-01 1.250e-01 1.250e-01 1.250e-01 1.250e-01 1.250e-01 1.250e-01 1.250e-01 1.250e-01 1.250e-01 1.250e-01 
+DEAL::Solve with diagonal preconditioner: 
+DEAL:cg::Starting value 3.87298
+DEAL:cg::Convergence step 1 value 2.09165
+Process #0
+Local range: [0, 15), global size: 15
+Vector data:
+1.250e-01 1.250e-01 1.250e-01 1.250e-01 1.250e-01 1.250e-01 1.250e-01 1.250e-01 1.250e-01 1.250e-01 1.250e-01 1.250e-01 1.250e-01 1.250e-01 1.250e-01 
+DEAL::Test solution accuracy with 2 iterations
+DEAL::Solve with PreconditionIdentity: 
+DEAL:cg::Starting value 3.87298
+DEAL:cg::Convergence step 2 value 1.41681
+Process #0
+Local range: [0, 15), global size: 15
+Vector data:
+3.309e-01 3.088e-01 2.868e-01 2.647e-01 2.426e-01 2.206e-01 1.985e-01 1.765e-01 1.544e-01 1.324e-01 1.103e-01 8.824e-02 6.618e-02 4.412e-02 2.206e-02 
+DEAL::Solve with diagonal preconditioner: 
+DEAL:cg::Starting value 3.87298
+DEAL:cg::Convergence step 2 value 1.41681
+Process #0
+Local range: [0, 15), global size: 15
+Vector data:
+3.309e-01 3.088e-01 2.868e-01 2.647e-01 2.426e-01 2.206e-01 1.985e-01 1.765e-01 1.544e-01 1.324e-01 1.103e-01 8.824e-02 6.618e-02 4.412e-02 2.206e-02 
+DEAL::Test solution accuracy with 3 iterations
+DEAL::Solve with PreconditionIdentity: 
+DEAL:cg::Starting value 3.87298
+DEAL:cg::Convergence step 3 value 0.977692
+Process #0
+Local range: [0, 15), global size: 15
+Vector data:
+5.539e-01 4.681e-01 3.905e-01 3.211e-01 2.598e-01 2.067e-01 1.618e-01 1.250e-01 9.641e-02 7.598e-02 6.373e-02 5.964e-02 6.373e-02 7.598e-02 9.641e-02 
+DEAL::Solve with diagonal preconditioner: 
+DEAL:cg::Starting value 3.87298
+DEAL:cg::Convergence step 3 value 0.977692
+Process #0
+Local range: [0, 15), global size: 15
+Vector data:
+5.539e-01 4.681e-01 3.905e-01 3.211e-01 2.598e-01 2.067e-01 1.618e-01 1.250e-01 9.641e-02 7.598e-02 6.373e-02 5.964e-02 6.373e-02 7.598e-02 9.641e-02 
+DEAL::Test solution accuracy with 4 iterations
+DEAL::Solve with PreconditionIdentity: 
+DEAL:cg::Starting value 3.87298
+DEAL:cg::Convergence step 4 value 0.656069
+Process #0
+Local range: [0, 15), global size: 15
+Vector data:
+7.417e-01 5.553e-01 4.081e-01 2.954e-01 2.128e-01 1.559e-01 1.200e-01 1.006e-01 9.331e-02 9.352e-02 9.675e-02 9.847e-02 9.417e-02 7.933e-02 4.945e-02 
+DEAL::Solve with diagonal preconditioner: 
+DEAL:cg::Starting value 3.87298
+DEAL:cg::Convergence step 4 value 0.656069
+Process #0
+Local range: [0, 15), global size: 15
+Vector data:
+7.417e-01 5.553e-01 4.081e-01 2.954e-01 2.128e-01 1.559e-01 1.200e-01 1.006e-01 9.331e-02 9.352e-02 9.675e-02 9.847e-02 9.417e-02 7.933e-02 4.945e-02 
+DEAL::Test solution accuracy with 5 iterations
+DEAL::Solve with PreconditionIdentity: 
+DEAL:cg::Starting value 3.87298
+DEAL:cg::Convergence step 5 value 0.419623
+Process #0
+Local range: [0, 15), global size: 15
+Vector data:
+8.709e-01 5.738e-01 3.754e-01 2.514e-01 1.806e-01 1.452e-01 1.304e-01 1.250e-01 1.208e-01 1.129e-01 9.972e-02 8.286e-02 6.721e-02 6.089e-02 7.528e-02 
+DEAL::Solve with diagonal preconditioner: 
+DEAL:cg::Starting value 3.87298
+DEAL:cg::Convergence step 5 value 0.419623
+Process #0
+Local range: [0, 15), global size: 15
+Vector data:
+8.709e-01 5.738e-01 3.754e-01 2.514e-01 1.806e-01 1.452e-01 1.304e-01 1.250e-01 1.208e-01 1.129e-01 9.972e-02 8.286e-02 6.721e-02 6.089e-02 7.528e-02 
+DEAL::Test solution accuracy with 6 iterations
+DEAL::Solve with PreconditionIdentity: 
+DEAL:cg::Starting value 3.87298
+DEAL:cg::Convergence step 6 value 0.252694
+Process #0
+Local range: [0, 15), global size: 15
+Vector data:
+9.447e-01 5.553e-01 3.404e-01 2.330e-01 1.848e-01 1.634e-01 1.498e-01 1.347e-01 1.165e-01 9.807e-02 8.398e-02 7.766e-02 7.856e-02 7.933e-02 6.298e-02 
+DEAL::Solve with diagonal preconditioner: 
+DEAL:cg::Starting value 3.87298
+DEAL:cg::Convergence step 6 value 0.252694
+Process #0
+Local range: [0, 15), global size: 15
+Vector data:
+9.447e-01 5.553e-01 3.404e-01 2.330e-01 1.848e-01 1.634e-01 1.498e-01 1.347e-01 1.165e-01 9.807e-02 8.398e-02 7.766e-02 7.856e-02 7.933e-02 6.298e-02 
+DEAL::Test solution accuracy with 7 iterations
+DEAL::Solve with PreconditionIdentity: 
+DEAL:cg::Starting value 3.87298
+DEAL:cg::Convergence step 7 value 0.141859
+Process #0
+Local range: [0, 15), global size: 15
+Vector data:
+9.799e-01 5.302e-01 3.246e-01 2.380e-01 1.991e-01 1.731e-01 1.484e-01 1.250e-01 1.068e-01 9.613e-02 9.133e-02 8.733e-02 7.895e-02 6.712e-02 6.801e-02 
+DEAL::Solve with diagonal preconditioner: 
+DEAL:cg::Starting value 3.87298
+DEAL:cg::Convergence step 7 value 0.141859
+Process #0
+Local range: [0, 15), global size: 15
+Vector data:
+9.799e-01 5.302e-01 3.246e-01 2.380e-01 1.991e-01 1.731e-01 1.484e-01 1.250e-01 1.068e-01 9.613e-02 9.133e-02 8.733e-02 7.895e-02 6.712e-02 6.801e-02 
+DEAL::Test solution accuracy with 8 iterations
+DEAL::Solve with PreconditionIdentity: 
+DEAL:cg::Starting value 3.87298
+DEAL:cg::Convergence step 8 value 0.0735126
+Process #0
+Local range: [0, 15), global size: 15
+Vector data:
+9.939e-01 5.127e-01 3.240e-01 2.474e-01 2.042e-01 1.698e-01 1.417e-01 1.221e-01 1.102e-01 1.019e-01 9.281e-02 8.245e-02 7.477e-02 7.324e-02 6.626e-02 
+DEAL::Solve with diagonal preconditioner: 
+DEAL:cg::Starting value 3.87298
+DEAL:cg::Convergence step 8 value 0.0735126
+Process #0
+Local range: [0, 15), global size: 15
+Vector data:
+9.939e-01 5.127e-01 3.240e-01 2.474e-01 2.042e-01 1.698e-01 1.417e-01 1.221e-01 1.102e-01 1.019e-01 9.281e-02 8.245e-02 7.477e-02 7.324e-02 6.626e-02 
+DEAL::Test solution accuracy with 9 iterations
+DEAL::Solve with PreconditionIdentity: 
+DEAL:cg::Starting value 3.87298
+DEAL:cg::Convergence step 9 value 0.0347680
+Process #0
+Local range: [0, 15), global size: 15
+Vector data:
+9.985e-01 5.042e-01 3.283e-01 2.514e-01 2.022e-01 1.660e-01 1.412e-01 1.250e-01 1.124e-01 1.004e-01 8.990e-02 8.285e-02 7.809e-02 7.083e-02 6.677e-02 
+DEAL::Solve with diagonal preconditioner: 
+DEAL:cg::Starting value 3.87298
+DEAL:cg::Convergence step 9 value 0.0347680
+Process #0
+Local range: [0, 15), global size: 15
+Vector data:
+9.985e-01 5.042e-01 3.283e-01 2.514e-01 2.022e-01 1.660e-01 1.412e-01 1.250e-01 1.124e-01 1.004e-01 8.990e-02 8.285e-02 7.809e-02 7.083e-02 6.677e-02 
+DEAL::Test solution accuracy with 10 iterations
+DEAL::Solve with PreconditionIdentity: 
+DEAL:cg::Starting value 3.87298
+DEAL:cg::Convergence step 10 value 0.0147898
+Process #0
+Local range: [0, 15), global size: 15
+Vector data:
+9.997e-01 5.010e-01 3.315e-01 2.515e-01 2.001e-01 1.658e-01 1.429e-01 1.256e-01 1.111e-01 9.947e-02 9.096e-02 8.382e-02 7.650e-02 7.158e-02 6.665e-02 
+DEAL::Solve with diagonal preconditioner: 
+DEAL:cg::Starting value 3.87298
+DEAL:cg::Convergence step 10 value 0.0147898
+Process #0
+Local range: [0, 15), global size: 15
+Vector data:
+9.997e-01 5.010e-01 3.315e-01 2.515e-01 2.001e-01 1.658e-01 1.429e-01 1.256e-01 1.111e-01 9.947e-02 9.096e-02 8.382e-02 7.650e-02 7.158e-02 6.665e-02 
+DEAL::Test solution accuracy with 11 iterations
+DEAL::Solve with PreconditionIdentity: 
+DEAL:cg::Starting value 3.87298
+DEAL:cg::Convergence step 11 value 0.00554307
+Process #0
+Local range: [0, 15), global size: 15
+Vector data:
+1.000e+00 5.002e-01 3.329e-01 2.506e-01 1.997e-01 1.665e-01 1.431e-01 1.250e-01 1.109e-01 1.001e-01 9.106e-02 8.313e-02 7.703e-02 7.140e-02 6.667e-02 
+DEAL::Solve with diagonal preconditioner: 
+DEAL:cg::Starting value 3.87298
+DEAL:cg::Convergence step 11 value 0.00554307
+Process #0
+Local range: [0, 15), global size: 15
+Vector data:
+1.000e+00 5.002e-01 3.329e-01 2.506e-01 1.997e-01 1.665e-01 1.431e-01 1.250e-01 1.109e-01 1.001e-01 9.106e-02 8.313e-02 7.703e-02 7.140e-02 6.667e-02