]> https://gitweb.dealii.org/ - dealii.git/commitdiff
MatrixFree pre and post operation for loop()
authorSebastian Proell <sebastian.proell@tum.de>
Sat, 8 Oct 2022 11:25:12 +0000 (13:25 +0200)
committerSebastian Proell <sebastian.proell@tum.de>
Tue, 11 Oct 2022 15:44:33 +0000 (17:44 +0200)
doc/news/changes/minor/20221011Proell [new file with mode: 0644]
include/deal.II/matrix_free/matrix_free.h
tests/matrix_free/pre_and_post_loops_07.cc [new file with mode: 0644]
tests/matrix_free/pre_and_post_loops_07.with_p4est=true.mpirun=1.output [new file with mode: 0644]
tests/matrix_free/pre_and_post_loops_07.with_p4est=true.mpirun=3.output [new file with mode: 0644]
tests/matrix_free/pre_and_post_loops_08.cc [new file with mode: 0644]
tests/matrix_free/pre_and_post_loops_08.with_p4est=true.mpirun=2.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20221011Proell b/doc/news/changes/minor/20221011Proell
new file mode 100644 (file)
index 0000000..11ea1f8
--- /dev/null
@@ -0,0 +1,5 @@
+New: A new overload for MatrixFree::loop allows
+to specify a pre- and post-operation, similar
+to MatrixFree::cell_loop.
+<br>
+(Sebastian Proell, 2022/10/11)
index 39fa06c6b46c16a9efaf95f9bb5290aa88388a49..aba89df07581b01c7572fdd5ee829fbed16deed3 100644 (file)
@@ -1245,6 +1245,223 @@ public:
        const DataAccessOnFaces src_vector_face_access =
          DataAccessOnFaces::unspecified) const;
 
