Add two new test cases for interleaved solver

diff --git a/tests/lac/solver_cg_interleave.cc b/tests/lac/solver_cg_interleave.cc
new file mode 100644 (file)
index 0000000..5047c73
--- /dev/null
+++ b/tests/lac/solver_cg_interleave.cc
@@ -0,0 +1,100 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 1998 - 2020 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. That feature is used for matrix-free loops with increased data
+// locality.
+
+
+#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;
+};
+
+
+
+SolverControl::State
+monitor_norm(const unsigned int iteration,
+             const double       check_value,
+             const LinearAlgebra::distributed::Vector<double> &)
+{
+  deallog << "   CG estimated residual at iteration " << iteration << ": "
+          << check_value << std::endl;
+  return SolverControl::success;
+}
+
+
+int
+main()
+{
+  initlog();
+
+  // Create diagonal matrix with entries between 1 and 30
+  DiagonalMatrix<LinearAlgebra::distributed::Vector<double>> unit_matrix;
+  unit_matrix.get_vector().reinit(30);
+  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.;
+
+  deallog << "Solve with PreconditionIdentity: " << std::endl;
+  SolverControl                                        control(30, 1e-4);
+  SolverCG<LinearAlgebra::distributed::Vector<double>> solver(control);
+  solver.connect(&monitor_norm);
+  solver.solve(matrix, sol, rhs, PreconditionIdentity());
+
+  deallog << "Solve with diagonal preconditioner: " << std::endl;
+  sol = 0;
+  solver.solve(matrix, sol, rhs, unit_matrix);
+}

diff --git a/tests/lac/solver_cg_interleave.output b/tests/lac/solver_cg_interleave.output
new file mode 100644 (file)
index 0000000..2485887
--- /dev/null
+++ b/tests/lac/solver_cg_interleave.output
@@ -0,0 +1,53 @@
+
+DEAL::Solve with PreconditionIdentity: 
+DEAL:cg::Starting value 5.47723
+DEAL:cg::   CG estimated residual at iteration 0: 5.47723
+DEAL:cg::   CG estimated residual at iteration 1: 3.05857
+DEAL:cg::   CG estimated residual at iteration 2: 2.21614
+DEAL:cg::   CG estimated residual at iteration 3: 1.69418
+DEAL:cg::   CG estimated residual at iteration 4: 1.30657
+DEAL:cg::   CG estimated residual at iteration 5: 0.998837
+DEAL:cg::   CG estimated residual at iteration 6: 0.750194
+DEAL:cg::   CG estimated residual at iteration 7: 0.550634
+DEAL:cg::   CG estimated residual at iteration 8: 0.393553
+DEAL:cg::   CG estimated residual at iteration 9: 0.273167
+DEAL:cg::   CG estimated residual at iteration 10: 0.183730
+DEAL:cg::   CG estimated residual at iteration 11: 0.119512
+DEAL:cg::   CG estimated residual at iteration 12: 0.0750441
+DEAL:cg::   CG estimated residual at iteration 13: 0.0454041
+DEAL:cg::   CG estimated residual at iteration 14: 0.0264187
+DEAL:cg::   CG estimated residual at iteration 15: 0.0147526
+DEAL:cg::   CG estimated residual at iteration 16: 0.00788820
+DEAL:cg::   CG estimated residual at iteration 17: 0.00402832
+DEAL:cg::   CG estimated residual at iteration 18: 0.00195897
+DEAL:cg::   CG estimated residual at iteration 19: 0.000904053
+DEAL:cg::   CG estimated residual at iteration 20: 0.000394320
+DEAL:cg::   CG estimated residual at iteration 21: 0.000161750
+DEAL:cg::Convergence step 22 value 6.20175e-05
+DEAL:cg::   CG estimated residual at iteration 22: 6.20175e-05
+DEAL::Solve with diagonal preconditioner: 
+DEAL:cg::Starting value 5.47723
+DEAL:cg::   CG estimated residual at iteration 0: 5.47723
+DEAL:cg::   CG estimated residual at iteration 1: 3.05857
+DEAL:cg::   CG estimated residual at iteration 2: 2.21614
+DEAL:cg::   CG estimated residual at iteration 3: 1.69418
+DEAL:cg::   CG estimated residual at iteration 4: 1.30657
+DEAL:cg::   CG estimated residual at iteration 5: 0.998837
+DEAL:cg::   CG estimated residual at iteration 6: 0.750194
+DEAL:cg::   CG estimated residual at iteration 7: 0.550634
+DEAL:cg::   CG estimated residual at iteration 8: 0.393553
+DEAL:cg::   CG estimated residual at iteration 9: 0.273167
+DEAL:cg::   CG estimated residual at iteration 10: 0.183730
+DEAL:cg::   CG estimated residual at iteration 11: 0.119512
+DEAL:cg::   CG estimated residual at iteration 12: 0.0750441
+DEAL:cg::   CG estimated residual at iteration 13: 0.0454041
+DEAL:cg::   CG estimated residual at iteration 14: 0.0264187
+DEAL:cg::   CG estimated residual at iteration 15: 0.0147526
+DEAL:cg::   CG estimated residual at iteration 16: 0.00788820
+DEAL:cg::   CG estimated residual at iteration 17: 0.00402832
+DEAL:cg::   CG estimated residual at iteration 18: 0.00195897
+DEAL:cg::   CG estimated residual at iteration 19: 0.000904053
+DEAL:cg::   CG estimated residual at iteration 20: 0.000394320
+DEAL:cg::   CG estimated residual at iteration 21: 0.000161750
+DEAL:cg::Convergence step 22 value 6.20175e-05
+DEAL:cg::   CG estimated residual at iteration 22: 6.20175e-05

