From: Sebastian Proell Date: Sat, 8 Oct 2022 11:25:12 +0000 (+0200) Subject: MatrixFree pre and post operation for loop() X-Git-Tag: v9.5.0-rc1~900^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=6d18ce85efe4cb711bd3703b66472c8f5da4058d;p=dealii.git MatrixFree pre and post operation for loop() --- diff --git a/doc/news/changes/minor/20221011Proell b/doc/news/changes/minor/20221011Proell new file mode 100644 index 0000000000..11ea1f8245 --- /dev/null +++ b/doc/news/changes/minor/20221011Proell @@ -0,0 +1,5 @@ +New: A new overload for MatrixFree::loop allows +to specify a pre- and post-operation, similar +to MatrixFree::cell_loop. +
+(Sebastian Proell, 2022/10/11) diff --git a/include/deal.II/matrix_free/matrix_free.h b/include/deal.II/matrix_free/matrix_free.h index 39fa06c6b4..aba89df075 100644 --- a/include/deal.II/matrix_free/matrix_free.h +++ b/include/deal.II/matrix_free/matrix_free.h @@ -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 cell_operation (const MatrixFree &, OutVector &, + * InVector &, std::pair &) 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 face_operation (const MatrixFree &, OutVector &, + * InVector &, std::pair &) 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 boundary_operation (const MatrixFree &, OutVector + * &, InVector &, std::pair &) 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 + void + loop( + void (CLASS::*cell_operation)(const MatrixFree &, + OutVector &, + const InVector &, + const std::pair &) + const, + void (CLASS::*face_operation)(const MatrixFree &, + OutVector &, + const InVector &, + const std::pair &) + const, + void (CLASS::*boundary_operation)( + const MatrixFree &, + OutVector &, + const InVector &, + const std::pair &) const, + const CLASS * owning_class, + OutVector & dst, + const InVector &src, + const std::function + &operation_before_loop, + const std::function + & 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 + void + loop(void (CLASS::*cell_operation)( + const MatrixFree &, + OutVector &, + const InVector &, + const std::pair &), + void (CLASS::*face_operation)( + const MatrixFree &, + OutVector &, + const InVector &, + const std::pair &), + void (CLASS::*boundary_operation)( + const MatrixFree &, + OutVector &, + const InVector &, + const std::pair &), + const CLASS * owning_class, + OutVector & dst, + const InVector &src, + const std::function + &operation_before_loop, + const std::function + & 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 + void + loop(const std::function< + void(const MatrixFree &, + OutVector &, + const InVector &, + const std::pair &)> &cell_operation, + const std::function< + void(const MatrixFree &, + OutVector &, + const InVector &, + const std::pair &)> &face_operation, + const std::function &, + OutVector &, + const InVector &, + const std::pair &)> &boundary_operation, + OutVector & dst, + const InVector & src, + const std::function + &operation_before_loop, + const std::function + & 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::loop( +template +template +inline void +MatrixFree::loop( + const std::function &, + OutVector &, + const InVector &, + const std::pair &)> + &cell_operation, + const std::function &, + OutVector &, + const InVector &, + const std::pair &)> + &face_operation, + const std::function &, + OutVector &, + const InVector &, + const std::pair &)> + & boundary_operation, + OutVector & dst, + const InVector &src, + const std::function + &operation_before_loop, + const std::function + & 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, + InVector, + OutVector>; + Wrapper wrap(cell_operation, face_operation, boundary_operation); + internal::MFWorker, + 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 +template +inline void +MatrixFree::loop( + void (CLASS::*cell_operation)(const MatrixFree &, + OutVector &, + const InVector &, + const std::pair &) + const, + void (CLASS::*face_operation)(const MatrixFree &, + OutVector &, + const InVector &, + const std::pair &) + const, + void (CLASS::*boundary_operation)( + const MatrixFree &, + OutVector &, + const InVector &, + const std::pair &) const, + const CLASS * owning_class, + OutVector & dst, + const InVector &src, + const std::function + &operation_before_loop, + const std::function + & 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, + 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 +template +inline void +MatrixFree::loop( + void (CLASS::*cell_operation)(const MatrixFree &, + OutVector &, + const InVector &, + const std::pair &), + void (CLASS::*face_operation)(const MatrixFree &, + OutVector &, + const InVector &, + const std::pair &), + void (CLASS::*boundary_operation)( + const MatrixFree &, + OutVector &, + const InVector &, + const std::pair &), + const CLASS * owning_class, + OutVector & dst, + const InVector &src, + const std::function + &operation_before_loop, + const std::function + & 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, + 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 template 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 index 0000000000..75fd553524 --- /dev/null +++ b/tests/matrix_free/pre_and_post_loops_07.cc @@ -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 + +#include +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include + +#include "../tests.h" + +#include "matrix_vector_mf.h" + + + +template +class Matrix +{ +public: + Matrix(const MatrixFree &data_in) + : data(data_in) + {} + + const std::function &, + LinearAlgebra::distributed::Vector &, + const LinearAlgebra::distributed::Vector &, + const std::pair &)> + do_nothing = [](const MatrixFree &, + LinearAlgebra::distributed::Vector &, + const LinearAlgebra::distributed::Vector &, + const std::pair &) {}; + + void + vmult(LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src) const + { + const std::function &, + LinearAlgebra::distributed::Vector &, + const LinearAlgebra::distributed::Vector &, + const std::pair &)> + wrap = helmholtz_operator, + 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 &dst, + LinearAlgebra::distributed::Vector &src, + LinearAlgebra::distributed::Vector &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 &dst, + LinearAlgebra::distributed::Vector &src, + LinearAlgebra::distributed::Vector &other_vector) const + { + const std::function &, + LinearAlgebra::distributed::Vector &, + const LinearAlgebra::distributed::Vector &, + const std::pair &)> + wrap = helmholtz_operator, + 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 &data; +}; + + + +template +void +test() +{ + parallel::distributed::Triangulation tria(MPI_COMM_WORLD); + GridGenerator::hyper_cube(tria); + tria.refine_global(7 - dim); + + FE_Q fe(fe_degree); + DoFHandler dof(tria); + dof.distribute_dofs(fe); + AffineConstraints constraints; + constraints.close(); + + deallog << "Testing " << dof.get_fe().get_name() << std::endl; + + MatrixFree mf_data; + { + const QGauss<1> quad(fe_degree + 1); + typename MatrixFree::AdditionalData data; + data.tasks_parallel_scheme = MatrixFree::AdditionalData::none; + data.mapping_update_flags_inner_faces = update_JxW_values; + mf_data.reinit(MappingQ1{}, dof, constraints, quad, data); + } + + Matrix mf(mf_data); + LinearAlgebra::distributed::Vector 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(); + vec1.local_element(i) = entry; + entry = 0.01 * random_value(); + vec2.local_element(i) = entry; + entry = 0.01 * random_value(); + vec3.local_element(i) = entry; + } + + LinearAlgebra::distributed::Vector ref1 = vec1; + LinearAlgebra::distributed::Vector ref2 = vec2; + LinearAlgebra::distributed::Vector 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 index 0000000000..57a015933a --- /dev/null +++ b/tests/matrix_free/pre_and_post_loops_07.with_p4est=true.mpirun=1.output @@ -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 index 0000000000..57a015933a --- /dev/null +++ b/tests/matrix_free/pre_and_post_loops_07.with_p4est=true.mpirun=3.output @@ -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 index 0000000000..21ddadb4ff --- /dev/null +++ b/tests/matrix_free/pre_and_post_loops_08.cc @@ -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 + +#include +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include + +#include "../tests.h" + +#include "matrix_vector_mf.h" + + + +template +class Matrix +{ +public: + Matrix(const MatrixFree &data_in) + : data(data_in) + {} + + void + do_nothing(const MatrixFree &, + LinearAlgebra::distributed::Vector &, + const LinearAlgebra::distributed::Vector &, + const std::pair &) const {}; + + void + do_boundary(const MatrixFree & mf, + LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src, + const std::pair & range) const + { + using FEFaceIntegrator = FEFaceEvaluation; + 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 & mf, + LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src, + const std::pair & range) const + { + using FEFaceIntegrator = FEFaceEvaluation; + 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 & dst, + const LinearAlgebra::distributed::Vector &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 &data; +}; + + + +template +void +test() +{ + parallel::distributed::Triangulation tria(MPI_COMM_WORLD); + GridGenerator::hyper_cube(tria); + tria.refine_global(7 - dim); + + FE_Q fe(fe_degree); + DoFHandler dof(tria); + dof.distribute_dofs(fe); + AffineConstraints constraints; + constraints.close(); + + deallog << "Testing " << dof.get_fe().get_name() << std::endl; + + MatrixFree mf_data; + { + const QGauss<1> quad(fe_degree + 1); + typename MatrixFree::AdditionalData data; + data.tasks_parallel_scheme = MatrixFree::AdditionalData::none; + data.mapping_update_flags_inner_faces = update_JxW_values; + mf_data.reinit(MappingQ1{}, dof, constraints, quad, data); + } + + Matrix mf(mf_data); + LinearAlgebra::distributed::Vector 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(); + vec1.local_element(i) = entry; + entry = 0.01 * random_value(); + 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 index 0000000000..af278fb74e --- /dev/null +++ b/tests/matrix_free/pre_and_post_loops_08.with_p4est=true.mpirun=2.output @@ -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