+  /**
+   * This function is similar to the loop method above, but adds two additional
+   * functors to execute some additional work before and after the cell, face
+   * and boundary integrals are computed.
+   *
+   * The two additional functors work on a range of degrees of freedom,
+   * expressed in terms of the degree-of-freedom numbering of the selected
+   * DoFHandler `dof_handler_index_pre_post` in MPI-local indices. The
+   * arguments to the functors represent a range of degrees of freedom at a
+   * granularity of
+   * internal::MatrixFreeFunctions::DoFInfo::chunk_size_zero_vector entries
+   * (except for the last chunk which is set to the number of locally owned
+   * entries) in the form `[first, last)`. The idea of these functors is to
+   * bring operations on vectors closer to the point where they accessed in a
+   * matrix-free loop, with the goal to increase cache hits by temporal
+   * locality. This loop guarantees that the `operation_before_loop` hits all
+   * relevant unknowns before they are first touched by any of the cell, face or
+   * boundary operations (including the MPI data exchange), allowing to execute
+   * some vector update that the `src` vector depends upon. The
+   * `operation_after_loop` is similar - it starts to execute on a range of DoFs
+   * once all DoFs in that range have been touched for the last time by the
+   * cell, face and boundary operations (including the MPI data exchange),
+   * allowing e.g. to compute some vector operations that depend on the result
+   * of the current cell loop in `dst` or want to modify `src`. The efficiency
+   * of caching depends on the numbering of the degrees of freedom because of
+   * the granularity of the ranges.
+   *
+   * @param cell_operation Pointer to member function of `CLASS` with the
+   * signature <tt>cell_operation (const MatrixFree<dim,Number> &, OutVector &,
+   * InVector &, std::pair<unsigned int,unsigned int> &)</tt> where the first
+   * argument passes the data of the calling class and the last argument
+   * defines the range of cells which should be worked on (typically more than
+   * one cell should be worked on in order to reduce overheads).
+   *
+   * @param face_operation Pointer to member function of `CLASS` with the
+   * signature <tt>face_operation (const MatrixFree<dim,Number> &, OutVector &,
+   * InVector &, std::pair<unsigned int,unsigned int> &)</tt> in analogy to
+   * `cell_operation`, but now the part associated to the work on interior
+   * faces. Note that the MatrixFree framework treats periodic faces as
+   * interior ones, so they will be assigned their correct neighbor after
+   * applying periodicity constraints within the face_operation calls.
+   *
+   * @param boundary_operation Pointer to member function of `CLASS` with the
+   * signature <tt>boundary_operation (const MatrixFree<dim,Number> &, OutVector
+   * &, InVector &, std::pair<unsigned int,unsigned int> &)</tt> in analogy to
+   * `cell_operation` and `face_operation`, but now the part associated to the
+   * work on boundary faces. Boundary faces are separated by their
+   * `boundary_id` and it is possible to query that id using
+   * MatrixFree::get_boundary_id(). Note that both interior and faces use the
+   * same numbering, and faces in the interior are assigned lower numbers than
+   * the boundary faces.
+   *
+   * @param owning_class The object which provides the `cell_operation`
+   * call. To be compatible with this interface, the class must allow to call
+   * `owning_class->cell_operation(...)`.
+   *
+   * @param dst Destination vector holding the result. If the vector is of
+   * type LinearAlgebra::distributed::Vector (or composite objects thereof
+   * such as LinearAlgebra::distributed::BlockVector), the loop calls
+   * LinearAlgebra::distributed::Vector::compress() at the end of the call
+   * internally. For other vectors, including parallel Trilinos or PETSc
+   * vectors, no such call is issued. Note that Trilinos/Epetra or PETSc
+   * vectors do currently not work in parallel because the present class uses
+   * MPI-local index addressing, as opposed to the global addressing implied
+   * by those external libraries.
+   *
+   * @param src Input vector. If the vector is of type
+   * LinearAlgebra::distributed::Vector (or composite objects thereof such as
+   * LinearAlgebra::distributed::BlockVector), the loop calls
+   * LinearAlgebra::distributed::Vector::update_ghost_values() at the start of
+   * the call internally to make sure all necessary data is locally
+   * available. Note, however, that the vector is reset to its original state
+   * at the end of the loop, i.e., if the vector was not ghosted upon entry of
+   * the loop, it will not be ghosted upon finishing the loop.
+   *
+   * @param operation_before_loop This functor can be used to perform an
+   * operation on entries of the `src` and `dst` vectors (or other vectors)
+   * before the operation on cells first touches a particular DoF according to
+   * the general description in the text above. This function is passed a
+   * range of the locally owned degrees of freedom on the selected
+   * `dof_handler_index_pre_post` (in MPI-local numbering).
+   *
+   * @param operation_after_loop This functor can be used to perform an
+   * operation on entries of the `src` and `dst` vectors (or other vectors)
+   * after the operation on cells last touches a particular DoF according to
+   * the general description in the text above. This function is passed a
+   * range of the locally owned degrees of freedom on the selected
+   * `dof_handler_index_pre_post` (in MPI-local numbering).
+   *
+   * @param dof_handler_index_pre_post Since MatrixFree can be initialized
+   * with a vector of DoFHandler objects, each of them will in general have
+   * different vector sizes and thus different ranges returned to
+   * `operation_before_loop` and `operation_after_loop`. Use this variable to
+   * specify which one of the DoFHandler objects the index range should be
+   * associated to. Defaults to the `dof_handler_index` 0.
+   *
+   * @param dst_vector_face_access Set the type of access into the vector
+   * `dst` that will happen inside the body of the @p face_operation
+   * function. As explained in the description of the DataAccessOnFaces
+   * struct, the purpose of this selection is to reduce the amount of data
+   * that must be exchanged over the MPI network (or via `memcpy` if within
+   * the shared memory region of a node) to gain performance. Note that there
+   * is no way to communicate this setting with the FEFaceEvaluation class,
+   * therefore this selection must be made at this site in addition to what is
+   * implemented inside the `face_operation` function. As a consequence, there
+   * is also no way to check that the setting passed to this call is
+   * consistent with what is later done by `FEFaceEvaluation`, and it is the
+   * user's responsibility to ensure correctness of data.
+   *
+   * @param src_vector_face_access Set the type of access into the vector
+   * `src` that will happen inside the body of the @p face_operation function,
+   * in analogy to `dst_vector_face_access`.
+   *
+   * @note The close locality of the `operation_before_loop` and
+   * `operation_after_loop` is currently only implemented for the MPI-only
+   * case. In case threading is enabled, the complete `operation_before_loop`
+   * is scheduled before the parallel loop, and `operation_after_loop` is
+   * scheduled strictly afterwards, due to the complicated dependencies.
+   */
+  template <typename CLASS, typename OutVector, typename InVector>
+  void
+  loop(
+    void (CLASS::*cell_operation)(const MatrixFree &,
+                                  OutVector &,
+                                  const InVector &,
+                                  const std::pair<unsigned int, unsigned int> &)
+      const,
+    void (CLASS::*face_operation)(const MatrixFree &,
+                                  OutVector &,
+                                  const InVector &,
+                                  const std::pair<unsigned int, unsigned int> &)
+      const,
+    void (CLASS::*boundary_operation)(
+      const MatrixFree &,
+      OutVector &,
+      const InVector &,
+      const std::pair<unsigned int, unsigned int> &) const,
+    const CLASS *   owning_class,
+    OutVector &     dst,
+    const InVector &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 unsigned int      dof_handler_index_pre_post = 0,
+    const DataAccessOnFaces dst_vector_face_access =
+      DataAccessOnFaces::unspecified,
+    const DataAccessOnFaces src_vector_face_access =
+      DataAccessOnFaces::unspecified) const;
+
+  /**
+   * Same as above, but for class member functions which are non-const.
+   */
+  template <typename CLASS, typename OutVector, typename InVector>
+  void
+  loop(void (CLASS::*cell_operation)(
+         const MatrixFree &,
+         OutVector &,
+         const InVector &,
+         const std::pair<unsigned int, unsigned int> &),
+       void (CLASS::*face_operation)(
+         const MatrixFree &,
+         OutVector &,
+         const InVector &,
+         const std::pair<unsigned int, unsigned int> &),
+       void (CLASS::*boundary_operation)(
+         const MatrixFree &,
+         OutVector &,
+         const InVector &,
+         const std::pair<unsigned int, unsigned int> &),
+       const CLASS *   owning_class,
+       OutVector &     dst,
+       const InVector &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 unsigned int      dof_handler_index_pre_post = 0,
+       const DataAccessOnFaces dst_vector_face_access =
+         DataAccessOnFaces::unspecified,
+       const DataAccessOnFaces src_vector_face_access =
+         DataAccessOnFaces::unspecified) const;
+
+  /**
+   * Same as above, but taking an `std::function` as the `cell_operation`,
+   * `face_operation` and `boundary_operation` rather than a class member
+   * function.
+   */
+  template <typename OutVector, typename InVector>
+  void
+  loop(const std::function<
+         void(const MatrixFree<dim, Number, VectorizedArrayType> &,
+              OutVector &,
+              const InVector &,
+              const std::pair<unsigned int, unsigned int> &)> &cell_operation,
+       const std::function<
+         void(const MatrixFree<dim, Number, VectorizedArrayType> &,
+              OutVector &,
+              const InVector &,
+              const std::pair<unsigned int, unsigned int> &)> &face_operation,
+       const std::function<void(
+         const MatrixFree<dim, Number, VectorizedArrayType> &,
+         OutVector &,
+         const InVector &,
+         const std::pair<unsigned int, unsigned int> &)> &boundary_operation,
+       OutVector &                                        dst,
+       const InVector &                                   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 unsigned int      dof_handler_index_pre_post = 0,
+       const DataAccessOnFaces dst_vector_face_access =
+         DataAccessOnFaces::unspecified,
+       const DataAccessOnFaces src_vector_face_access =
+         DataAccessOnFaces::unspecified) const;
+
   /**
    * This method runs the loop over all cells (in parallel) similarly as
    * cell_loop() does. However, this function is intended to be used
@@ -5081,6 +5298,168 @@ MatrixFree<dim, Number, VectorizedArrayType>::loop(
 
 
 
+template <int dim, typename Number, typename VectorizedArrayType>
+template <typename OutVector, typename InVector>
+inline void
+MatrixFree<dim, Number, VectorizedArrayType>::loop(
+  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
+                           OutVector &,
+                           const InVector &,
+                           const std::pair<unsigned int, unsigned int> &)>
+    &cell_operation,
+  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
+                           OutVector &,
+                           const InVector &,
+                           const std::pair<unsigned int, unsigned int> &)>
+    &face_operation,
+  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
+                           OutVector &,
+                           const InVector &,
+                           const std::pair<unsigned int, unsigned int> &)>
+    &             boundary_operation,
+  OutVector &     dst,
+  const InVector &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 unsigned int      dof_handler_index_pre_post,
+  const DataAccessOnFaces dst_vector_face_access,
+  const DataAccessOnFaces src_vector_face_access) const
+{
+  using Wrapper =
+    internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
+                             InVector,
+                             OutVector>;
+  Wrapper wrap(cell_operation, face_operation, boundary_operation);
+  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
+                     InVector,
+                     OutVector,
+                     Wrapper,
+                     true>
+    worker(*this,
+           src,
+           dst,
+           false,
+           wrap,
+           &Wrapper::cell_integrator,
+           &Wrapper::face_integrator,
+           &Wrapper::boundary_integrator,
+           src_vector_face_access,
+           dst_vector_face_access,
+           operation_before_loop,
+           operation_after_loop,
+           dof_handler_index_pre_post);
+
+  task_info.loop(worker);
+}
+
+
+
+template <int dim, typename Number, typename VectorizedArrayType>
+template <typename CLASS, typename OutVector, typename InVector>
+inline void
+MatrixFree<dim, Number, VectorizedArrayType>::loop(
+  void (CLASS::*cell_operation)(const MatrixFree &,
+                                OutVector &,
+                                const InVector &,
+                                const std::pair<unsigned int, unsigned int> &)
+    const,
+  void (CLASS::*face_operation)(const MatrixFree &,
+                                OutVector &,
+                                const InVector &,
+                                const std::pair<unsigned int, unsigned int> &)
+    const,
+  void (CLASS::*boundary_operation)(
+    const MatrixFree &,
+    OutVector &,
+    const InVector &,
+    const std::pair<unsigned int, unsigned int> &) const,
+  const CLASS *   owning_class,
+  OutVector &     dst,
+  const InVector &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 unsigned int      dof_handler_index_pre_post,
+  const DataAccessOnFaces dst_vector_face_access,
+  const DataAccessOnFaces src_vector_face_access) const
+{
+  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
+                     InVector,
+                     OutVector,
+                     CLASS,
+                     true>
+    worker(*this,
+           src,
+           dst,
+           false,
+           *owning_class,
+           cell_operation,
+           face_operation,
+           boundary_operation,
+           src_vector_face_access,
+           dst_vector_face_access,
+           operation_before_loop,
+           operation_after_loop,
+           dof_handler_index_pre_post);
+  task_info.loop(worker);
+}
+
+
+
+template <int dim, typename Number, typename VectorizedArrayType>
+template <typename CLASS, typename OutVector, typename InVector>
+inline void
+MatrixFree<dim, Number, VectorizedArrayType>::loop(
+  void (CLASS::*cell_operation)(const MatrixFree &,
+                                OutVector &,
+                                const InVector &,
+                                const std::pair<unsigned int, unsigned int> &),
+  void (CLASS::*face_operation)(const MatrixFree &,
+                                OutVector &,
+                                const InVector &,
+                                const std::pair<unsigned int, unsigned int> &),
+  void (CLASS::*boundary_operation)(
+    const MatrixFree &,
+    OutVector &,
+    const InVector &,
+    const std::pair<unsigned int, unsigned int> &),
+  const CLASS *   owning_class,
+  OutVector &     dst,
+  const InVector &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 unsigned int      dof_handler_index_pre_post,
+  const DataAccessOnFaces dst_vector_face_access,
+  const DataAccessOnFaces src_vector_face_access) const
+{
+  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
+                     InVector,
+                     OutVector,
+                     CLASS,
+                     false>
+    worker(*this,
+           src,
+           dst,
+           false,
+           *owning_class,
+           cell_operation,
+           face_operation,
+           boundary_operation,
+           src_vector_face_access,
+           dst_vector_face_access,
+           operation_before_loop,
+           operation_after_loop,
+           dof_handler_index_pre_post);
+  task_info.loop(worker);
+}
+
+
+
 template <int dim, typename Number, typename VectorizedArrayType>
 template <typename CLASS, typename OutVector, typename InVector>
 inline void
diff --git a/tests/matrix_free/pre_and_post_loops_07.cc b/tests/matrix_free/pre_and_post_loops_07.cc
new file mode 100644 (file)
index 0000000..75fd553
--- /dev/null
@@ -0,0 +1,224 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 - 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// test the correctness of the matrix-free loop with additional
+// operation_before_loop and operation_after_loop lambdas
+
+// similar to pre_and_post_loop_01 but uses loop instead of cell_loop with
+// do-nothing face and boundary operations
+
+#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)
+  {}
+
+  const std::function<void(const MatrixFree<dim, Number> &,
+                           LinearAlgebra::distributed::Vector<Number> &,
+                           const LinearAlgebra::distributed::Vector<Number> &,
+                           const std::pair<unsigned int, unsigned int> &)>
+    do_nothing = [](const MatrixFree<dim, Number> &,
+                    LinearAlgebra::distributed::Vector<Number> &,
+                    const LinearAlgebra::distributed::Vector<Number> &,
+                    const std::pair<unsigned int, unsigned int> &) {};
+
+  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.loop(wrap, do_nothing, do_nothing, 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.loop(
+      wrap,
+      do_nothing,
+      do_nothing,
+      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;
+    data.mapping_update_flags_inner_faces = update_JxW_values;
+    mf_data.reinit(MappingQ1<dim>{}, 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)
+    {
+      // Multiply by 0.01 to make float error with roundoff less than the
+      // numdiff absolute tolerance
+      double entry          = 0.01 * random_value<double>();
+      vec1.local_element(i) = entry;
+      entry                 = 0.01 * random_value<double>();
+      vec2.local_element(i) = entry;
+      entry                 = 0.01 * 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_merged(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_merged(vec1, vec2, vec3);
+
+  mf.vmult_with_update_basic(ref1, ref2, ref3);
+  mf.vmult_with_update_merged(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_07.with_p4est=true.mpirun=1.output b/tests/matrix_free/pre_and_post_loops_07.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_07.with_p4est=true.mpirun=3.output b/tests/matrix_free/pre_and_post_loops_07.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_08.cc b/tests/matrix_free/pre_and_post_loops_08.cc
new file mode 100644 (file)
index 0000000..21ddadb
--- /dev/null
@@ -0,0 +1,198 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 - 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// perform some boundary and face operation but fully overwrite the result
+// within the post operation
+
+#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
+  do_nothing(const MatrixFree<dim, Number> &,
+             LinearAlgebra::distributed::Vector<Number> &,
+             const LinearAlgebra::distributed::Vector<Number> &,
+             const std::pair<unsigned int, unsigned int> &) const {};
+
+  void
+  do_boundary(const MatrixFree<dim, Number> &                   mf,
+              LinearAlgebra::distributed::Vector<Number> &      dst,
+              const LinearAlgebra::distributed::Vector<Number> &src,
+              const std::pair<unsigned int, unsigned int> &     range) const
+  {
+    using FEFaceIntegrator = FEFaceEvaluation<dim, -1, 0, 1, Number>;
+    FEFaceIntegrator eval(mf, true);
+
+    for (unsigned int face = range.first; face < range.second; ++face)
+      {
+        eval.reinit(face);
+        eval.gather_evaluate(src, EvaluationFlags::values);
+        for (unsigned int q = 0; q < eval.n_q_points; ++q)
+          {
+            eval.submit_value(1.23 * eval.get_value(q), q);
+          }
+        eval.integrate_scatter(EvaluationFlags::values, dst);
+      }
+  };
+
+  void
+  do_faces(const MatrixFree<dim, Number> &                   mf,
+           LinearAlgebra::distributed::Vector<Number> &      dst,
+           const LinearAlgebra::distributed::Vector<Number> &src,
+           const std::pair<unsigned int, unsigned int> &     range) const
+  {
+    using FEFaceIntegrator = FEFaceEvaluation<dim, -1, 0, 1, Number>;
+    FEFaceIntegrator eval(mf, true);
+
+    for (unsigned int face = range.first; face < range.second; ++face)
+      {
+        eval.reinit(face);
+        eval.gather_evaluate(src, EvaluationFlags::values);
+        for (unsigned int q = 0; q < eval.n_q_points; ++q)
+          {
+            eval.submit_value(1.23 * eval.get_value(q), q);
+          }
+        eval.integrate_scatter(EvaluationFlags::values, dst);
+      }
+  };
+
+  void
+  vmult_pre_post_op_only(
+    LinearAlgebra::distributed::Vector<Number> &      dst,
+    const LinearAlgebra::distributed::Vector<Number> &src) const
+  {
+    data.loop(
+      &Matrix::do_nothing,
+      &Matrix::do_faces,
+      &Matrix::do_boundary,
+      this,
+      dst,
+      src,
+      // operation before
+      [&](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;
+          }
+      },
+      // operation after
+      [&](const unsigned int start_range, const unsigned int end_range) {
+        for (unsigned int i = start_range; i < end_range; ++i)
+          {
+            dst.local_element(i) = 2 * src.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;
+    data.mapping_update_flags_inner_faces = update_JxW_values;
+    mf_data.reinit(MappingQ1<dim>{}, 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);
+
+  for (unsigned int i = 0; i < vec1.local_size(); ++i)
+    {
+      // Multiply by 0.01 to make float error with roundoff less than the
+      // numdiff absolute tolerance
+      double entry          = 0.01 * random_value<double>();
+      vec1.local_element(i) = entry;
+      entry                 = 0.01 * random_value<double>();
+      vec2.local_element(i) = entry;
+    }
+
+  mf.vmult_pre_post_op_only(vec2, vec1);
+
+  auto ref = vec1;
+  ref *= 2.0;
+  vec2 -= ref;
+  deallog << "Error: " << vec2.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_08.with_p4est=true.mpirun=2.output b/tests/matrix_free/pre_and_post_loops_08.with_p4est=true.mpirun=2.output
new file mode 100644 (file)
index 0000000..af278fb
--- /dev/null
@@ -0,0 +1,7 @@
+
+DEAL::Testing FE_Q<2>(3)
+DEAL::Error: 0.00
+DEAL::Testing FE_Q<2>(3)
+DEAL::Error: 0.00
+DEAL::Testing FE_Q<3>(3)
+DEAL::Error: 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.