diff --git a/tests/matrix_free/solver_cg_interleave.cc b/tests/matrix_free/solver_cg_interleave.cc
new file mode 100644 (file)
index 0000000..0063740
--- /dev/null
+++ b/tests/matrix_free/solver_cg_interleave.cc
@@ -0,0 +1,234 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2022 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// this tests the correctness of using the SolverCG class with interleaved
+// vector operations, making use of the data locality options available by
+// MatrixFree loops.
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/manifold_lib.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/diagonal_matrix.h>
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/solver_cg.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+
+#include <iostream>
+
+#include "../tests.h"
+
+#include "matrix_vector_mf.h"
+
+
+
+template <int dim, typename number = double>
+class HelmholtzOperator : public Subscriptor
+{
+public:
+  using value_type = number;
+  using VectorType = LinearAlgebra::distributed::Vector<number>;
+
+  HelmholtzOperator(const MatrixFree<dim, number> &matrix_free)
+    : n_calls_vmult(0)
+    , data(matrix_free)
+  {}
+
+  ~HelmholtzOperator()
+  {
+    if (n_calls_vmult > 0)
+      deallog << "Number of calls to special vmult: " << n_calls_vmult
+              << std::endl;
+  }
+
+  void
+  initialize(const Mapping<dim> &mapping, DoFHandler<dim> &dof_handler)
+  {
+    n_calls_vmult = 0;
+  }
+
+  void
+  vmult(VectorType &dst, const VectorType &src) const
+  {
+    const std::function<
+      void(const MatrixFree<dim, typename VectorType::value_type> &,
+           VectorType &,
+           const VectorType &,
+           const std::pair<unsigned int, unsigned int> &)>
+      function = helmholtz_operator<dim, -1, VectorType, 0>;
+    data.cell_loop(function, dst, src, true);
+    for (unsigned int i : data.get_constrained_dofs())
+      dst.local_element(i) = src.local_element(i);
+  }
+
+  void
+  vmult(VectorType &      dst,
+        const VectorType &src,
+        const std::function<void(const unsigned int, const unsigned int)>
+          &operation_before_loop,
+        const std::function<void(const unsigned int, const unsigned int)>
+          &operation_after_loop) const
+  {
+    ++n_calls_vmult;
+    const std::function<
+      void(const MatrixFree<dim, typename VectorType::value_type> &,
+           VectorType &,
+           const VectorType &,
+           const std::pair<unsigned int, unsigned int> &)>
+      function = helmholtz_operator<dim, -1, VectorType, 0>;
+    data.cell_loop(
+      function, dst, src, operation_before_loop, operation_after_loop);
+    for (unsigned int i : data.get_constrained_dofs())
+      dst.local_element(i) = src.local_element(i);
+  }
+
+private:
+  mutable unsigned int    n_calls_vmult;
+  MatrixFree<dim, number> data;
+};
+
+
+template <typename Number>
+struct MyDiagonalMatrix
+{
+  MyDiagonalMatrix(const LinearAlgebra::distributed::Vector<Number> &vec)
+    : vec(vec)
+  {}
+
+  void
+  vmult(LinearAlgebra::distributed::Vector<Number> &      dst,
+        const LinearAlgebra::distributed::Vector<Number> &src) const
+  {
+    dst = src;
+    dst.scale(vec);
+  }
+
+  const LinearAlgebra::distributed::Vector<Number> &vec;
+};
+
+
+
+template <int dim, typename number>
+void
+test(const unsigned int fe_degree)
+{
+  using VectorType = LinearAlgebra::distributed::Vector<number>;
+
+  parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
+  GridGenerator::hyper_cube(tria);
+  tria.refine_global(6 - dim);
+
+  MappingQ1<dim>  mapping;
+  FE_Q<dim>       fe(fe_degree);
+  DoFHandler<dim> dof(tria);
+  dof.distribute_dofs(fe);
+
+  IndexSet                  owned_set = dof.locally_owned_dofs();
+  IndexSet                  relevant_set;
+  AffineConstraints<double> constraints(relevant_set);
+  typename MatrixFree<dim, number>::AdditionalData addit_data;
+  addit_data.tasks_parallel_scheme =
+    MatrixFree<dim, number>::AdditionalData::none;
+
+  {
+    DoFTools::extract_locally_relevant_dofs(dof, relevant_set);
+    constraints.close();
+
+    DoFRenumbering::matrix_free_data_locality(dof, constraints, addit_data);
+  }
+
+  constraints.clear();
+  DoFTools::extract_locally_relevant_dofs(dof, relevant_set);
+  constraints.close();
+
+  const QGauss<1>         quadrature(dof.get_fe().degree + 1);
+  MatrixFree<dim, number> mf_data;
+  mf_data.reinit(mapping, dof, constraints, quadrature, addit_data);
+
+  MatrixFreeTest<dim, -1, number, VectorType> mf(mf_data);
+  VectorType                                  rhs, sol;
+  mf_data.initialize_dof_vector(rhs);
+  mf_data.initialize_dof_vector(sol);
+  rhs = 1.;
+  DiagonalMatrix<VectorType> preconditioner;
+  mf_data.initialize_dof_vector(preconditioner.get_vector());
+  mf.vmult(preconditioner.get_vector(), rhs);
+  for (number &a : preconditioner.get_vector())
+    if (a != 0.)
+      a = 1. / a;
+    else
+      a = 0.;
+
+  // Step 1: solve with CG solver for a matrix that does not support the
+  // interleaving operation
+  {
+    SolverControl        control(200, 1e-2 * rhs.l2_norm());
+    SolverCG<VectorType> solver(control);
+    solver.solve(mf, sol, rhs, preconditioner);
+    deallog << "Norm of the solution: " << sol.l2_norm() << std::endl;
+  }
+
+  // Step 2: solve with CG solver for a matrix that does support the
+  // interleaving operation
+  {
+    sol = 0;
+    HelmholtzOperator<dim, number> matrix(mf_data);
+    SolverControl                  control(200, 1e-2 * rhs.l2_norm());
+    SolverCG<VectorType>           solver(control);
+    solver.solve(matrix, sol, rhs, preconditioner);
+    deallog << "Norm of the solution: " << sol.l2_norm() << std::endl;
+  }
+
+  // Step 3: solve with CG solver for a matrix that does support but a
+  // preconditioner that does not support the interleaving operation
+  {
+    sol = 0;
+    HelmholtzOperator<dim, number> matrix(mf_data);
+    MyDiagonalMatrix               simple_diagonal(preconditioner.get_vector());
+    SolverControl                  control(200, 1e-2 * rhs.l2_norm());
+    SolverCG<VectorType>           solver(control);
+    solver.solve(matrix, sol, rhs, simple_diagonal);
+    deallog << "Norm of the solution: " << sol.l2_norm() << std::endl;
+  }
+}
+
+
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+  mpi_initlog();
+
+  test<2, double>(3);
+  test<3, double>(4);
+  test<3, double>(3);
+}

