]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Test cases 8728/head
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Thu, 19 Sep 2019 16:02:26 +0000 (18:02 +0200)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Thu, 19 Sep 2019 16:33:21 +0000 (18:33 +0200)
tests/matrix_free/pre_and_post_loops_01.cc [new file with mode: 0644]
tests/matrix_free/pre_and_post_loops_01.with_mpi=true.with_p4est=true.mpirun=1.output [new file with mode: 0644]
tests/matrix_free/pre_and_post_loops_01.with_mpi=true.with_p4est=true.mpirun=3.output [new file with mode: 0644]
tests/matrix_free/pre_and_post_loops_02.cc [new file with mode: 0644]
tests/matrix_free/pre_and_post_loops_02.with_mpi=true.with_p4est=true.mpirun=3.output [new file with mode: 0644]
tests/matrix_free/pre_and_post_loops_03.cc [new file with mode: 0644]
tests/matrix_free/pre_and_post_loops_03.with_mpi=true.with_p4est=true.mpirun=3.output [new file with mode: 0644]
tests/matrix_free/pre_and_post_loops_04.cc [new file with mode: 0644]
tests/matrix_free/pre_and_post_loops_04.with_mpi=true.with_p4est=true.mpirun=2.output [new file with mode: 0644]
tests/matrix_free/pre_and_post_loops_05.cc [new file with mode: 0644]
tests/matrix_free/pre_and_post_loops_05.with_mpi=true.with_p4est=true.mpirun=2.output [new file with mode: 0644]