diff --git a/tests/matrix_free/solver_cg_interleave.with_p4est=true.mpirun=3.output b/tests/matrix_free/solver_cg_interleave.with_p4est=true.mpirun=3.output
new file mode 100644 (file)
index 0000000..7c60fc4
--- /dev/null
+++ b/tests/matrix_free/solver_cg_interleave.with_p4est=true.mpirun=3.output
@@ -0,0 +1,31 @@
+
+DEAL:cg::Starting value 49.0000
+DEAL:cg::Convergence step 51 value 0.429721
+DEAL::Norm of the solution: 11776.3
+DEAL:cg::Starting value 49.0000
+DEAL:cg::Convergence step 51 value 0.429721
+DEAL::Norm of the solution: 11776.3
+DEAL::Number of calls to special vmult: 51
+DEAL:cg::Starting value 49.0000
+DEAL:cg::Convergence step 51 value 0.429721
+DEAL::Norm of the solution: 11776.3
+DEAL:cg::Starting value 189.571
+DEAL:cg::Convergence step 61 value 1.71195
+DEAL::Norm of the solution: 684335.
+DEAL:cg::Starting value 189.571
+DEAL:cg::Convergence step 61 value 1.71197
+DEAL::Norm of the solution: 684335.
+DEAL::Number of calls to special vmult: 61
+DEAL:cg::Starting value 189.571
+DEAL:cg::Convergence step 61 value 1.71195
+DEAL::Norm of the solution: 684335.
+DEAL:cg::Starting value 125.000
+DEAL:cg::Convergence step 39 value 1.10898
+DEAL::Norm of the solution: 196549.
+DEAL:cg::Starting value 125.000
+DEAL:cg::Convergence step 39 value 1.10898
+DEAL::Norm of the solution: 196549.
+DEAL::Number of calls to special vmult: 39
+DEAL:cg::Starting value 125.000
+DEAL:cg::Convergence step 39 value 1.10898
+DEAL::Norm of the solution: 196549.