diff --git a/tests/matrix_free/pre_and_post_loops_01.cc b/tests/matrix_free/pre_and_post_loops_01.cc
new file mode 100644 (file)
index 0000000..7551a41
--- /dev/null
@@ -0,0 +1,205 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// test the correctness of the matrix-free cell loop with additional
+// operation_before_loop and operation_after_loop lambdas
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/sparse_matrix.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+#include "matrix_vector_mf.h"
+
+
+
+template <int dim,
+          int fe_degree,
+          int n_q_points_1d = fe_degree + 1,
+          typename Number   = double>
+class Matrix
+{
+public:
+  Matrix(const MatrixFree<dim, Number> &data_in)
+    : data(data_in)
+  {}
+
+  void
+  vmult(LinearAlgebra::distributed::Vector<Number> &      dst,
+        const LinearAlgebra::distributed::Vector<Number> &src) const
+  {
+    const std::function<void(const MatrixFree<dim, Number> &,
+                             LinearAlgebra::distributed::Vector<Number> &,
+                             const LinearAlgebra::distributed::Vector<Number> &,
+                             const std::pair<unsigned int, unsigned int> &)>
+      wrap = helmholtz_operator<dim,
+                                fe_degree,
+                                LinearAlgebra::distributed::Vector<Number>,
+                                n_q_points_1d>;
+    dst    = 0;
+    data.cell_loop(wrap, dst, src);
+    for (auto i : data.get_constrained_dofs())
+      dst.local_element(i) = src.local_element(i);
+  }
+
+  void
+  vmult_with_update_basic(
+    LinearAlgebra::distributed::Vector<Number> &dst,
+    LinearAlgebra::distributed::Vector<Number> &src,
+    LinearAlgebra::distributed::Vector<Number> &other_vector) const
+  {
+    src.add(1.5, other_vector);
+    other_vector.add(0.7, dst);
+    dst = 0;
+    vmult(dst, src);
+    dst.sadd(0.63, -1.3, other_vector);
+  }
+
+  void
+  vmult_with_update_merged(
+    LinearAlgebra::distributed::Vector<Number> &dst,
+    LinearAlgebra::distributed::Vector<Number> &src,
+    LinearAlgebra::distributed::Vector<Number> &other_vector) const
+  {
+    const std::function<void(const MatrixFree<dim, Number> &,
+                             LinearAlgebra::distributed::Vector<Number> &,
+                             const LinearAlgebra::distributed::Vector<Number> &,
+                             const std::pair<unsigned int, unsigned int> &)>
+      wrap = helmholtz_operator<dim,
+                                fe_degree,
+                                LinearAlgebra::distributed::Vector<Number>,
+                                n_q_points_1d>;
+    data.cell_loop(
+      wrap,
+      dst,
+      src,
+      // operation before cell operation
+      [&](const unsigned int start_range, const unsigned int end_range) {
+        for (unsigned int i = start_range; i < end_range; ++i)
+          {
+            src.local_element(i) += 1.5 * other_vector.local_element(i);
+            other_vector.local_element(i) += 0.7 * dst.local_element(i);
+            dst.local_element(i) = 0;
+          }
+      },
+      // operation after cell operation
+      [&](const unsigned int start_range, const unsigned int end_range) {
+        for (unsigned int i = start_range; i < end_range; ++i)
+          {
+            dst.local_element(i) =
+              0.63 * dst.local_element(i) - 1.3 * other_vector.local_element(i);
+          }
+      });
+  }
+
+private:
+  const MatrixFree<dim, Number> &data;
+};
+
+
+
+template <int dim, int fe_degree, typename number>
+void
+test()
+{
+  parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
+  GridGenerator::hyper_cube(tria);
+  tria.refine_global(7 - dim);
+
+  FE_Q<dim>       fe(fe_degree);
+  DoFHandler<dim> dof(tria);
+  dof.distribute_dofs(fe);
+  AffineConstraints<double> constraints;
+  constraints.close();
+
+  deallog << "Testing " << dof.get_fe().get_name() << std::endl;
+
+  MatrixFree<dim, number> mf_data;
+  {
+    const QGauss<1>                                  quad(fe_degree + 1);
+    typename MatrixFree<dim, number>::AdditionalData data;
+    data.tasks_parallel_scheme = MatrixFree<dim, number>::AdditionalData::none;
+    mf_data.reinit(dof, constraints, quad, data);
+  }
+
+  Matrix<dim, fe_degree, fe_degree + 1, number> mf(mf_data);
+  LinearAlgebra::distributed::Vector<number>    vec1, vec2, vec3;
+  mf_data.initialize_dof_vector(vec1);
+  mf_data.initialize_dof_vector(vec2);
+  mf_data.initialize_dof_vector(vec3);
+
+  for (unsigned int i = 0; i < vec1.local_size(); ++i)
+    {
+      double entry          = random_value<double>();
+      vec1.local_element(i) = entry;
+      entry                 = random_value<double>();
+      vec2.local_element(i) = entry;
+      entry                 = random_value<double>();
+      vec3.local_element(i) = entry;
+    }
+
+  LinearAlgebra::distributed::Vector<number> ref1 = vec1;
+  LinearAlgebra::distributed::Vector<number> ref2 = vec2;
+  LinearAlgebra::distributed::Vector<number> ref3 = vec3;
+
+  mf.vmult_with_update_basic(ref3, ref2, ref1);
+  mf.vmult_with_update_basic(vec3, vec2, vec1);
+
+  vec3 -= ref3;
+  deallog << "Error in 1x merged operation: " << vec3.linfty_norm()
+          << std::endl;
+
+  ref3 = 0;
+
+  mf.vmult_with_update_basic(ref1, ref2, ref3);
+  mf.vmult_with_update_basic(vec1, vec2, vec3);
+
+  mf.vmult_with_update_basic(ref1, ref2, ref3);
+  mf.vmult_with_update_basic(vec1, vec2, vec3);
+
+  vec3 -= ref3;
+  deallog << "Error in 2x merged operation: " << vec3.linfty_norm()
+          << std::endl;
+}
+
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1);
+
+  mpi_initlog();
+  deallog << std::setprecision(3);
+
+  test<2, 3, double>();
+  test<2, 3, float>();
+  test<3, 3, double>();
+}
diff --git a/tests/matrix_free/pre_and_post_loops_01.with_mpi=true.with_p4est=true.mpirun=1.output b/tests/matrix_free/pre_and_post_loops_01.with_mpi=true.with_p4est=true.mpirun=1.output
new file mode 100644 (file)
index 0000000..57a0159
--- /dev/null
@@ -0,0 +1,10 @@
+
+DEAL::Testing FE_Q<2>(3)
+DEAL::Error in 1x merged operation: 0.00
+DEAL::Error in 2x merged operation: 0.00
+DEAL::Testing FE_Q<2>(3)
+DEAL::Error in 1x merged operation: 0.00
+DEAL::Error in 2x merged operation: 0.00
+DEAL::Testing FE_Q<3>(3)
+DEAL::Error in 1x merged operation: 0.00
+DEAL::Error in 2x merged operation: 0.00
diff --git a/tests/matrix_free/pre_and_post_loops_01.with_mpi=true.with_p4est=true.mpirun=3.output b/tests/matrix_free/pre_and_post_loops_01.with_mpi=true.with_p4est=true.mpirun=3.output
new file mode 100644 (file)
index 0000000..57a0159
--- /dev/null
@@ -0,0 +1,10 @@
+
+DEAL::Testing FE_Q<2>(3)
+DEAL::Error in 1x merged operation: 0.00
+DEAL::Error in 2x merged operation: 0.00
+DEAL::Testing FE_Q<2>(3)
+DEAL::Error in 1x merged operation: 0.00
+DEAL::Error in 2x merged operation: 0.00
+DEAL::Testing FE_Q<3>(3)
+DEAL::Error in 1x merged operation: 0.00
+DEAL::Error in 2x merged operation: 0.00
diff --git a/tests/matrix_free/pre_and_post_loops_02.cc b/tests/matrix_free/pre_and_post_loops_02.cc
new file mode 100644 (file)
index 0000000..e999c60
--- /dev/null
@@ -0,0 +1,218 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// test the correctness of the matrix-free cell loop with additional
+// operation_before_loop and operation_after_loop lambdas. similar to
+// pre_and_post_loops_01 but using Dirichlet boundary conditions as
+// constraints
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/sparse_matrix.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+#include "matrix_vector_mf.h"
+
+
+
+template <int dim,
+          int fe_degree,
+          int n_q_points_1d = fe_degree + 1,
+          typename Number   = double>
+class Matrix
+{
+public:
+  Matrix(const MatrixFree<dim, Number> &data_in)
+    : data(data_in)
+  {}
+
+  void
+  vmult(LinearAlgebra::distributed::Vector<Number> &      dst,
+        const LinearAlgebra::distributed::Vector<Number> &src) const
+  {
+    const std::function<void(const MatrixFree<dim, Number> &,
+                             LinearAlgebra::distributed::Vector<Number> &,
+                             const LinearAlgebra::distributed::Vector<Number> &,
+                             const std::pair<unsigned int, unsigned int> &)>
+      wrap = helmholtz_operator<dim,
+                                fe_degree,
+                                LinearAlgebra::distributed::Vector<Number>,
+                                n_q_points_1d>;
+    dst    = 0;
+    data.cell_loop(wrap, dst, src);
+    for (auto i : data.get_constrained_dofs())
+      dst.local_element(i) = src.local_element(i);
+  }
+
+  void
+  vmult_with_update_basic(
+    LinearAlgebra::distributed::Vector<Number> &dst,
+    LinearAlgebra::distributed::Vector<Number> &src,
+    LinearAlgebra::distributed::Vector<Number> &other_vector) const
+  {
+    src.add(1.5, other_vector);
+    other_vector.add(0.7, dst);
+    dst = 0;
+    vmult(dst, src);
+    dst.sadd(0.63, -1.3, other_vector);
+  }
+
+  void
+  vmult_with_update_merged(
+    LinearAlgebra::distributed::Vector<Number> &dst,
+    LinearAlgebra::distributed::Vector<Number> &src,
+    LinearAlgebra::distributed::Vector<Number> &other_vector) const
+  {
+    const std::function<void(const MatrixFree<dim, Number> &,
+                             LinearAlgebra::distributed::Vector<Number> &,
+                             const LinearAlgebra::distributed::Vector<Number> &,
+                             const std::pair<unsigned int, unsigned int> &)>
+      wrap = helmholtz_operator<dim,
+                                fe_degree,
+                                LinearAlgebra::distributed::Vector<Number>,
+                                n_q_points_1d>;
+    data.cell_loop(
+      wrap,
+      dst,
+      src,
+      // operation before cell operation
+      [&](const unsigned int start_range, const unsigned int end_range) {
+        for (unsigned int i = start_range; i < end_range; ++i)
+          {
+            src.local_element(i) += 1.5 * other_vector.local_element(i);
+            other_vector.local_element(i) += 0.7 * dst.local_element(i);
+            dst.local_element(i) = 0;
+          }
+      },
+      // operation after cell operation
+      [&](const unsigned int start_range, const unsigned int end_range) {
+        for (unsigned int i = start_range; i < end_range; ++i)
+          {
+            dst.local_element(i) =
+              0.63 * dst.local_element(i) - 1.3 * other_vector.local_element(i);
+          }
+      });
+  }
+
+private:
+  const MatrixFree<dim, Number> &data;
+};
+
+
+
+template <int dim, int fe_degree, typename number>
+void
+test()
+{
+  parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
+  GridGenerator::hyper_cube(tria);
+  tria.refine_global(7 - dim);
+
+  FE_Q<dim>       fe(fe_degree);
+  DoFHandler<dim> dof(tria);
+  dof.distribute_dofs(fe);
+  AffineConstraints<double> constraints;
+  IndexSet                  relevant_dofs;
+  DoFTools::extract_locally_relevant_dofs(dof, relevant_dofs);
+  constraints.reinit(relevant_dofs);
+  VectorTools::interpolate_boundary_values(dof,
+                                           0,
+                                           ZeroFunction<dim>(),
+                                           constraints);
+  constraints.close();
+
+  deallog << "Testing " << dof.get_fe().get_name() << std::endl;
+
+  MatrixFree<dim, number> mf_data;
+  {
+    const QGauss<1>                                  quad(fe_degree + 1);
+    typename MatrixFree<dim, number>::AdditionalData data;
+    data.tasks_parallel_scheme = MatrixFree<dim, number>::AdditionalData::none;
+    mf_data.reinit(dof, constraints, quad, data);
+  }
+
+  Matrix<dim, fe_degree, fe_degree + 1, number> mf(mf_data);
+  LinearAlgebra::distributed::Vector<number>    vec1, vec2, vec3;
+  mf_data.initialize_dof_vector(vec1);
+  mf_data.initialize_dof_vector(vec2);
+  mf_data.initialize_dof_vector(vec3);
+
+  for (unsigned int i = 0; i < vec1.local_size(); ++i)
+    {
+      if (constraints.is_constrained(
+            vec1.get_partitioner()->local_to_global(i)))
+        continue;
+
+      double entry          = random_value<double>();
+      vec1.local_element(i) = entry;
+      entry                 = random_value<double>();
+      vec2.local_element(i) = entry;
+      entry                 = random_value<double>();
+      vec3.local_element(i) = entry;
+    }
+
+  LinearAlgebra::distributed::Vector<number> ref1 = vec1;
+  LinearAlgebra::distributed::Vector<number> ref2 = vec2;
+  LinearAlgebra::distributed::Vector<number> ref3 = vec3;
+
+  mf.vmult_with_update_basic(ref3, ref2, ref1);
+  mf.vmult_with_update_basic(vec3, vec2, vec1);
+
+  vec3 -= ref3;
+  deallog << "Error in 1x merged operation: " << vec3.linfty_norm()
+          << std::endl;
+
+  ref3 = 0;
+
+  mf.vmult_with_update_basic(ref1, ref2, ref3);
+  mf.vmult_with_update_basic(vec1, vec2, vec3);
+
+  mf.vmult_with_update_basic(ref1, ref2, ref3);
+  mf.vmult_with_update_basic(vec1, vec2, vec3);
+
+  vec3 -= ref3;
+  deallog << "Error in 2x merged operation: " << vec3.linfty_norm()
+          << std::endl;
+}
+
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1);
+
+  mpi_initlog();
+  deallog << std::setprecision(3);
+
+  test<2, 3, double>();
+  test<2, 3, float>();
+  test<3, 3, double>();
+}
diff --git a/tests/matrix_free/pre_and_post_loops_02.with_mpi=true.with_p4est=true.mpirun=3.output b/tests/matrix_free/pre_and_post_loops_02.with_mpi=true.with_p4est=true.mpirun=3.output
new file mode 100644 (file)
index 0000000..57a0159
--- /dev/null
@@ -0,0 +1,10 @@
+
+DEAL::Testing FE_Q<2>(3)
+DEAL::Error in 1x merged operation: 0.00
+DEAL::Error in 2x merged operation: 0.00
+DEAL::Testing FE_Q<2>(3)
+DEAL::Error in 1x merged operation: 0.00
+DEAL::Error in 2x merged operation: 0.00
+DEAL::Testing FE_Q<3>(3)
+DEAL::Error in 1x merged operation: 0.00
+DEAL::Error in 2x merged operation: 0.00
diff --git a/tests/matrix_free/pre_and_post_loops_03.cc b/tests/matrix_free/pre_and_post_loops_03.cc
new file mode 100644 (file)
index 0000000..5973585
--- /dev/null
@@ -0,0 +1,207 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// test the correctness of the matrix-free cell loop with additional
+// operation_before_loop and operation_after_loop lambdas. Similar to
+// pre_and_post_loops_01 but working on a hyper cube where we have some
+// different orientations, stressing different indices
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/sparse_matrix.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+#include "matrix_vector_mf.h"
+
+
+
+template <int dim,
+          int fe_degree,
+          int n_q_points_1d = fe_degree + 1,
+          typename Number   = double>
+class Matrix
+{
+public:
+  Matrix(const MatrixFree<dim, Number> &data_in)
+    : data(data_in)
+  {}
+
+  void
+  vmult(LinearAlgebra::distributed::Vector<Number> &      dst,
+        const LinearAlgebra::distributed::Vector<Number> &src) const
+  {
+    const std::function<void(const MatrixFree<dim, Number> &,
+                             LinearAlgebra::distributed::Vector<Number> &,
+                             const LinearAlgebra::distributed::Vector<Number> &,
+                             const std::pair<unsigned int, unsigned int> &)>
+      wrap = helmholtz_operator<dim,
+                                fe_degree,
+                                LinearAlgebra::distributed::Vector<Number>,
+                                n_q_points_1d>;
+    dst    = 0;
+    data.cell_loop(wrap, dst, src);
+    for (auto i : data.get_constrained_dofs())
+      dst.local_element(i) = src.local_element(i);
+  }
+
+  void
+  vmult_with_update_basic(
+    LinearAlgebra::distributed::Vector<Number> &dst,
+    LinearAlgebra::distributed::Vector<Number> &src,
+    LinearAlgebra::distributed::Vector<Number> &other_vector) const
+  {
+    src.add(1.5, other_vector);
+    other_vector.add(0.7, dst);
+    dst = 0;
+    vmult(dst, src);
+    dst.sadd(0.63, -1.3, other_vector);
+  }
+
+  void
+  vmult_with_update_merged(
+    LinearAlgebra::distributed::Vector<Number> &dst,
+    LinearAlgebra::distributed::Vector<Number> &src,
+    LinearAlgebra::distributed::Vector<Number> &other_vector) const
+  {
+    const std::function<void(const MatrixFree<dim, Number> &,
+                             LinearAlgebra::distributed::Vector<Number> &,
+                             const LinearAlgebra::distributed::Vector<Number> &,
+                             const std::pair<unsigned int, unsigned int> &)>
+      wrap = helmholtz_operator<dim,
+                                fe_degree,
+                                LinearAlgebra::distributed::Vector<Number>,
+                                n_q_points_1d>;
+    data.cell_loop(
+      wrap,
+      dst,
+      src,
+      // operation before cell operation
+      [&](const unsigned int start_range, const unsigned int end_range) {
+        for (unsigned int i = start_range; i < end_range; ++i)
+          {
+            src.local_element(i) += 1.5 * other_vector.local_element(i);
+            other_vector.local_element(i) += 0.7 * dst.local_element(i);
+            dst.local_element(i) = 0;
+          }
+      },
+      // operation after cell operation
+      [&](const unsigned int start_range, const unsigned int end_range) {
+        for (unsigned int i = start_range; i < end_range; ++i)
+          {
+            dst.local_element(i) =
+              0.63 * dst.local_element(i) - 1.3 * other_vector.local_element(i);
+          }
+      });
+  }
+
+private:
+  const MatrixFree<dim, Number> &data;
+};
+
+
+
+template <int dim, int fe_degree, typename number>
+void
+test()
+{
+  parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
+  GridGenerator::hyper_ball(tria);
+  tria.refine_global(6 - dim);
+
+  FE_Q<dim>       fe(fe_degree);
+  DoFHandler<dim> dof(tria);
+  dof.distribute_dofs(fe);
+  AffineConstraints<double> constraints;
+  constraints.close();
+
+  deallog << "Testing " << dof.get_fe().get_name() << std::endl;
+
+  MatrixFree<dim, number> mf_data;
+  {
+    const QGauss<1>                                  quad(fe_degree + 1);
+    typename MatrixFree<dim, number>::AdditionalData data;
+    data.tasks_parallel_scheme = MatrixFree<dim, number>::AdditionalData::none;
+    mf_data.reinit(dof, constraints, quad, data);
+  }
+
+  Matrix<dim, fe_degree, fe_degree + 1, number> mf(mf_data);
+  LinearAlgebra::distributed::Vector<number>    vec1, vec2, vec3;
+  mf_data.initialize_dof_vector(vec1);
+  mf_data.initialize_dof_vector(vec2);
+  mf_data.initialize_dof_vector(vec3);
+
+  for (unsigned int i = 0; i < vec1.local_size(); ++i)
+    {
+      double entry          = random_value<double>();
+      vec1.local_element(i) = entry;
+      entry                 = random_value<double>();
+      vec2.local_element(i) = entry;
+      entry                 = random_value<double>();
+      vec3.local_element(i) = entry;
+    }
+
+  LinearAlgebra::distributed::Vector<number> ref1 = vec1;
+  LinearAlgebra::distributed::Vector<number> ref2 = vec2;
+  LinearAlgebra::distributed::Vector<number> ref3 = vec3;
+
+  mf.vmult_with_update_basic(ref3, ref2, ref1);
+  mf.vmult_with_update_basic(vec3, vec2, vec1);
+
+  vec3 -= ref3;
+  deallog << "Error in 1x merged operation: " << vec3.linfty_norm()
+          << std::endl;
+
+  ref3 = 0;
+
+  mf.vmult_with_update_basic(ref1, ref2, ref3);
+  mf.vmult_with_update_basic(vec1, vec2, vec3);
+
+  mf.vmult_with_update_basic(ref1, ref2, ref3);
+  mf.vmult_with_update_basic(vec1, vec2, vec3);
+
+  vec3 -= ref3;
+  deallog << "Error in 2x merged operation: " << vec3.linfty_norm()
+          << std::endl;
+}
+
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1);
+
+  mpi_initlog();
+  deallog << std::setprecision(3);
+
+  test<2, 3, double>();
+  test<2, 3, float>();
+  test<3, 3, double>();
+}
diff --git a/tests/matrix_free/pre_and_post_loops_03.with_mpi=true.with_p4est=true.mpirun=3.output b/tests/matrix_free/pre_and_post_loops_03.with_mpi=true.with_p4est=true.mpirun=3.output
new file mode 100644 (file)
index 0000000..57a0159
--- /dev/null
@@ -0,0 +1,10 @@
+
+DEAL::Testing FE_Q<2>(3)
+DEAL::Error in 1x merged operation: 0.00
+DEAL::Error in 2x merged operation: 0.00
+DEAL::Testing FE_Q<2>(3)
+DEAL::Error in 1x merged operation: 0.00
+DEAL::Error in 2x merged operation: 0.00
+DEAL::Testing FE_Q<3>(3)
+DEAL::Error in 1x merged operation: 0.00
+DEAL::Error in 2x merged operation: 0.00
diff --git a/tests/matrix_free/pre_and_post_loops_04.cc b/tests/matrix_free/pre_and_post_loops_04.cc
new file mode 100644 (file)
index 0000000..c0d71b7
--- /dev/null
@@ -0,0 +1,202 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// test the correctness of the matrix-free cell loop with additional
+// operation_before_loop and operation_after_loop lambdas. Similar to
+// pre_and_post_loops_01 but using a pointer to a member function in
+// MatrixFree::cell_loop rather than a separate function
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/sparse_matrix.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+#include "matrix_vector_mf.h"
+
+
+
+template <int dim,
+          int fe_degree,
+          int n_q_points_1d = fe_degree + 1,
+          typename Number   = double>
+class Matrix
+{
+public:
+  Matrix(const MatrixFree<dim, Number> &data_in)
+    : data(data_in)
+  {}
+
+  void
+  vmult(LinearAlgebra::distributed::Vector<Number> &      dst,
+        const LinearAlgebra::distributed::Vector<Number> &src) const
+  {
+    dst = 0;
+    data.cell_loop(&Matrix::local_apply, this, dst, src);
+    for (auto i : data.get_constrained_dofs())
+      dst.local_element(i) = src.local_element(i);
+  }
+
+  void
+  vmult_with_update_basic(
+    LinearAlgebra::distributed::Vector<Number> &dst,
+    LinearAlgebra::distributed::Vector<Number> &src,
+    LinearAlgebra::distributed::Vector<Number> &other_vector) const
+  {
+    src.add(1.5, other_vector);
+    other_vector.add(0.7, dst);
+    dst = 0;
+    vmult(dst, src);
+    dst.sadd(0.63, -1.3, other_vector);
+  }
+
+  void
+  vmult_with_update_merged(
+    LinearAlgebra::distributed::Vector<Number> &dst,
+    LinearAlgebra::distributed::Vector<Number> &src,
+    LinearAlgebra::distributed::Vector<Number> &other_vector) const
+  {
+    data.cell_loop(
+      &Matrix::local_apply,
+      this,
+      dst,
+      src,
+      // operation before cell operation
+      [&](const unsigned int start_range, const unsigned int end_range) {
+        for (unsigned int i = start_range; i < end_range; ++i)
+          {
+            src.local_element(i) += 1.5 * other_vector.local_element(i);
+            other_vector.local_element(i) += 0.7 * dst.local_element(i);
+            dst.local_element(i) = 0;
+          }
+      },
+      // operation after cell operation
+      [&](const unsigned int start_range, const unsigned int end_range) {
+        for (unsigned int i = start_range; i < end_range; ++i)
+          {
+            dst.local_element(i) =
+              0.63 * dst.local_element(i) - 1.3 * other_vector.local_element(i);
+          }
+      });
+  }
+
+private:
+  const MatrixFree<dim, Number> &data;
+
+  void
+  local_apply(const MatrixFree<dim, Number> &,
+              LinearAlgebra::distributed::Vector<Number> &      dst,
+              const LinearAlgebra::distributed::Vector<Number> &src,
+              const std::pair<unsigned int, unsigned int> &     range) const
+  {
+    helmholtz_operator<dim,
+                       fe_degree,
+                       LinearAlgebra::distributed::Vector<Number>,
+                       n_q_points_1d>(data, dst, src, range);
+  }
+};
+
+
+
+template <int dim, int fe_degree, typename number>
+void
+test()
+{
+  parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
+  GridGenerator::hyper_cube(tria);
+  tria.refine_global(7 - dim);
+
+  FE_Q<dim>       fe(fe_degree);
+  DoFHandler<dim> dof(tria);
+  dof.distribute_dofs(fe);
+  AffineConstraints<double> constraints;
+  constraints.close();
+
+  deallog << "Testing " << dof.get_fe().get_name() << std::endl;
+
+  MatrixFree<dim, number> mf_data;
+  {
+    const QGauss<1>                                  quad(fe_degree + 1);
+    typename MatrixFree<dim, number>::AdditionalData data;
+    data.tasks_parallel_scheme = MatrixFree<dim, number>::AdditionalData::none;
+    mf_data.reinit(dof, constraints, quad, data);
+  }
+
+  Matrix<dim, fe_degree, fe_degree + 1, number> mf(mf_data);
+  LinearAlgebra::distributed::Vector<number>    vec1, vec2, vec3;
+  mf_data.initialize_dof_vector(vec1);
+  mf_data.initialize_dof_vector(vec2);
+  mf_data.initialize_dof_vector(vec3);
+
+  for (unsigned int i = 0; i < vec1.local_size(); ++i)
+    {
+      double entry          = random_value<double>();
+      vec1.local_element(i) = entry;
+      entry                 = random_value<double>();
+      vec2.local_element(i) = entry;
+      entry                 = random_value<double>();
+      vec3.local_element(i) = entry;
+    }
+
+  LinearAlgebra::distributed::Vector<number> ref1 = vec1;
+  LinearAlgebra::distributed::Vector<number> ref2 = vec2;
+  LinearAlgebra::distributed::Vector<number> ref3 = vec3;
+
+  mf.vmult_with_update_basic(ref3, ref2, ref1);
+  mf.vmult_with_update_basic(vec3, vec2, vec1);
+
+  vec3 -= ref3;
+  deallog << "Error in 1x merged operation: " << vec3.linfty_norm()
+          << std::endl;
+
+  ref3 = 0;
+
+  mf.vmult_with_update_basic(ref1, ref2, ref3);
+  mf.vmult_with_update_basic(vec1, vec2, vec3);
+
+  mf.vmult_with_update_basic(ref1, ref2, ref3);
+  mf.vmult_with_update_basic(vec1, vec2, vec3);
+
+  vec3 -= ref3;
+  deallog << "Error in 2x merged operation: " << vec3.linfty_norm()
+          << std::endl;
+}
+
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1);
+
+  mpi_initlog();
+  deallog << std::setprecision(3);
+
+  test<2, 3, double>();
+}
diff --git a/tests/matrix_free/pre_and_post_loops_04.with_mpi=true.with_p4est=true.mpirun=2.output b/tests/matrix_free/pre_and_post_loops_04.with_mpi=true.with_p4est=true.mpirun=2.output
new file mode 100644 (file)
index 0000000..bb542bb
--- /dev/null
@@ -0,0 +1,4 @@
+
+DEAL::Testing FE_Q<2>(3)
+DEAL::Error in 1x merged operation: 0.00
+DEAL::Error in 2x merged operation: 0.00
diff --git a/tests/matrix_free/pre_and_post_loops_05.cc b/tests/matrix_free/pre_and_post_loops_05.cc
new file mode 100644 (file)
index 0000000..8196125
--- /dev/null
@@ -0,0 +1,221 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// test the correctness of the matrix-free cell loop with additional
+// operation_before_loop and operation_after_loop lambdas. Similar to
+// pre_and_post_loops_01 but using a system of two DoFHandler objects,
+// operating on the second one
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/sparse_matrix.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+#include "matrix_vector_mf.h"
+
+
+
+template <int dim,
+          int fe_degree,
+          int n_q_points_1d = fe_degree + 1,
+          typename Number   = double>
+class Matrix
+{
+public:
+  Matrix(const MatrixFree<dim, Number> &data_in)
+    : data(data_in)
+  {}
+
+  void
+  vmult(LinearAlgebra::distributed::Vector<Number> &      dst,
+        const LinearAlgebra::distributed::Vector<Number> &src) const
+  {
+    dst = 0;
+    data.cell_loop(&Matrix::local_apply, this, dst, src);
+    for (auto i : data.get_constrained_dofs())
+      dst.local_element(i) = src.local_element(i);
+  }
+
+  void
+  vmult_with_update_basic(
+    LinearAlgebra::distributed::Vector<Number> &dst,
+    LinearAlgebra::distributed::Vector<Number> &src,
+    LinearAlgebra::distributed::Vector<Number> &other_vector) const
+  {
+    src.add(1.5, other_vector);
+    other_vector.add(0.7, dst);
+    dst = 0;
+    vmult(dst, src);
+    dst.sadd(0.63, -1.3, other_vector);
+  }
+
+  void
+  vmult_with_update_merged(
+    LinearAlgebra::distributed::Vector<Number> &dst,
+    LinearAlgebra::distributed::Vector<Number> &src,
+    LinearAlgebra::distributed::Vector<Number> &other_vector) const
+  {
+    data.cell_loop(
+      &Matrix::local_apply,
+      this,
+      dst,
+      src,
+      // operation before cell operation
+      [&](const unsigned int start_range, const unsigned int end_range) {
+        for (unsigned int i = start_range; i < end_range; ++i)
+          {
+            src.local_element(i) += 1.5 * other_vector.local_element(i);
+            other_vector.local_element(i) += 0.7 * dst.local_element(i);
+            dst.local_element(i) = 0;
+          }
+      },
+      // operation after cell operation
+      [&](const unsigned int start_range, const unsigned int end_range) {
+        for (unsigned int i = start_range; i < end_range; ++i)
+          {
+            dst.local_element(i) =
+              0.63 * dst.local_element(i) - 1.3 * other_vector.local_element(i);
+          }
+      },
+      // dof handler component for the two lambdas
+      1);
+  }
+
+private:
+  const MatrixFree<dim, Number> &data;
+
+  void
+  local_apply(const MatrixFree<dim, Number> &,
+              LinearAlgebra::distributed::Vector<Number> &      dst,
+              const LinearAlgebra::distributed::Vector<Number> &src,
+              const std::pair<unsigned int, unsigned int> &     range) const
+  {
+    FEEvaluation<dim, fe_degree, n_q_points_1d, 1, Number> fe_eval(data, 1);
+
+    for (unsigned int cell = range.first; cell < range.second; ++cell)
+      {
+        fe_eval.reinit(cell);
+        fe_eval.read_dof_values(src);
+        fe_eval.evaluate(false, true);
+        for (unsigned int q = 0; q < fe_eval.n_q_points; ++q)
+          fe_eval.submit_gradient(fe_eval.get_gradient(q), q);
+        fe_eval.integrate(false, true);
+        fe_eval.distribute_local_to_global(dst);
+      }
+  }
+};
+
+
+
+template <int dim, int fe_degree, typename number>
+void
+test()
+{
+  parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
+  GridGenerator::hyper_cube(tria);
+  tria.refine_global(7 - dim);
+
+  FE_Q<dim>       fe(fe_degree);
+  DoFHandler<dim> dof(tria);
+  dof.distribute_dofs(fe);
+  FE_Q<dim>       fe2(1);
+  DoFHandler<dim> dof2(tria);
+  dof2.distribute_dofs(fe2);
+  AffineConstraints<double> constraints;
+  constraints.close();
+
+  deallog << "Testing " << dof.get_fe().get_name() << std::endl;
+
+  MatrixFree<dim, number> mf_data;
+  {
+    const QGauss<1>                                  quad(fe_degree + 1);
+    typename MatrixFree<dim, number>::AdditionalData data;
+    data.tasks_parallel_scheme = MatrixFree<dim, number>::AdditionalData::none;
+    mf_data.reinit(std::vector<const DoFHandler<dim> *>({&dof2, &dof}),
+                   std::vector<const AffineConstraints<double> *>(
+                     {&constraints, &constraints}),
+                   std::vector<Quadrature<1>>({quad}),
+                   data);
+  }
+
+  Matrix<dim, fe_degree, fe_degree + 1, number> mf(mf_data);
+  LinearAlgebra::distributed::Vector<number>    vec1, vec2, vec3;
+  mf_data.initialize_dof_vector(vec1, 1);
+  mf_data.initialize_dof_vector(vec2, 1);
+  mf_data.initialize_dof_vector(vec3, 1);
+
+  for (unsigned int i = 0; i < vec1.local_size(); ++i)
+    {
+      double entry          = random_value<double>();
+      vec1.local_element(i) = entry;
+      entry                 = random_value<double>();
+      vec2.local_element(i) = entry;
+      entry                 = random_value<double>();
+      vec3.local_element(i) = entry;
+    }
+
+  LinearAlgebra::distributed::Vector<number> ref1 = vec1;
+  LinearAlgebra::distributed::Vector<number> ref2 = vec2;
+  LinearAlgebra::distributed::Vector<number> ref3 = vec3;
+
+  mf.vmult_with_update_basic(ref3, ref2, ref1);
+  mf.vmult_with_update_basic(vec3, vec2, vec1);
+
+  vec3 -= ref3;
+  deallog << "Error in 1x merged operation: " << vec3.linfty_norm()
+          << std::endl;
+
+  ref3 = 0;
+
+  mf.vmult_with_update_basic(ref1, ref2, ref3);
+  mf.vmult_with_update_basic(vec1, vec2, vec3);
+
+  mf.vmult_with_update_basic(ref1, ref2, ref3);
+  mf.vmult_with_update_basic(vec1, vec2, vec3);
+
+  vec3 -= ref3;
+  deallog << "Error in 2x merged operation: " << vec3.linfty_norm()
+          << std::endl;
+}
+
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1);
+
+  mpi_initlog();
+  deallog << std::setprecision(3);
+
+  test<2, 3, double>();
+  test<2, 3, float>();
+  test<3, 2, double>();
+}
diff --git a/tests/matrix_free/pre_and_post_loops_05.with_mpi=true.with_p4est=true.mpirun=2.output b/tests/matrix_free/pre_and_post_loops_05.with_mpi=true.with_p4est=true.mpirun=2.output
new file mode 100644 (file)
index 0000000..02b75aa
--- /dev/null
@@ -0,0 +1,10 @@
+
+DEAL::Testing FE_Q<2>(3)
+DEAL::Error in 1x merged operation: 0.00
+DEAL::Error in 2x merged operation: 0.00
+DEAL::Testing FE_Q<2>(3)
+DEAL::Error in 1x merged operation: 0.00
+DEAL::Error in 2x merged operation: 0.00
+DEAL::Testing FE_Q<3>(2)
+DEAL::Error in 1x merged operation: 0.00
+DEAL::Error in 2x merged operation: 0.00

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.