matrix_free.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2011 - 2018 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.
//
// ---------------------------------------------------------------------


#ifndef dealii_matrix_free_h
#define dealii_matrix_free_h

#include <deal.II/base/aligned_vector.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/quadrature.h>
#include <deal.II/base/template_constraints.h>
#include <deal.II/base/thread_local_storage.h>
#include <deal.II/base/vectorization.h>

#include <deal.II/dofs/dof_handler.h>

#include <deal.II/fe/fe.h>
#include <deal.II/fe/mapping.h>
#include <deal.II/fe/mapping_q1.h>

#include <deal.II/grid/grid_tools.h>

#include <deal.II/hp/dof_handler.h>
#include <deal.II/hp/q_collection.h>

#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/block_vector_base.h>
#include <deal.II/lac/la_parallel_vector.h>
#include <deal.II/lac/vector_operation.h>

#include <deal.II/matrix_free/dof_info.h>
#include <deal.II/matrix_free/mapping_info.h>
#include <deal.II/matrix_free/shape_info.h>
#include <deal.II/matrix_free/task_info.h>

#include <cstdlib>
#include <limits>
#include <list>
#include <memory>


DEAL_II_NAMESPACE_OPEN



template <int dim, typename Number = double>
class MatrixFree : public Subscriptor
{
public:
  using value_type = Number;

  static const unsigned int dimension = dim;

  struct AdditionalData
  {
    enum TasksParallelScheme
    {
      none = internal::MatrixFreeFunctions::TaskInfo::none,
      partition_partition =
        internal::MatrixFreeFunctions::TaskInfo::partition_partition,
      partition_color =
        internal::MatrixFreeFunctions::TaskInfo::partition_color,
      color = internal::MatrixFreeFunctions::TaskInfo::color
    };

    AdditionalData(
      const TasksParallelScheme tasks_parallel_scheme = partition_partition,
      const unsigned int        tasks_block_size      = 0,
      const UpdateFlags         mapping_update_flags  = update_gradients |
                                               update_JxW_values,
      const UpdateFlags  mapping_update_flags_boundary_faces = update_default,
      const UpdateFlags  mapping_update_flags_inner_faces    = update_default,
      const UpdateFlags  mapping_update_flags_faces_by_cells = update_default,
      const unsigned int level_mg_handler    = numbers::invalid_unsigned_int,
      const bool         store_plain_indices = true,
      const bool         initialize_indices  = true,
      const bool         initialize_mapping  = true,
      const bool         overlap_communication_computation    = true,
      const bool         hold_all_faces_to_owned_cells        = false,
      const bool         cell_vectorization_categories_strict = false)
      : tasks_parallel_scheme(tasks_parallel_scheme)
      , tasks_block_size(tasks_block_size)
      , mapping_update_flags(mapping_update_flags)
      , mapping_update_flags_boundary_faces(mapping_update_flags_boundary_faces)
      , mapping_update_flags_inner_faces(mapping_update_flags_inner_faces)
      , mapping_update_flags_faces_by_cells(mapping_update_flags_faces_by_cells)
      , level_mg_handler(level_mg_handler)
      , store_plain_indices(store_plain_indices)
      , initialize_indices(initialize_indices)
      , initialize_mapping(initialize_mapping)
      , overlap_communication_computation(overlap_communication_computation)
      , hold_all_faces_to_owned_cells(hold_all_faces_to_owned_cells)
      , cell_vectorization_categories_strict(
          cell_vectorization_categories_strict)
    {}

    TasksParallelScheme tasks_parallel_scheme;

    unsigned int tasks_block_size;

    UpdateFlags mapping_update_flags;

    UpdateFlags mapping_update_flags_boundary_faces;

    UpdateFlags mapping_update_flags_inner_faces;

    UpdateFlags mapping_update_flags_faces_by_cells;

    unsigned int level_mg_handler;

    bool store_plain_indices;

    bool initialize_indices;

    bool initialize_mapping;

    bool overlap_communication_computation;

    bool hold_all_faces_to_owned_cells;

    std::vector<unsigned int> cell_vectorization_category;

    bool cell_vectorization_categories_strict;
  };

  MatrixFree();

  MatrixFree(const MatrixFree<dim, Number> &other);

  ~MatrixFree() override = default;

  template <typename DoFHandlerType, typename QuadratureType, typename number2>
  void
  reinit(const Mapping<dim> &              mapping,
         const DoFHandlerType &            dof_handler,
         const AffineConstraints<number2> &constraint,
         const QuadratureType &            quad,
         const AdditionalData              additional_data = AdditionalData());

  template <typename DoFHandlerType, typename QuadratureType, typename number2>
  void
  reinit(const DoFHandlerType &            dof_handler,
         const AffineConstraints<number2> &constraint,
         const QuadratureType &            quad,
         const AdditionalData              additional_data = AdditionalData());

  template <typename DoFHandlerType, typename QuadratureType, typename number2>
  DEAL_II_DEPRECATED void
  reinit(const Mapping<dim> &              mapping,
         const DoFHandlerType &            dof_handler,
         const AffineConstraints<number2> &constraint,
         const IndexSet &                  locally_owned_dofs,
         const QuadratureType &            quad,
         const AdditionalData              additional_data = AdditionalData());

  template <typename DoFHandlerType, typename QuadratureType, typename number2>
  void
  reinit(const Mapping<dim> &                                   mapping,
         const std::vector<const DoFHandlerType *> &            dof_handler,
         const std::vector<const AffineConstraints<number2> *> &constraint,
         const std::vector<QuadratureType> &                    quad,
         const AdditionalData additional_data = AdditionalData());

  template <typename DoFHandlerType, typename QuadratureType, typename number2>
  void
  reinit(const std::vector<const DoFHandlerType *> &            dof_handler,
         const std::vector<const AffineConstraints<number2> *> &constraint,
         const std::vector<QuadratureType> &                    quad,
         const AdditionalData additional_data = AdditionalData());

  template <typename DoFHandlerType, typename QuadratureType, typename number2>
  DEAL_II_DEPRECATED void
  reinit(const Mapping<dim> &                                   mapping,
         const std::vector<const DoFHandlerType *> &            dof_handler,
         const std::vector<const AffineConstraints<number2> *> &constraint,
         const std::vector<IndexSet> &      locally_owned_set,
         const std::vector<QuadratureType> &quad,
         const AdditionalData               additional_data = AdditionalData());

  template <typename DoFHandlerType, typename QuadratureType, typename number2>
  void
  reinit(const Mapping<dim> &                                   mapping,
         const std::vector<const DoFHandlerType *> &            dof_handler,
         const std::vector<const AffineConstraints<number2> *> &constraint,
         const QuadratureType &                                 quad,
         const AdditionalData additional_data = AdditionalData());

  template <typename DoFHandlerType, typename QuadratureType, typename number2>
  void
  reinit(const std::vector<const DoFHandlerType *> &            dof_handler,
         const std::vector<const AffineConstraints<number2> *> &constraint,
         const QuadratureType &                                 quad,
         const AdditionalData additional_data = AdditionalData());

  void
  copy_from(const MatrixFree<dim, Number> &matrix_free_base);

  void
  clear();


  enum class DataAccessOnFaces
  {
    none,

    values,

    gradients,

    unspecified
  };

  template <typename OutVector, typename InVector>
  void
  cell_loop(
    const std::function<void(const MatrixFree<dim, Number> &,
                             OutVector &,
                             const InVector &,
                             const std::pair<unsigned int, unsigned int> &)>
      &             cell_operation,
    OutVector &     dst,
    const InVector &src,
    const bool      zero_dst_vector = false) const;

  template <typename CLASS, typename OutVector, typename InVector>
  void
  cell_loop(void (CLASS::*cell_operation)(
              const MatrixFree &,
              OutVector &,
              const InVector &,
              const std::pair<unsigned int, unsigned int> &) const,
            const CLASS *   owning_class,
            OutVector &     dst,
            const InVector &src,
            const bool      zero_dst_vector = false) const;

  template <typename CLASS, typename OutVector, typename InVector>
  void
  cell_loop(void (CLASS::*cell_operation)(
              const MatrixFree &,
              OutVector &,
              const InVector &,
              const std::pair<unsigned int, unsigned int> &),
            CLASS *         owning_class,
            OutVector &     dst,
            const InVector &src,
            const bool      zero_dst_vector = false) const;

  template <typename OutVector, typename InVector>
  void
  loop(const std::function<void(const MatrixFree<dim, Number> &,
                                OutVector &,
                                const InVector &,
                                const std::pair<unsigned int, unsigned int> &)>
         &cell_operation,
       const std::function<void(const MatrixFree<dim, Number> &,
                                OutVector &,
                                const InVector &,
                                const std::pair<unsigned int, unsigned int> &)>
         &face_operation,
       const std::function<void(const MatrixFree<dim, Number> &,
                                OutVector &,
                                const InVector &,
                                const std::pair<unsigned int, unsigned int> &)>
         &                     boundary_operation,
       OutVector &             dst,
       const InVector &        src,
       const bool              zero_dst_vector = false,
       const DataAccessOnFaces dst_vector_face_access =
         DataAccessOnFaces::unspecified,
       const DataAccessOnFaces src_vector_face_access =
         DataAccessOnFaces::unspecified) 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> &)
      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 bool              zero_dst_vector = false,
    const DataAccessOnFaces dst_vector_face_access =
      DataAccessOnFaces::unspecified,
    const DataAccessOnFaces src_vector_face_access =
      DataAccessOnFaces::unspecified) 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> &),
       CLASS *                 owning_class,
       OutVector &             dst,
       const InVector &        src,
       const bool              zero_dst_vector = false,
       const DataAccessOnFaces dst_vector_face_access =
         DataAccessOnFaces::unspecified,
       const DataAccessOnFaces src_vector_face_access =
         DataAccessOnFaces::unspecified) const;

  std::pair<unsigned int, unsigned int>
  create_cell_subrange_hp(const std::pair<unsigned int, unsigned int> &range,
                          const unsigned int fe_degree,
                          const unsigned int dof_handler_index = 0) const;

  std::pair<unsigned int, unsigned int>
  create_cell_subrange_hp_by_index(
    const std::pair<unsigned int, unsigned int> &range,
    const unsigned int                           fe_index,
    const unsigned int                           dof_handler_index = 0) const;


  template <typename VectorType>
  void
  initialize_dof_vector(VectorType &       vec,
                        const unsigned int dof_handler_index = 0) const;

  template <typename Number2>
  void
  initialize_dof_vector(LinearAlgebra::distributed::Vector<Number2> &vec,
                        const unsigned int dof_handler_index = 0) const;

  const std::shared_ptr<const Utilities::MPI::Partitioner> &
  get_vector_partitioner(const unsigned int dof_handler_index = 0) const;

  const IndexSet &
  get_locally_owned_set(const unsigned int dof_handler_index = 0) const;

  const IndexSet &
  get_ghost_set(const unsigned int dof_handler_index = 0) const;

  const std::vector<unsigned int> &
  get_constrained_dofs(const unsigned int dof_handler_index = 0) const;

  void
  renumber_dofs(std::vector<types::global_dof_index> &renumbering,
                const unsigned int                    dof_handler_index = 0);


  template <int spacedim>
  static bool
  is_supported(const FiniteElement<dim, spacedim> &fe);

  unsigned int
  n_components() const;

  unsigned int
  n_base_elements(const unsigned int dof_handler_index) const;

  unsigned int
  n_physical_cells() const;

  unsigned int
  n_macro_cells() const;

  unsigned int
  n_cell_batches() const;

  unsigned int
  n_ghost_cell_batches() const;

  unsigned int
  n_inner_face_batches() const;

  unsigned int
  n_boundary_face_batches() const;

  unsigned int
  n_ghost_inner_face_batches() const;

  types::boundary_id
  get_boundary_id(const unsigned int macro_face) const;

  std::array<types::boundary_id, VectorizedArray<Number>::n_array_elements>
  get_faces_by_cells_boundary_id(const unsigned int macro_cell,
                                 const unsigned int face_number) const;

  const DoFHandler<dim> &
  get_dof_handler(const unsigned int dof_handler_index = 0) const;

  typename DoFHandler<dim>::cell_iterator
  get_cell_iterator(const unsigned int macro_cell_number,
                    const unsigned int vector_number,
                    const unsigned int fe_component = 0) const;

  typename hp::DoFHandler<dim>::active_cell_iterator
  get_hp_cell_iterator(const unsigned int macro_cell_number,
                       const unsigned int vector_number,
                       const unsigned int dof_handler_index = 0) const;

  bool
  at_irregular_cell(const unsigned int macro_cell_number) const;

  unsigned int
  n_components_filled(const unsigned int cell_batch_number) const;

  unsigned int
  n_active_entries_per_cell_batch(const unsigned int cell_batch_number) const;

  unsigned int
  n_active_entries_per_face_batch(const unsigned int face_batch_number) const;

  unsigned int
  get_dofs_per_cell(const unsigned int dof_handler_index  = 0,
                    const unsigned int hp_active_fe_index = 0) const;

  unsigned int
  get_n_q_points(const unsigned int quad_index         = 0,
                 const unsigned int hp_active_fe_index = 0) const;

  unsigned int
  get_dofs_per_face(const unsigned int fe_component       = 0,
                    const unsigned int hp_active_fe_index = 0) const;

  unsigned int
  get_n_q_points_face(const unsigned int quad_index         = 0,
                      const unsigned int hp_active_fe_index = 0) const;

  const Quadrature<dim> &
  get_quadrature(const unsigned int quad_index         = 0,
                 const unsigned int hp_active_fe_index = 0) const;

  const Quadrature<dim - 1> &
  get_face_quadrature(const unsigned int quad_index         = 0,
                      const unsigned int hp_active_fe_index = 0) const;

  unsigned int
  get_cell_category(const unsigned int macro_cell) const;

  std::pair<unsigned int, unsigned int>
  get_face_category(const unsigned int macro_face) const;

  bool
  indices_initialized() const;

  bool
  mapping_initialized() const;

  std::size_t
  memory_consumption() const;

  template <typename StreamType>
  void
  print_memory_consumption(StreamType &out) const;

  void
  print(std::ostream &out) const;


  const internal::MatrixFreeFunctions::TaskInfo &
  get_task_info() const;

  DEAL_II_DEPRECATED
  const internal::MatrixFreeFunctions::TaskInfo &
  get_size_info() const;

  /*
   * Return geometry-dependent information on the cells.
   */
  const internal::MatrixFreeFunctions::MappingInfo<dim, Number> &
  get_mapping_info() const;

  const internal::MatrixFreeFunctions::DoFInfo &
  get_dof_info(const unsigned int dof_handler_index_component = 0) const;

  unsigned int
  n_constraint_pool_entries() const;

  const Number *
  constraint_pool_begin(const unsigned int pool_index) const;

  const Number *
  constraint_pool_end(const unsigned int pool_index) const;

  const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &
  get_shape_info(const unsigned int dof_handler_index_component = 0,
                 const unsigned int quad_index                  = 0,
                 const unsigned int fe_base_element             = 0,
                 const unsigned int hp_active_fe_index          = 0,
                 const unsigned int hp_active_quad_index        = 0) const;

  const internal::MatrixFreeFunctions::FaceToCellTopology<
    VectorizedArray<Number>::n_array_elements> &
  get_face_info(const unsigned int face_batch_number) const;

  AlignedVector<VectorizedArray<Number>> *
  acquire_scratch_data() const;

  void
  release_scratch_data(
    const AlignedVector<VectorizedArray<Number>> *memory) const;

  AlignedVector<Number> *
  acquire_scratch_data_non_threadsafe() const;

  void
  release_scratch_data_non_threadsafe(
    const AlignedVector<Number> *memory) const;


private:
  template <typename number2>
  void
  internal_reinit(
    const Mapping<dim> &                                   mapping,
    const std::vector<const DoFHandler<dim> *> &           dof_handler,
    const std::vector<const AffineConstraints<number2> *> &constraint,
    const std::vector<IndexSet> &                          locally_owned_set,
    const std::vector<hp::QCollection<1>> &                quad,
    const AdditionalData &                                 additional_data);

  template <typename number2>
  void
  internal_reinit(
    const Mapping<dim> &                                   mapping,
    const std::vector<const hp::DoFHandler<dim> *> &       dof_handler,
    const std::vector<const AffineConstraints<number2> *> &constraint,
    const std::vector<IndexSet> &                          locally_owned_set,
    const std::vector<hp::QCollection<1>> &                quad,
    const AdditionalData &                                 additional_data);

  template <typename number2>
  void
  initialize_indices(
    const std::vector<const AffineConstraints<number2> *> &constraint,
    const std::vector<IndexSet> &                          locally_owned_set,
    const AdditionalData &                                 additional_data);

  void
  initialize_dof_handlers(
    const std::vector<const DoFHandler<dim> *> &dof_handlers,
    const AdditionalData &                      additional_data);

  void
  initialize_dof_handlers(
    const std::vector<const hp::DoFHandler<dim> *> &dof_handlers,
    const AdditionalData &                          additional_data);

  void
  make_connectivity_graph_faces(DynamicSparsityPattern &connectivity);

  struct DoFHandlers
  {
    DoFHandlers()
      : active_dof_handler(usual)
      , n_dof_handlers(0)
      , level(numbers::invalid_unsigned_int)
    {}

    std::vector<SmartPointer<const DoFHandler<dim>>>     dof_handler;
    std::vector<SmartPointer<const hp::DoFHandler<dim>>> hp_dof_handler;
    enum ActiveDoFHandler
    {
      usual,
      hp
    } active_dof_handler;
    unsigned int n_dof_handlers;
    unsigned int level;
  };

  DoFHandlers dof_handlers;

  std::vector<internal::MatrixFreeFunctions::DoFInfo> dof_info;

  std::vector<Number> constraint_pool_data;

  std::vector<unsigned int> constraint_pool_row_index;

  internal::MatrixFreeFunctions::MappingInfo<dim, Number> mapping_info;

  Table<4, internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>>>
    shape_info;

  std::vector<std::pair<unsigned int, unsigned int>> cell_level_index;


  unsigned int cell_level_index_end_local;

  internal::MatrixFreeFunctions::TaskInfo task_info;

  internal::MatrixFreeFunctions::FaceInfo<
    VectorizedArray<Number>::n_array_elements>
    face_info;

  bool indices_are_initialized;

  bool mapping_is_initialized;

  mutable Threads::ThreadLocalStorage<
    std::list<std::pair<bool, AlignedVector<VectorizedArray<Number>>>>>
    scratch_pad;

  mutable std::list<std::pair<bool, AlignedVector<Number>>>
    scratch_pad_non_threadsafe;
};



/*----------------------- Inline functions ----------------------------------*/

#ifndef DOXYGEN



template <int dim, typename Number>
template <typename VectorType>
inline void
MatrixFree<dim, Number>::initialize_dof_vector(VectorType &       vec,
                                               const unsigned int comp) const
{
  AssertIndexRange(comp, n_components());
  vec.reinit(dof_info[comp].vector_partitioner->size());
}



template <int dim, typename Number>
template <typename Number2>
inline void
MatrixFree<dim, Number>::initialize_dof_vector(
  LinearAlgebra::distributed::Vector<Number2> &vec,
  const unsigned int                           comp) const
{
  AssertIndexRange(comp, n_components());
  vec.reinit(dof_info[comp].vector_partitioner);
}



template <int dim, typename Number>
inline const std::shared_ptr<const Utilities::MPI::Partitioner> &
MatrixFree<dim, Number>::get_vector_partitioner(const unsigned int comp) const
{
  AssertIndexRange(comp, n_components());
  return dof_info[comp].vector_partitioner;
}



template <int dim, typename Number>
inline const std::vector<unsigned int> &
MatrixFree<dim, Number>::get_constrained_dofs(const unsigned int comp) const
{
  AssertIndexRange(comp, n_components());
  return dof_info[comp].constrained_dofs;
}



template <int dim, typename Number>
inline unsigned int
MatrixFree<dim, Number>::n_components() const
{
  AssertDimension(dof_handlers.n_dof_handlers, dof_info.size());
  return dof_handlers.n_dof_handlers;
}



template <int dim, typename Number>
inline unsigned int
MatrixFree<dim, Number>::n_base_elements(const unsigned int dof_no) const
{
  AssertDimension(dof_handlers.n_dof_handlers, dof_info.size());
  AssertIndexRange(dof_no, dof_handlers.n_dof_handlers);
  return dof_handlers.dof_handler[dof_no]->get_fe().n_base_elements();
}



template <int dim, typename Number>
inline const internal::MatrixFreeFunctions::TaskInfo &
MatrixFree<dim, Number>::get_task_info() const
{
  return task_info;
}



template <int dim, typename Number>
inline const internal::MatrixFreeFunctions::TaskInfo &
MatrixFree<dim, Number>::get_size_info() const
{
  return task_info;
}



template <int dim, typename Number>
inline unsigned int
MatrixFree<dim, Number>::n_macro_cells() const
{
  return *(task_info.cell_partition_data.end() - 2);
}



template <int dim, typename Number>
inline unsigned int
MatrixFree<dim, Number>::n_physical_cells() const
{
  return task_info.n_active_cells;
}



template <int dim, typename Number>
inline unsigned int
MatrixFree<dim, Number>::n_cell_batches() const
{
  return *(task_info.cell_partition_data.end() - 2);
}



template <int dim, typename Number>
inline unsigned int
MatrixFree<dim, Number>::n_ghost_cell_batches() const
{
  return *(task_info.cell_partition_data.end() - 1) -
         *(task_info.cell_partition_data.end() - 2);
}



template <int dim, typename Number>
inline unsigned int
MatrixFree<dim, Number>::n_inner_face_batches() const
{
  if (task_info.face_partition_data.size() == 0)
    return 0;
  return task_info.face_partition_data.back();
}



template <int dim, typename Number>
inline unsigned int
MatrixFree<dim, Number>::n_boundary_face_batches() const
{
  if (task_info.face_partition_data.size() == 0)
    return 0;
  return task_info.boundary_partition_data.back() -
         task_info.face_partition_data.back();
}



template <int dim, typename Number>
inline unsigned int
MatrixFree<dim, Number>::n_ghost_inner_face_batches() const
{
  if (task_info.face_partition_data.size() == 0)
    return 0;
  return face_info.faces.size() - task_info.boundary_partition_data.back();
}



template <int dim, typename Number>
inline types::boundary_id
MatrixFree<dim, Number>::get_boundary_id(const unsigned int macro_face) const
{
  Assert(macro_face >= task_info.boundary_partition_data[0] &&
           macro_face < task_info.boundary_partition_data.back(),
         ExcIndexRange(macro_face,
                       task_info.boundary_partition_data[0],
                       task_info.boundary_partition_data.back()));
  return types::boundary_id(face_info.faces[macro_face].exterior_face_no);
}



template <int dim, typename Number>
inline std::array<types::boundary_id, VectorizedArray<Number>::n_array_elements>
MatrixFree<dim, Number>::get_faces_by_cells_boundary_id(
  const unsigned int macro_cell,
  const unsigned int face_number) const
{
  AssertIndexRange(macro_cell, n_macro_cells());
  AssertIndexRange(face_number, GeometryInfo<dim>::faces_per_cell);
  Assert(face_info.cell_and_face_boundary_id.size(0) >= n_macro_cells(),
         ExcNotInitialized());
  std::array<types::boundary_id, VectorizedArray<Number>::n_array_elements>
    result;
  result.fill(numbers::invalid_boundary_id);
  for (unsigned int v = 0; v < n_active_entries_per_cell_batch(macro_cell); ++v)
    result[v] = face_info.cell_and_face_boundary_id(macro_cell, face_number, v);
  return result;
}



template <int dim, typename Number>
inline const internal::MatrixFreeFunctions::MappingInfo<dim, Number> &
MatrixFree<dim, Number>::get_mapping_info() const
{
  return mapping_info;
}



template <int dim, typename Number>
inline const internal::MatrixFreeFunctions::DoFInfo &
MatrixFree<dim, Number>::get_dof_info(const unsigned int dof_index) const
{
  AssertIndexRange(dof_index, n_components());
  return dof_info[dof_index];
}



template <int dim, typename Number>
inline unsigned int
MatrixFree<dim, Number>::n_constraint_pool_entries() const
{
  return constraint_pool_row_index.size() - 1;
}



template <int dim, typename Number>
inline const Number *
MatrixFree<dim, Number>::constraint_pool_begin(const unsigned int row) const
{
  AssertIndexRange(row, constraint_pool_row_index.size() - 1);
  return constraint_pool_data.empty() ?
           nullptr :
           constraint_pool_data.data() + constraint_pool_row_index[row];
}



template <int dim, typename Number>
inline const Number *
MatrixFree<dim, Number>::constraint_pool_end(const unsigned int row) const
{
  AssertIndexRange(row, constraint_pool_row_index.size() - 1);
  return constraint_pool_data.empty() ?
           nullptr :
           constraint_pool_data.data() + constraint_pool_row_index[row + 1];
}



template <int dim, typename Number>
inline std::pair<unsigned int, unsigned int>
MatrixFree<dim, Number>::create_cell_subrange_hp(
  const std::pair<unsigned int, unsigned int> &range,
  const unsigned int                           degree,
  const unsigned int                           dof_handler_component) const
{
  if (dof_info[dof_handler_component].cell_active_fe_index.empty())
    {
      AssertDimension(
        dof_info[dof_handler_component].fe_index_conversion.size(), 1);
      AssertDimension(
        dof_info[dof_handler_component].fe_index_conversion[0].size(), 1);
      if (dof_info[dof_handler_component].fe_index_conversion[0][0] == degree)
        return range;
      else
        return std::pair<unsigned int, unsigned int>(range.second,
                                                     range.second);
    }

  const unsigned int fe_index =
    dof_info[dof_handler_component].fe_index_from_degree(0, degree);
  if (fe_index >= dof_info[dof_handler_component].max_fe_index)
    return std::pair<unsigned int, unsigned int>(range.second, range.second);
  else
    return create_cell_subrange_hp_by_index(range,
                                            fe_index,
                                            dof_handler_component);
}



template <int dim, typename Number>
inline bool
MatrixFree<dim, Number>::at_irregular_cell(const unsigned int macro_cell) const
{
  AssertIndexRange(macro_cell, task_info.cell_partition_data.back());
  return VectorizedArray<Number>::n_array_elements > 1 &&
         cell_level_index[(macro_cell + 1) *
                            VectorizedArray<Number>::n_array_elements -
                          1] ==
           cell_level_index[(macro_cell + 1) *
                              VectorizedArray<Number>::n_array_elements -
                            2];
}



template <int dim, typename Number>
inline unsigned int
MatrixFree<dim, Number>::n_components_filled(
  const unsigned int cell_batch_number) const
{
  return n_active_entries_per_cell_batch(cell_batch_number);
}



template <int dim, typename Number>
inline unsigned int
MatrixFree<dim, Number>::n_active_entries_per_cell_batch(
  const unsigned int cell_batch_number) const
{
  AssertIndexRange(cell_batch_number, task_info.cell_partition_data.back());
  unsigned int n_components = VectorizedArray<Number>::n_array_elements;
  while (n_components > 1 &&
         cell_level_index[cell_batch_number *
                            VectorizedArray<Number>::n_array_elements +
                          n_components - 1] ==
           cell_level_index[cell_batch_number *
                              VectorizedArray<Number>::n_array_elements +
                            n_components - 2])
    --n_components;
  AssertIndexRange(n_components - 1, VectorizedArray<Number>::n_array_elements);
  return n_components;
}



template <int dim, typename Number>
inline unsigned int
MatrixFree<dim, Number>::n_active_entries_per_face_batch(
  const unsigned int face_batch_number) const
{
  AssertIndexRange(face_batch_number, face_info.faces.size());
  unsigned int n_components = VectorizedArray<Number>::n_array_elements;
  while (n_components > 1 &&
         face_info.faces[face_batch_number].cells_interior[n_components - 1] ==
           numbers::invalid_unsigned_int)
    --n_components;
  AssertIndexRange(n_components - 1, VectorizedArray<Number>::n_array_elements);
  return n_components;
}



template <int dim, typename Number>
inline unsigned int
MatrixFree<dim, Number>::get_dofs_per_cell(
  const unsigned int dof_handler_index,
  const unsigned int active_fe_index) const
{
  return dof_info[dof_handler_index].dofs_per_cell[active_fe_index];
}



template <int dim, typename Number>
inline unsigned int
MatrixFree<dim, Number>::get_n_q_points(
  const unsigned int quad_index,
  const unsigned int active_fe_index) const
{
  AssertIndexRange(quad_index, mapping_info.cell_data.size());
  return mapping_info.cell_data[quad_index]
    .descriptor[active_fe_index]
    .n_q_points;
}



template <int dim, typename Number>
inline unsigned int
MatrixFree<dim, Number>::get_dofs_per_face(
  const unsigned int dof_handler_index,
  const unsigned int active_fe_index) const
{
  return dof_info[dof_handler_index].dofs_per_face[active_fe_index];
}



template <int dim, typename Number>
inline unsigned int
MatrixFree<dim, Number>::get_n_q_points_face(
  const unsigned int quad_index,
  const unsigned int active_fe_index) const
{
  AssertIndexRange(quad_index, mapping_info.face_data.size());
  return mapping_info.face_data[quad_index]
    .descriptor[active_fe_index]
    .n_q_points;
}



template <int dim, typename Number>
inline const IndexSet &
MatrixFree<dim, Number>::get_locally_owned_set(
  const unsigned int dof_handler_index) const
{
  return dof_info[dof_handler_index].vector_partitioner->locally_owned_range();
}



template <int dim, typename Number>
inline const IndexSet &
MatrixFree<dim, Number>::get_ghost_set(
  const unsigned int dof_handler_index) const
{
  return dof_info[dof_handler_index].vector_partitioner->ghost_indices();
}



template <int dim, typename Number>
inline const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &
MatrixFree<dim, Number>::get_shape_info(
  const unsigned int dof_handler_index,
  const unsigned int index_quad,
  const unsigned int index_fe,
  const unsigned int active_fe_index,
  const unsigned int active_quad_index) const
{
  AssertIndexRange(dof_handler_index, dof_info.size());
  const unsigned int ind =
    dof_info[dof_handler_index].global_base_element_offset + index_fe;
  AssertIndexRange(ind, shape_info.size(0));
  AssertIndexRange(index_quad, shape_info.size(1));
  AssertIndexRange(active_fe_index, shape_info.size(2));
  AssertIndexRange(active_quad_index, shape_info.size(3));
  return shape_info(ind, index_quad, active_fe_index, active_quad_index);
}



template <int dim, typename Number>
inline const internal::MatrixFreeFunctions::FaceToCellTopology<
  VectorizedArray<Number>::n_array_elements> &
MatrixFree<dim, Number>::get_face_info(const unsigned int macro_face) const
{
  AssertIndexRange(macro_face, face_info.faces.size());
  return face_info.faces[macro_face];
}



template <int dim, typename Number>
inline const Quadrature<dim> &
MatrixFree<dim, Number>::get_quadrature(
  const unsigned int quad_index,
  const unsigned int active_fe_index) const
{
  AssertIndexRange(quad_index, mapping_info.cell_data.size());
  return mapping_info.cell_data[quad_index]
    .descriptor[active_fe_index]
    .quadrature;
}



template <int dim, typename Number>
inline const Quadrature<dim - 1> &
MatrixFree<dim, Number>::get_face_quadrature(
  const unsigned int quad_index,
  const unsigned int active_fe_index) const
{
  AssertIndexRange(quad_index, mapping_info.face_data.size());
  return mapping_info.face_data[quad_index]
    .descriptor[active_fe_index]
    .quadrature;
}



template <int dim, typename Number>
inline unsigned int
MatrixFree<dim, Number>::get_cell_category(const unsigned int macro_cell) const
{
  AssertIndexRange(0, dof_info.size());
  AssertIndexRange(macro_cell, dof_info[0].cell_active_fe_index.size());
  if (dof_info[0].cell_active_fe_index.empty())
    return 0;
  else
    return dof_info[0].cell_active_fe_index[macro_cell];
}



template <int dim, typename Number>
inline std::pair<unsigned int, unsigned int>
MatrixFree<dim, Number>::get_face_category(const unsigned int macro_face) const
{
  AssertIndexRange(macro_face, face_info.faces.size());
  if (dof_info[0].cell_active_fe_index.empty())
    return std::make_pair(0U, 0U);

  std::pair<unsigned int, unsigned int> result;
  for (unsigned int v = 0; v < VectorizedArray<Number>::n_array_elements &&
                           face_info.faces[macro_face].cells_interior[v] !=
                             numbers::invalid_unsigned_int;
       ++v)
    result.first = std::max(
      result.first,
      dof_info[0]
        .cell_active_fe_index[face_info.faces[macro_face].cells_interior[v]]);
  if (face_info.faces[macro_face].cells_exterior[0] !=
      numbers::invalid_unsigned_int)
    for (unsigned int v = 0; v < VectorizedArray<Number>::n_array_elements &&
                             face_info.faces[macro_face].cells_exterior[v] !=
                               numbers::invalid_unsigned_int;
         ++v)
      result.second = std::max(
        result.first,
        dof_info[0]
          .cell_active_fe_index[face_info.faces[macro_face].cells_exterior[v]]);
  else
    result.second = numbers::invalid_unsigned_int;
  return result;
}



template <int dim, typename Number>
inline bool
MatrixFree<dim, Number>::indices_initialized() const
{
  return indices_are_initialized;
}



template <int dim, typename Number>
inline bool
MatrixFree<dim, Number>::mapping_initialized() const
{
  return mapping_is_initialized;
}



template <int dim, typename Number>
AlignedVector<VectorizedArray<Number>> *
MatrixFree<dim, Number>::acquire_scratch_data() const
{
  using list_type =
    std::list<std::pair<bool, AlignedVector<VectorizedArray<Number>>>>;
  list_type &data = scratch_pad.get();
  for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
    if (it->first == false)
      {
        it->first = true;
        return &it->second;
      }
  data.push_front(
    std::make_pair(true, AlignedVector<VectorizedArray<Number>>()));
  return &data.front().second;
}



template <int dim, typename Number>
void
MatrixFree<dim, Number>::release_scratch_data(
  const AlignedVector<VectorizedArray<Number>> *scratch) const
{
  using list_type =
    std::list<std::pair<bool, AlignedVector<VectorizedArray<Number>>>>;
  list_type &data = scratch_pad.get();
  for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
    if (&it->second == scratch)
      {
        Assert(it->first == true, ExcInternalError());
        it->first = false;
        return;
      }
  AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
}



template <int dim, typename Number>
AlignedVector<Number> *
MatrixFree<dim, Number>::acquire_scratch_data_non_threadsafe() const
{
  for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
         scratch_pad_non_threadsafe.begin();
       it != scratch_pad_non_threadsafe.end();
       ++it)
    if (it->first == false)
      {
        it->first = true;
        return &it->second;
      }
  scratch_pad_non_threadsafe.push_front(
    std::make_pair(true, AlignedVector<Number>()));
  return &scratch_pad_non_threadsafe.front().second;
}



template <int dim, typename Number>
void
MatrixFree<dim, Number>::release_scratch_data_non_threadsafe(
  const AlignedVector<Number> *scratch) const
{
  for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
         scratch_pad_non_threadsafe.begin();
       it != scratch_pad_non_threadsafe.end();
       ++it)
    if (&it->second == scratch)
      {
        Assert(it->first == true, ExcInternalError());
        it->first = false;
        return;
      }
  AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
}



// ------------------------------ reinit functions ---------------------------

namespace internal
{
  namespace MatrixFreeImplementation
  {
    template <typename DoFHandlerType>
    inline std::vector<IndexSet>
    extract_locally_owned_index_sets(
      const std::vector<const DoFHandlerType *> &dofh,
      const unsigned int                         level)
    {
      std::vector<IndexSet> locally_owned_set;
      locally_owned_set.reserve(dofh.size());
      for (unsigned int j = 0; j < dofh.size(); j++)
        if (level == numbers::invalid_unsigned_int)
          locally_owned_set.push_back(dofh[j]->locally_owned_dofs());
        else
          AssertThrow(false, ExcNotImplemented());
      return locally_owned_set;
    }

    template <int dim, int spacedim>
    inline std::vector<IndexSet>
    extract_locally_owned_index_sets(
      const std::vector<const ::DoFHandler<dim, spacedim> *> &dofh,
      const unsigned int                                              level)
    {
      std::vector<IndexSet> locally_owned_set;
      locally_owned_set.reserve(dofh.size());
      for (unsigned int j = 0; j < dofh.size(); j++)
        if (level == numbers::invalid_unsigned_int)
          locally_owned_set.push_back(dofh[j]->locally_owned_dofs());
        else
          locally_owned_set.push_back(dofh[j]->locally_owned_mg_dofs(level));
      return locally_owned_set;
    }
  } // namespace MatrixFreeImplementation
} // namespace internal



template <int dim, typename Number>
template <typename DoFHandlerType, typename QuadratureType, typename number2>
void
MatrixFree<dim, Number>::reinit(
  const DoFHandlerType &                                 dof_handler,
  const AffineConstraints<number2> &                     constraints_in,
  const QuadratureType &                                 quad,
  const typename MatrixFree<dim, Number>::AdditionalData additional_data)
{
  std::vector<const DoFHandlerType *>             dof_handlers;
  std::vector<const AffineConstraints<number2> *> constraints;
  std::vector<QuadratureType>                     quads;

  dof_handlers.push_back(&dof_handler);
  constraints.push_back(&constraints_in);
  quads.push_back(quad);

  std::vector<IndexSet> locally_owned_sets =
    internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
      dof_handlers, additional_data.level_mg_handler);

  std::vector<hp::QCollection<1>> quad_hp;
  quad_hp.emplace_back(quad);

  internal_reinit(StaticMappingQ1<dim>::mapping,
                  dof_handlers,
                  constraints,
                  locally_owned_sets,
                  quad_hp,
                  additional_data);
}



template <int dim, typename Number>
template <typename DoFHandlerType, typename QuadratureType, typename number2>
void
MatrixFree<dim, Number>::reinit(
  const Mapping<dim> &                                   mapping,
  const DoFHandlerType &                                 dof_handler,
  const AffineConstraints<number2> &                     constraints_in,
  const QuadratureType &                                 quad,
  const typename MatrixFree<dim, Number>::AdditionalData additional_data)
{
  std::vector<const DoFHandlerType *>             dof_handlers;
  std::vector<const AffineConstraints<number2> *> constraints;

  dof_handlers.push_back(&dof_handler);
  constraints.push_back(&constraints_in);

  std::vector<IndexSet> locally_owned_sets =
    internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
      dof_handlers, additional_data.level_mg_handler);

  std::vector<hp::QCollection<1>> quad_hp;
  quad_hp.emplace_back(quad);

  internal_reinit(mapping,
                  dof_handlers,
                  constraints,
                  locally_owned_sets,
                  quad_hp,
                  additional_data);
}



template <int dim, typename Number>
template <typename DoFHandlerType, typename QuadratureType, typename number2>
void
MatrixFree<dim, Number>::reinit(
  const std::vector<const DoFHandlerType *> &            dof_handler,
  const std::vector<const AffineConstraints<number2> *> &constraint,
  const std::vector<QuadratureType> &                    quad,
  const typename MatrixFree<dim, Number>::AdditionalData additional_data)
{
  std::vector<IndexSet> locally_owned_set =
    internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
      dof_handler, additional_data.level_mg_handler);
  std::vector<hp::QCollection<1>> quad_hp;
  for (unsigned int q = 0; q < quad.size(); ++q)
    quad_hp.emplace_back(quad[q]);
  internal_reinit(StaticMappingQ1<dim>::mapping,
                  dof_handler,
                  constraint,
                  locally_owned_set,
                  quad_hp,
                  additional_data);
}



template <int dim, typename Number>
template <typename DoFHandlerType, typename QuadratureType, typename number2>
void
MatrixFree<dim, Number>::reinit(
  const std::vector<const DoFHandlerType *> &            dof_handler,
  const std::vector<const AffineConstraints<number2> *> &constraint,
  const QuadratureType &                                 quad,
  const typename MatrixFree<dim, Number>::AdditionalData additional_data)
{
  std::vector<IndexSet> locally_owned_set =
    internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
      dof_handler, additional_data.level_mg_handler);
  std::vector<hp::QCollection<1>> quad_hp;
  quad_hp.emplace_back(quad);
  internal_reinit(StaticMappingQ1<dim>::mapping,
                  dof_handler,
                  constraint,
                  locally_owned_set,
                  quad_hp,
                  additional_data);
}



template <int dim, typename Number>
template <typename DoFHandlerType, typename QuadratureType, typename number2>
void
MatrixFree<dim, Number>::reinit(
  const Mapping<dim> &                                   mapping,
  const std::vector<const DoFHandlerType *> &            dof_handler,
  const std::vector<const AffineConstraints<number2> *> &constraint,
  const QuadratureType &                                 quad,
  const typename MatrixFree<dim, Number>::AdditionalData additional_data)
{
  std::vector<IndexSet> locally_owned_set =
    internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
      dof_handler, additional_data.level_mg_handler);
  std::vector<hp::QCollection<1>> quad_hp;
  quad_hp.emplace_back(quad);
  internal_reinit(mapping,
                  dof_handler,
                  constraint,
                  locally_owned_set,
                  quad_hp,
                  additional_data);
}



template <int dim, typename Number>
template <typename DoFHandlerType, typename QuadratureType, typename number2>
void
MatrixFree<dim, Number>::reinit(
  const Mapping<dim> &                                   mapping,
  const std::vector<const DoFHandlerType *> &            dof_handler,
  const std::vector<const AffineConstraints<number2> *> &constraint,
  const std::vector<QuadratureType> &                    quad,
  const typename MatrixFree<dim, Number>::AdditionalData additional_data)
{
  std::vector<IndexSet> locally_owned_set =
    internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
      dof_handler, additional_data.level_mg_handler);
  std::vector<hp::QCollection<1>> quad_hp;
  for (unsigned int q = 0; q < quad.size(); ++q)
    quad_hp.emplace_back(quad[q]);
  internal_reinit(mapping,
                  dof_handler,
                  constraint,
                  locally_owned_set,
                  quad_hp,
                  additional_data);
}



template <int dim, typename Number>
template <typename DoFHandlerType, typename QuadratureType, typename number2>
void
MatrixFree<dim, Number>::reinit(
  const Mapping<dim> &                                   mapping,
  const std::vector<const DoFHandlerType *> &            dof_handler,
  const std::vector<const AffineConstraints<number2> *> &constraint,
  const std::vector<IndexSet> &                          locally_owned_set,
  const std::vector<QuadratureType> &                    quad,
  const typename MatrixFree<dim, Number>::AdditionalData additional_data)
{
  // find out whether we use a hp Quadrature or a standard quadrature
  std::vector<hp::QCollection<1>> quad_hp;
  for (unsigned int q = 0; q < quad.size(); ++q)
    quad_hp.emplace_back(quad[q]);
  internal_reinit(mapping,
                  dof_handler,
                  constraint,
                  locally_owned_set,
                  quad_hp,
                  additional_data);
}



// ------------------------------ implementation of loops --------------------

// internal helper functions that define how to call MPI data exchange
// functions: for generic vectors, do nothing at all. For distributed vectors,
// call update_ghost_values_start function and so on. If we have collections
// of vectors, just do the individual functions of the components. In order to
// keep ghost values consistent (whether we are in read or write mode), we
// also reset the values at the end. the whole situation is a bit complicated
// by the fact that we need to treat block vectors differently, which use some
// additional helper functions to select the blocks and template magic.
namespace internal
{
  template <int dim, typename Number>
  struct VectorDataExchange
  {
    // An arbitrary shift for communication to reduce the risk for accidental
    // interaction with other open communications that a user program might
    // set up
    static constexpr unsigned int channel_shift = 103;

    VectorDataExchange(
      const ::MatrixFree<dim, Number> &matrix_free,
      const typename ::MatrixFree<dim, Number>::DataAccessOnFaces
                         vector_face_access,
      const unsigned int n_components)
      : matrix_free(matrix_free)
      , vector_face_access(
          matrix_free.get_task_info().face_partition_data.empty() ?
            ::MatrixFree<dim, Number>::DataAccessOnFaces::unspecified :
            vector_face_access)
      , ghosts_were_set(false)
#  ifdef DEAL_II_WITH_MPI
      , tmp_data(n_components)
      , requests(n_components)
#  endif
    {
      (void)n_components;
      if (this->vector_face_access !=
          ::MatrixFree<dim, Number>::DataAccessOnFaces::unspecified)
        for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
          AssertDimension(
            matrix_free.get_dof_info(c).vector_partitioner_face_variants.size(),
            3);
    }

    ~VectorDataExchange()
    {
#  ifdef DEAL_II_WITH_MPI
      for (unsigned int i = 0; i < tmp_data.size(); ++i)
        if (tmp_data[i] != nullptr)
          matrix_free.release_scratch_data_non_threadsafe(tmp_data[i]);
#  endif
    }

    unsigned int
    find_vector_in_mf(const LinearAlgebra::distributed::Vector<Number> &vec,
                      const bool check_global_compatibility = true) const
    {
      unsigned int mf_component = numbers::invalid_unsigned_int;
      (void)check_global_compatibility;
      for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
        if (
#  ifdef DEBUG
          check_global_compatibility ?
            vec.get_partitioner()->is_globally_compatible(
              *matrix_free.get_dof_info(c).vector_partitioner) :
#  endif
            vec.get_partitioner()->is_compatible(
              *matrix_free.get_dof_info(c).vector_partitioner))
          {
            mf_component = c;
            break;
          }
      return mf_component;
    }

    const Utilities::MPI::Partitioner &
    get_partitioner(const unsigned int mf_component) const
    {
      AssertDimension(matrix_free.get_dof_info(mf_component)
                        .vector_partitioner_face_variants.size(),
                      3);
      if (vector_face_access ==
          ::MatrixFree<dim, Number>::DataAccessOnFaces::none)
        return *matrix_free.get_dof_info(mf_component)
                  .vector_partitioner_face_variants[0];
      else if (vector_face_access ==
               ::MatrixFree<dim, Number>::DataAccessOnFaces::values)
        return *matrix_free.get_dof_info(mf_component)
                  .vector_partitioner_face_variants[1];
      else
        return *matrix_free.get_dof_info(mf_component)
                  .vector_partitioner_face_variants[2];
    }

    void
    update_ghost_values_start(
      const unsigned int component_in_block_vector,
      const LinearAlgebra::distributed::Vector<Number> &vec)
    {
      (void)component_in_block_vector;
      bool ghosts_set = vec.has_ghost_elements();
      if (ghosts_set)
        ghosts_were_set = true;
      if (vector_face_access ==
            ::MatrixFree<dim, Number>::DataAccessOnFaces::unspecified ||
          vec.size() == 0)
        vec.update_ghost_values_start(component_in_block_vector +
                                      channel_shift);
      else
        {
#  ifdef DEAL_II_WITH_MPI
          const unsigned int mf_component = find_vector_in_mf(vec);
          if (&get_partitioner(mf_component) ==
              matrix_free.get_dof_info(mf_component).vector_partitioner.get())
            {
              vec.update_ghost_values_start(component_in_block_vector +
                                            channel_shift);
              return;
            }

          const Utilities::MPI::Partitioner &part =
            get_partitioner(mf_component);
          if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0)
            return;

          tmp_data[component_in_block_vector] =
            matrix_free.acquire_scratch_data_non_threadsafe();
          tmp_data[component_in_block_vector]->resize_fast(
            part.n_import_indices());
          AssertDimension(requests.size(), tmp_data.size());

          part.export_to_ghosted_array_start(
            component_in_block_vector + channel_shift,
            ArrayView<const Number>(vec.begin(), part.local_size()),
            ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
                              part.n_import_indices()),
            ArrayView<Number>(const_cast<Number *>(vec.begin()) +
                                vec.get_partitioner()->local_size(),
                              vec.get_partitioner()->n_ghost_indices()),
            this->requests[component_in_block_vector]);
#  endif
        }
    }

    void
    update_ghost_values_finish(
      const unsigned int component_in_block_vector,
      const LinearAlgebra::distributed::Vector<Number> &vec)
    {
      (void)component_in_block_vector;
      if (vector_face_access ==
            ::MatrixFree<dim, Number>::DataAccessOnFaces::unspecified ||
          vec.size() == 0)
        vec.update_ghost_values_finish();
      else
        {
#  ifdef DEAL_II_WITH_MPI

          AssertIndexRange(component_in_block_vector, tmp_data.size());
          AssertDimension(requests.size(), tmp_data.size());

          const unsigned int mf_component = find_vector_in_mf(vec);
          const Utilities::MPI::Partitioner &part =
            get_partitioner(mf_component);
          if (&part ==
              matrix_free.get_dof_info(mf_component).vector_partitioner.get())
            {
              vec.update_ghost_values_finish();
              return;
            }

          if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0)
            return;

          part.export_to_ghosted_array_finish(
            ArrayView<Number>(const_cast<Number *>(vec.begin()) +
                                vec.get_partitioner()->local_size(),
                              vec.get_partitioner()->n_ghost_indices()),
            this->requests[component_in_block_vector]);

          matrix_free.release_scratch_data_non_threadsafe(
            tmp_data[component_in_block_vector]);
          tmp_data[component_in_block_vector] = nullptr;
#  endif
        }
    }

    void
    compress_start(const unsigned int component_in_block_vector,
                   LinearAlgebra::distributed::Vector<Number> &vec)
    {
      (void)component_in_block_vector;
      Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
      if (vector_face_access ==
            ::MatrixFree<dim, Number>::DataAccessOnFaces::unspecified ||
          vec.size() == 0)
        vec.compress_start(component_in_block_vector + channel_shift);
      else
        {
#  ifdef DEAL_II_WITH_MPI

          const unsigned int mf_component = find_vector_in_mf(vec);
          const Utilities::MPI::Partitioner &part =
            get_partitioner(mf_component);
          if (&part ==
              matrix_free.get_dof_info(mf_component).vector_partitioner.get())
            {
              vec.compress_start(component_in_block_vector + channel_shift);
              return;
            }

          if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0)
            return;

          tmp_data[component_in_block_vector] =
            matrix_free.acquire_scratch_data_non_threadsafe();
          tmp_data[component_in_block_vector]->resize_fast(
            part.n_import_indices());
          AssertDimension(requests.size(), tmp_data.size());

          part.import_from_ghosted_array_start(
            ::VectorOperation::add,
            component_in_block_vector + channel_shift,
            ArrayView<Number>(vec.begin() + vec.get_partitioner()->local_size(),
                              vec.get_partitioner()->n_ghost_indices()),
            ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
                              part.n_import_indices()),
            this->requests[component_in_block_vector]);
#  endif
        }
    }

    void
    compress_finish(const unsigned int component_in_block_vector,
                    LinearAlgebra::distributed::Vector<Number> &vec)
    {
      (void)component_in_block_vector;
      if (vector_face_access ==
            ::MatrixFree<dim, Number>::DataAccessOnFaces::unspecified ||
          vec.size() == 0)
        vec.compress_finish(::VectorOperation::add);
      else
        {
#  ifdef DEAL_II_WITH_MPI
          AssertIndexRange(component_in_block_vector, tmp_data.size());
          AssertDimension(requests.size(), tmp_data.size());

          const unsigned int mf_component = find_vector_in_mf(vec);

          const Utilities::MPI::Partitioner &part =
            get_partitioner(mf_component);
          if (&part ==
              matrix_free.get_dof_info(mf_component).vector_partitioner.get())
            {
              vec.compress_finish(::VectorOperation::add);
              return;
            }

          if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0)
            return;

          part.import_from_ghosted_array_finish(
            VectorOperation::add,
            ArrayView<const Number>(
              tmp_data[component_in_block_vector]->begin(),
              part.n_import_indices()),
            ArrayView<Number>(vec.begin(), part.local_size()),
            ArrayView<Number>(vec.begin() + vec.get_partitioner()->local_size(),
                              vec.get_partitioner()->n_ghost_indices()),
            this->requests[component_in_block_vector]);

          matrix_free.release_scratch_data_non_threadsafe(
            tmp_data[component_in_block_vector]);
          tmp_data[component_in_block_vector] = nullptr;
#  endif
        }
    }

    void
    reset_ghost_values(
      const LinearAlgebra::distributed::Vector<Number> &vec) const
    {
      if (ghosts_were_set == true)
        return;

      if (vector_face_access ==
            ::MatrixFree<dim, Number>::DataAccessOnFaces::unspecified ||
          vec.size() == 0)
        vec.zero_out_ghosts();
      else
        {
#  ifdef DEAL_II_WITH_MPI
          AssertDimension(requests.size(), tmp_data.size());

          const unsigned int mf_component = find_vector_in_mf(vec);
          const Utilities::MPI::Partitioner &part =
            get_partitioner(mf_component);
          if (&part ==
              matrix_free.get_dof_info(mf_component).vector_partitioner.get())
            vec.zero_out_ghosts();
          else if (part.n_ghost_indices() > 0)
            {
              for (std::vector<std::pair<unsigned int, unsigned int>>::
                     const_iterator my_ghosts =
                       part.ghost_indices_within_larger_ghost_set().begin();
                   my_ghosts !=
                   part.ghost_indices_within_larger_ghost_set().end();
                   ++my_ghosts)
                for (unsigned int j = my_ghosts->first; j < my_ghosts->second;
                     j++)
                  {
                    const_cast<LinearAlgebra::distributed::Vector<Number> &>(
                      vec)
                      .local_element(j + part.local_size()) = 0.;
                  }
            }
#  endif
        }
    }

    void
    zero_vector_region(const unsigned int                          range_index,
                       LinearAlgebra::distributed::Vector<Number> &vec) const
    {
      if (range_index == numbers::invalid_unsigned_int)
        vec = Number();
      else
        {
          const unsigned int mf_component = find_vector_in_mf(vec, false);
          const internal::MatrixFreeFunctions::DoFInfo &dof_info =
            matrix_free.get_dof_info(mf_component);
          Assert(dof_info.vector_zero_range_list_index.empty() == false,
                 ExcNotInitialized());

          Assert(vec.partitioners_are_compatible(*dof_info.vector_partitioner),
                 ExcInternalError());
          AssertIndexRange(range_index,
                           dof_info.vector_zero_range_list_index.size() - 1);
          for (unsigned int id =
                 dof_info.vector_zero_range_list_index[range_index];
               id != dof_info.vector_zero_range_list_index[range_index + 1];
               ++id)
            {
              const unsigned int start_pos =
                dof_info.vector_zero_range_list[id] *
                internal::MatrixFreeFunctions::DoFInfo::chunk_size_zero_vector;
              const unsigned int end_pos =
                std::min((dof_info.vector_zero_range_list[id] + 1) *
                           internal::MatrixFreeFunctions::DoFInfo::
                             chunk_size_zero_vector,
                         dof_info.vector_partitioner->local_size() +
                           dof_info.vector_partitioner->n_ghost_indices());
              std::memset(vec.begin() + start_pos,
                          0,
                          (end_pos - start_pos) * sizeof(Number));
            }
        }
    }

    const ::MatrixFree<dim, Number> &matrix_free;
    const typename ::MatrixFree<dim, Number>::DataAccessOnFaces
         vector_face_access;
    bool ghosts_were_set;
#  ifdef DEAL_II_WITH_MPI
    std::vector<AlignedVector<Number> *>  tmp_data;
    std::vector<std::vector<MPI_Request>> requests;
#  endif
  };

  template <typename VectorStruct>
  unsigned int
  n_components(const VectorStruct &vec);

  template <typename VectorStruct>
  unsigned int
  n_components_block(const VectorStruct &vec,
                     std::integral_constant<bool, true>)
  {
    unsigned int components = 0;
    for (unsigned int bl = 0; bl < vec.n_blocks(); ++bl)
      components += n_components(vec.block(bl));
    return components;
  }

  template <typename VectorStruct>
  unsigned int
  n_components_block(const VectorStruct &, std::integral_constant<bool, false>)
  {
    return 1;
  }

  template <typename VectorStruct>
  unsigned int
  n_components(const VectorStruct &vec)
  {
    return n_components_block(
      vec, std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
  }

  template <typename VectorStruct>
  inline unsigned int
  n_components(const std::vector<VectorStruct> &vec)
  {
    unsigned int components = 0;
    for (unsigned int comp = 0; comp < vec.size(); comp++)
      components += n_components_block(
        vec[comp],
        std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
    return components;
  }

  template <typename VectorStruct>
  inline unsigned int
  n_components(const std::vector<VectorStruct *> &vec)
  {
    unsigned int components = 0;
    for (unsigned int comp = 0; comp < vec.size(); comp++)
      components += n_components_block(
        *vec[comp],
        std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
    return components;
  }

  template <int dim, typename VectorStruct, typename Number>
  void
  update_ghost_values_start_block(const VectorStruct &vec,
                                  const unsigned int  channel,
                                  std::integral_constant<bool, true>,
                                  VectorDataExchange<dim, Number> &exchanger);
  template <int dim, typename VectorStruct, typename Number>
  void
  reset_ghost_values_block(const VectorStruct &vec,
                           std::integral_constant<bool, true>,
                           VectorDataExchange<dim, Number> &exchanger);
  template <int dim, typename VectorStruct, typename Number>
  void
  update_ghost_values_finish_block(const VectorStruct &vec,
                                   const unsigned int  channel,
                                   std::integral_constant<bool, true>,
                                   VectorDataExchange<dim, Number> &exchanger);
  template <int dim, typename VectorStruct, typename Number>
  void
  compress_start_block(const VectorStruct &vec,
                       const unsigned int  channel,
                       std::integral_constant<bool, true>,
                       VectorDataExchange<dim, Number> &exchanger);
  template <int dim, typename VectorStruct, typename Number>
  void
  compress_finish_block(const VectorStruct &vec,
                        const unsigned int  channel,
                        std::integral_constant<bool, true>,
                        VectorDataExchange<dim, Number> &exchanger);
  template <int dim, typename VectorStruct, typename Number>
  void
  zero_vector_region_block(const unsigned int range_index,
                           VectorStruct &,
                           std::integral_constant<bool, true>,
                           VectorDataExchange<dim, Number> &);

  template <int dim, typename VectorStruct, typename Number>
  void
  update_ghost_values_start_block(const VectorStruct &,
                                  const unsigned int,
                                  std::integral_constant<bool, false>,
                                  VectorDataExchange<dim, Number> &)
  {}
  template <int dim, typename VectorStruct, typename Number>
  void
  reset_ghost_values_block(const VectorStruct &,
                           std::integral_constant<bool, false>,
                           VectorDataExchange<dim, Number> &)
  {}
  template <int dim, typename VectorStruct, typename Number>
  void
  update_ghost_values_finish_block(const VectorStruct &,
                                   const unsigned int,
                                   std::integral_constant<bool, false>,
                                   VectorDataExchange<dim, Number> &)
  {}
  template <int dim, typename VectorStruct, typename Number>
  void
  compress_start_block(const VectorStruct &,
                       const unsigned int,
                       std::integral_constant<bool, false>,
                       VectorDataExchange<dim, Number> &)
  {}
  template <int dim, typename VectorStruct, typename Number>
  void
  compress_finish_block(const VectorStruct &,
                        const unsigned int,
                        std::integral_constant<bool, false>,
                        VectorDataExchange<dim, Number> &)
  {}
  template <int dim, typename VectorStruct, typename Number>
  void
  zero_vector_region_block(const unsigned int range_index,
                           VectorStruct &     vec,
                           std::integral_constant<bool, false>,
                           VectorDataExchange<dim, Number> &)
  {
    if (range_index == 0 || range_index == numbers::invalid_unsigned_int)
      vec = 0;
  }



  // if the input vector did not have ghosts imported, clear them here again
  // in order to avoid subsequent operations e.g. in linear solvers to work
  // with ghosts all the time
  template <int dim, typename VectorStruct, typename Number>
  inline void
  reset_ghost_values(const VectorStruct &             vec,
                     VectorDataExchange<dim, Number> &exchanger)
  {
    reset_ghost_values_block(
      vec,
      std::integral_constant<bool, IsBlockVector<VectorStruct>::value>(),
      exchanger);
  }



  template <int dim, typename Number, typename Number2>
  inline void
  reset_ghost_values(const LinearAlgebra::distributed::Vector<Number> &vec,
                     VectorDataExchange<dim, Number2> &exchanger)
  {
    exchanger.reset_ghost_values(vec);
  }



  template <int dim, typename VectorStruct, typename Number>
  inline void
  reset_ghost_values(const std::vector<VectorStruct> &vec,
                     VectorDataExchange<dim, Number> &exchanger)
  {
    // return immediately if there is nothing to do.
    if (exchanger.ghosts_were_set == true)
      return;

    for (unsigned int comp = 0; comp < vec.size(); comp++)
      reset_ghost_values(vec[comp], exchanger);
  }



  template <int dim, typename VectorStruct, typename Number>
  inline void
  reset_ghost_values(const std::vector<VectorStruct *> &vec,
                     VectorDataExchange<dim, Number> &  exchanger)
  {
    // return immediately if there is nothing to do.
    if (exchanger.ghosts_were_set == true)
      return;

    for (unsigned int comp = 0; comp < vec.size(); comp++)
      reset_ghost_values(*vec[comp], exchanger);
  }



  template <int dim, typename VectorStruct, typename Number>
  inline void
  reset_ghost_values_block(const VectorStruct &vec,
                           std::integral_constant<bool, true>,
                           VectorDataExchange<dim, Number> &exchanger)
  {
    // return immediately if there is nothing to do.
    if (exchanger.ghosts_were_set == true)
      return;

    for (unsigned int i = 0; i < vec.n_blocks(); ++i)
      reset_ghost_values(vec.block(i), exchanger);
  }



  // A helper function to identify block vectors with many components where we
  // should not try to overlap computations and communcation because there
  // would be too many outstanding communcation requests. This is the base case
  // for generic vectors
  template <typename VectorStruct>
  constexpr unsigned int
  get_communication_block_size(const VectorStruct &)
  {
    return numbers::invalid_unsigned_int;
  }



  // Specialized case for the block vector which as the additional member
  // variable
  template <typename Number>
  constexpr unsigned int
  get_communication_block_size(
    const LinearAlgebra::distributed::BlockVector<Number> &)
  {
    return LinearAlgebra::distributed::BlockVector<
      Number>::communication_block_size;
  }



  template <int dim, typename VectorStruct, typename Number>
  inline void
  update_ghost_values_start(const VectorStruct &             vec,
                            VectorDataExchange<dim, Number> &exchanger,
                            const unsigned int               channel = 0)
  {
    update_ghost_values_start_block(
      vec,
      channel,
      std::integral_constant<bool, IsBlockVector<VectorStruct>::value>(),
      exchanger);
  }



  template <int dim, typename Number, typename Number2>
  inline void
  update_ghost_values_start(
    const LinearAlgebra::distributed::Vector<Number> &vec,
    VectorDataExchange<dim, Number2> &                exchanger,
    const unsigned int                                channel = 0)
  {
    exchanger.update_ghost_values_start(channel, vec);
  }



  template <int dim, typename VectorStruct, typename Number>
  inline void
  update_ghost_values_start(const std::vector<VectorStruct> &vec,
                            VectorDataExchange<dim, Number> &exchanger)
  {
    unsigned int component_index = 0;
    for (unsigned int comp = 0; comp < vec.size(); comp++)
      {
        update_ghost_values_start(vec[comp], exchanger, component_index);
        component_index += n_components(vec[comp]);
      }
  }



  template <int dim, typename VectorStruct, typename Number>
  inline void
  update_ghost_values_start(const std::vector<VectorStruct *> &vec,
                            VectorDataExchange<dim, Number> &  exchanger)
  {
    unsigned int component_index = 0;
    for (unsigned int comp = 0; comp < vec.size(); comp++)
      {
        update_ghost_values_start(*vec[comp], exchanger, component_index);
        component_index += n_components(*vec[comp]);
      }
  }



  template <int dim, typename VectorStruct, typename Number>
  inline void
  update_ghost_values_start_block(const VectorStruct &vec,
                                  const unsigned int  channel,
                                  std::integral_constant<bool, true>,
                                  VectorDataExchange<dim, Number> &exchanger)
  {
    if (get_communication_block_size(vec) < vec.n_blocks())
      {
        // don't forget to set ghosts_were_set, that otherwise happens
        // inside VectorDataExchange::update_ghost_values_start()
        exchanger.ghosts_were_set = vec.has_ghost_elements();
        vec.update_ghost_values();
      }
    else
      {
        for (unsigned int i = 0; i < vec.n_blocks(); ++i)
          update_ghost_values_start(vec.block(i), exchanger, channel + i);
      }
  }



  template <int dim, typename VectorStruct, typename Number>
  inline void
  update_ghost_values_finish(const VectorStruct &             vec,
                             VectorDataExchange<dim, Number> &exchanger,
                             const unsigned int               channel = 0)
  {
    update_ghost_values_finish_block(
      vec,
      channel,
      std::integral_constant<bool, IsBlockVector<VectorStruct>::value>(),
      exchanger);
  }



  template <int dim, typename Number, typename Number2>
  inline void
  update_ghost_values_finish(
    const LinearAlgebra::distributed::Vector<Number> &vec,
    VectorDataExchange<dim, Number2> &                exchanger,
    const unsigned int                                channel = 0)
  {
    exchanger.update_ghost_values_finish(channel, vec);
  }



  template <int dim, typename VectorStruct, typename Number>
  inline void
  update_ghost_values_finish(const std::vector<VectorStruct> &vec,
                             VectorDataExchange<dim, Number> &exchanger)
  {
    unsigned int component_index = 0;
    for (unsigned int comp = 0; comp < vec.size(); comp++)
      {
        update_ghost_values_finish(vec[comp], exchanger, component_index);
        component_index += n_components(vec[comp]);
      }
  }



  template <int dim, typename VectorStruct, typename Number>
  inline void
  update_ghost_values_finish(const std::vector<VectorStruct *> &vec,
                             VectorDataExchange<dim, Number> &  exchanger)
  {
    unsigned int component_index = 0;
    for (unsigned int comp = 0; comp < vec.size(); comp++)
      {
        update_ghost_values_finish(*vec[comp], exchanger, component_index);
        component_index += n_components(*vec[comp]);
      }
  }



  template <int dim, typename VectorStruct, typename Number>
  inline void
  update_ghost_values_finish_block(const VectorStruct &vec,
                                   const unsigned int  channel,
                                   std::integral_constant<bool, true>,
                                   VectorDataExchange<dim, Number> &exchanger)
  {
    if (get_communication_block_size(vec) < vec.n_blocks())
      {
        // do nothing, everything has already been completed in the _start()
        // call
      }
    else
      for (unsigned int i = 0; i < vec.n_blocks(); ++i)
        update_ghost_values_finish(vec.block(i), exchanger, channel + i);
  }



  template <int dim, typename VectorStruct, typename Number>
  inline void
  compress_start(VectorStruct &                   vec,
                 VectorDataExchange<dim, Number> &exchanger,
                 const unsigned int               channel = 0)
  {
    compress_start_block(
      vec,
      channel,
      std::integral_constant<bool, IsBlockVector<VectorStruct>::value>(),
      exchanger);
  }



  template <int dim, typename Number, typename Number2>
  inline void
  compress_start(LinearAlgebra::distributed::Vector<Number> &vec,
                 VectorDataExchange<dim, Number2> &          exchanger,
                 const unsigned int                          channel = 0)
  {
    exchanger.compress_start(channel, vec);
  }



  template <int dim, typename VectorStruct, typename Number>
  inline void
  compress_start(std::vector<VectorStruct> &      vec,
                 VectorDataExchange<dim, Number> &exchanger)
  {
    unsigned int component_index = 0;
    for (unsigned int comp = 0; comp < vec.size(); comp++)
      {
        compress_start(vec[comp], exchanger, component_index);
        component_index += n_components(vec[comp]);
      }
  }



  template <int dim, typename VectorStruct, typename Number>
  inline void
  compress_start(std::vector<VectorStruct *> &    vec,
                 VectorDataExchange<dim, Number> &exchanger)
  {
    unsigned int component_index = 0;
    for (unsigned int comp = 0; comp < vec.size(); comp++)
      {
        compress_start(*vec[comp], exchanger, component_index);
        component_index += n_components(*vec[comp]);
      }
  }



  template <int dim, typename VectorStruct, typename Number>
  inline void
  compress_start_block(VectorStruct &     vec,
                       const unsigned int channel,
                       std::integral_constant<bool, true>,
                       VectorDataExchange<dim, Number> &exchanger)
  {
    if (get_communication_block_size(vec) < vec.n_blocks())
      vec.compress(::VectorOperation::add);
    else
      for (unsigned int i = 0; i < vec.n_blocks(); ++i)
        compress_start(vec.block(i), exchanger, channel + i);
  }



  template <int dim, typename VectorStruct, typename Number>
  inline void
  compress_finish(VectorStruct &                   vec,
                  VectorDataExchange<dim, Number> &exchanger,
                  const unsigned int               channel = 0)
  {
    compress_finish_block(
      vec,
      channel,
      std::integral_constant<bool, IsBlockVector<VectorStruct>::value>(),
      exchanger);
  }



  template <int dim, typename Number, typename Number2>
  inline void
  compress_finish(LinearAlgebra::distributed::Vector<Number> &vec,
                  VectorDataExchange<dim, Number2> &          exchanger,
                  const unsigned int                          channel = 0)
  {
    exchanger.compress_finish(channel, vec);
  }



  template <int dim, typename VectorStruct, typename Number>
  inline void
  compress_finish(std::vector<VectorStruct> &      vec,
                  VectorDataExchange<dim, Number> &exchanger)
  {
    unsigned int component_index = 0;
    for (unsigned int comp = 0; comp < vec.size(); comp++)
      {
        compress_finish(vec[comp], exchanger, component_index);
        component_index += n_components(vec[comp]);
      }
  }



  template <int dim, typename VectorStruct, typename Number>
  inline void
  compress_finish(std::vector<VectorStruct *> &    vec,
                  VectorDataExchange<dim, Number> &exchanger)
  {
    unsigned int component_index = 0;
    for (unsigned int comp = 0; comp < vec.size(); comp++)
      {
        compress_finish(*vec[comp], exchanger, component_index);
        component_index += n_components(*vec[comp]);
      }
  }



  template <int dim, typename VectorStruct, typename Number>
  inline void
  compress_finish_block(VectorStruct &     vec,
                        const unsigned int channel,
                        std::integral_constant<bool, true>,
                        VectorDataExchange<dim, Number> &exchanger)
  {
    if (get_communication_block_size(vec) < vec.n_blocks())
      {
        // do nothing, everything has already been completed in the _start()
        // call
      }
    else
      for (unsigned int i = 0; i < vec.n_blocks(); ++i)
        compress_finish(vec.block(i), exchanger, channel + i);
  }



  template <int dim, typename VectorStruct, typename Number>
  inline void
  zero_vector_region(const unsigned int               range_index,
                     VectorStruct &                   vec,
                     VectorDataExchange<dim, Number> &exchanger)
  {
    zero_vector_region_block(
      range_index,
      vec,
      std::integral_constant<bool, IsBlockVector<VectorStruct>::value>(),
      exchanger);
  }



  template <int dim, typename Number, typename Number2>
  inline void
  zero_vector_region(const unsigned int                          range_index,
                     LinearAlgebra::distributed::Vector<Number> &vec,
                     VectorDataExchange<dim, Number2> &          exchanger)
  {
    exchanger.zero_vector_region(range_index, vec);
  }



  template <int dim, typename VectorStruct, typename Number>
  inline void
  zero_vector_region(const unsigned int               range_index,
                     std::vector<VectorStruct> &      vec,
                     VectorDataExchange<dim, Number> &exchanger)
  {
    for (unsigned int comp = 0; comp < vec.size(); comp++)
      zero_vector_region(range_index, vec[comp], exchanger);
  }



  template <int dim, typename VectorStruct, typename Number>
  inline void
  zero_vector_region(const unsigned int               range_index,
                     std::vector<VectorStruct *> &    vec,
                     VectorDataExchange<dim, Number> &exchanger)
  {
    for (unsigned int comp = 0; comp < vec.size(); comp++)
      zero_vector_region(range_index, *vec[comp], exchanger);
  }



  template <int dim, typename VectorStruct, typename Number>
  inline void
  zero_vector_region_block(const unsigned int range_index,
                           VectorStruct &     vec,
                           std::integral_constant<bool, true>,
                           VectorDataExchange<dim, Number> &exchanger)
  {
    for (unsigned int i = 0; i < vec.n_blocks(); ++i)
      zero_vector_region(range_index, vec.block(i), exchanger);
  }



  namespace MatrixFreeFunctions
  {
    // struct to select between a const interface and a non-const interface
    // for MFWorker
    template <typename, typename, typename, typename, bool>
    struct InterfaceSelector
    {};

    // Version for constant functions
    template <typename MF,
              typename InVector,
              typename OutVector,
              typename Container>
    struct InterfaceSelector<MF, InVector, OutVector, Container, true>
    {
      using function_type = void (Container::*)(
        const MF &,
        OutVector &,
        const InVector &,
        const std::pair<unsigned int, unsigned int> &) const;
    };

    // Version for non-constant functions
    template <typename MF,
              typename InVector,
              typename OutVector,
              typename Container>
    struct InterfaceSelector<MF, InVector, OutVector, Container, false>
    {
      using function_type =
        void (Container::*)(const MF &,
                            OutVector &,
                            const InVector &,
                            const std::pair<unsigned int, unsigned int> &);
    };
  } // namespace MatrixFreeFunctions



  // A implementation class for the worker object that runs the various
  // operations we want to perform during the matrix-free loop
  template <typename MF,
            typename InVector,
            typename OutVector,
            typename Container,
            bool is_constant>
  class MFWorker : public MFWorkerInterface
  {
  public:
    // An alias to make the arguments further down more readable
    using function_type = typename MatrixFreeFunctions::
      InterfaceSelector<MF, InVector, OutVector, Container, is_constant>::
        function_type;

    // constructor, binds all the arguments to this class
    MFWorker(const MF &                           matrix_free,
             const InVector &                     src,
             OutVector &                          dst,
             const bool                           zero_dst_vector_setting,
             const Container &                    container,
             function_type                        cell_function,
             function_type                        face_function,
             function_type                        boundary_function,
             const typename MF::DataAccessOnFaces src_vector_face_access =
               MF::DataAccessOnFaces::none,
             const typename MF::DataAccessOnFaces dst_vector_face_access =
               MF::DataAccessOnFaces::none)
      : matrix_free(matrix_free)
      , container(const_cast<Container &>(container))
      , cell_function(cell_function)
      , face_function(face_function)
      , boundary_function(boundary_function)
      , src(src)
      , dst(dst)
      , src_data_exchanger(matrix_free,
                           src_vector_face_access,
                           n_components(src))
      , dst_data_exchanger(matrix_free,
                           dst_vector_face_access,
                           n_components(dst))
      , src_and_dst_are_same(PointerComparison::equal(&src, &dst))
      , zero_dst_vector_setting(zero_dst_vector_setting &&
                                !src_and_dst_are_same)
    {}

    // Runs the cell work. If no function is given, nothing is done
    virtual void
    cell(const std::pair<unsigned int, unsigned int> &cell_range) override
    {
      if (cell_function != nullptr && cell_range.second > cell_range.first)
        (container.*
         cell_function)(matrix_free, this->dst, this->src, cell_range);
    }

    // Runs the assembler on interior faces. If no function is given, nothing
    // is done
    virtual void
    face(const std::pair<unsigned int, unsigned int> &face_range) override
    {
      if (face_function != nullptr && face_range.second > face_range.first)
        (container.*
         face_function)(matrix_free, this->dst, this->src, face_range);
    }

    // Runs the assembler on boundary faces. If no function is given, nothing
    // is done
    virtual void
    boundary(const std::pair<unsigned int, unsigned int> &face_range) override
    {
      if (boundary_function != nullptr && face_range.second > face_range.first)
        (container.*
         boundary_function)(matrix_free, this->dst, this->src, face_range);
    }

    // Starts the communication for the update ghost values operation. We
    // cannot call this update if ghost and destination are the same because
    // that would introduce spurious entries in the destination (there is also
    // the problem that reading from a vector that we also write to is usually
    // not intended in case there is overlap, but this is up to the
    // application code to decide and we cannot catch this case here).
    virtual void
    vector_update_ghosts_start() override
    {
      if (!src_and_dst_are_same)
        internal::update_ghost_values_start(src, src_data_exchanger);
    }

    // Finishes the communication for the update ghost values operation
    virtual void
    vector_update_ghosts_finish() override
    {
      if (!src_and_dst_are_same)
        internal::update_ghost_values_finish(src, src_data_exchanger);
    }

    // Starts the communication for the vector compress operation
    virtual void
    vector_compress_start() override
    {
      internal::compress_start(dst, dst_data_exchanger);
    }

    // Finishes the communication for the vector compress operation
    virtual void
    vector_compress_finish() override
    {
      internal::compress_finish(dst, dst_data_exchanger);
      if (!src_and_dst_are_same)
        internal::reset_ghost_values(src, src_data_exchanger);
    }

    // Zeros the given input vector
    virtual void
    zero_dst_vector_range(const unsigned int range_index) override
    {
      if (zero_dst_vector_setting)
        internal::zero_vector_region(range_index, dst, dst_data_exchanger);
    }

  private:
    const MF &    matrix_free;
    Container &   container;
    function_type cell_function;
    function_type face_function;
    function_type boundary_function;

    const InVector &src;
    OutVector &     dst;
    VectorDataExchange<MF::dimension, typename MF::value_type>
      src_data_exchanger;
    VectorDataExchange<MF::dimension, typename MF::value_type>
               dst_data_exchanger;
    const bool src_and_dst_are_same;
    const bool zero_dst_vector_setting;
  };



  template <class MF, typename InVector, typename OutVector>
  struct MFClassWrapper
  {
    using function_type =
      std::function<void(const MF &,
                         OutVector &,
                         const InVector &,
                         const std::pair<unsigned int, unsigned int> &)>;

    MFClassWrapper(const function_type cell,
                   const function_type face,
                   const function_type boundary)
      : cell(cell)
      , face(face)
      , boundary(boundary)
    {}

    void
    cell_integrator(const MF &                                   mf,
                    OutVector &                                  dst,
                    const InVector &                             src,
                    const std::pair<unsigned int, unsigned int> &range) const
    {
      if (cell)
        cell(mf, dst, src, range);
    }

    void
    face_integrator(const MF &                                   mf,
                    OutVector &                                  dst,
                    const InVector &                             src,
                    const std::pair<unsigned int, unsigned int> &range) const
    {
      if (face)
        face(mf, dst, src, range);
    }

    void
    boundary_integrator(
      const MF &                                   mf,
      OutVector &                                  dst,
      const InVector &                             src,
      const std::pair<unsigned int, unsigned int> &range) const
    {
      if (boundary)
        boundary(mf, dst, src, range);
    }

    const function_type cell;
    const function_type face;
    const function_type boundary;
  };

} // end of namespace internal



template <int dim, typename Number>
template <typename OutVector, typename InVector>
inline void
MatrixFree<dim, Number>::cell_loop(
  const std::function<void(const MatrixFree<dim, Number> &,
                           OutVector &,
                           const InVector &,
                           const std::pair<unsigned int, unsigned int> &)>
    &             cell_operation,
  OutVector &     dst,
  const InVector &src,
  const bool      zero_dst_vector) const
{
  using Wrapper =
    internal::MFClassWrapper<MatrixFree<dim, Number>, InVector, OutVector>;
  Wrapper wrap(cell_operation, nullptr, nullptr);
  internal::
    MFWorker<MatrixFree<dim, Number>, InVector, OutVector, Wrapper, true>
      worker(*this,
             src,
             dst,
             zero_dst_vector,
             wrap,
             &Wrapper::cell_integrator,
             &Wrapper::face_integrator,
             &Wrapper::boundary_integrator);

  task_info.loop(worker);
}



template <int dim, typename Number>
template <typename OutVector, typename InVector>
inline void
MatrixFree<dim, Number>::loop(
  const std::function<void(const MatrixFree<dim, Number> &,
                           OutVector &,
                           const InVector &,
                           const std::pair<unsigned int, unsigned int> &)>
    &cell_operation,
  const std::function<void(const MatrixFree<dim, Number> &,
                           OutVector &,
                           const InVector &,
                           const std::pair<unsigned int, unsigned int> &)>
    &face_operation,
  const std::function<void(const MatrixFree<dim, Number> &,
                           OutVector &,
                           const InVector &,
                           const std::pair<unsigned int, unsigned int> &)>
    &                     boundary_operation,
  OutVector &             dst,
  const InVector &        src,
  const bool              zero_dst_vector,
  const DataAccessOnFaces dst_vector_face_access,
  const DataAccessOnFaces src_vector_face_access) const
{
  using Wrapper =
    internal::MFClassWrapper<MatrixFree<dim, Number>, InVector, OutVector>;
  Wrapper wrap(cell_operation, face_operation, boundary_operation);
  internal::
    MFWorker<MatrixFree<dim, Number>, InVector, OutVector, Wrapper, true>
      worker(*this,
             src,
             dst,
             zero_dst_vector,
             wrap,
             &Wrapper::cell_integrator,
             &Wrapper::face_integrator,
             &Wrapper::boundary_integrator,
             src_vector_face_access,
             dst_vector_face_access);

  task_info.loop(worker);
}



template <int dim, typename Number>
template <typename CLASS, typename OutVector, typename InVector>
inline void
MatrixFree<dim, Number>::cell_loop(
  void (CLASS::*function_pointer)(const MatrixFree<dim, Number> &,
                                  OutVector &,
                                  const InVector &,
                                  const std::pair<unsigned int, unsigned int> &)
    const,
  const CLASS *   owning_class,
  OutVector &     dst,
  const InVector &src,
  const bool      zero_dst_vector) const
{
  internal::MFWorker<MatrixFree<dim, Number>, InVector, OutVector, CLASS, true>
    worker(*this,
           src,
           dst,
           zero_dst_vector,
           *owning_class,
           function_pointer,
           nullptr,
           nullptr);
  task_info.loop(worker);
}



template <int dim, typename Number>
template <typename CLASS, typename OutVector, typename InVector>
inline void
MatrixFree<dim, Number>::loop(
  void (CLASS::*cell_operation)(const MatrixFree<dim, Number> &,
                                OutVector &,
                                const InVector &,
                                const std::pair<unsigned int, unsigned int> &)
    const,
  void (CLASS::*face_operation)(const MatrixFree<dim, Number> &,
                                OutVector &,
                                const InVector &,
                                const std::pair<unsigned int, unsigned int> &)
    const,
  void (CLASS::*boundary_operation)(
    const MatrixFree<dim, Number> &,
    OutVector &,
    const InVector &,
    const std::pair<unsigned int, unsigned int> &) const,
  const CLASS *           owning_class,
  OutVector &             dst,
  const InVector &        src,
  const bool              zero_dst_vector,
  const DataAccessOnFaces dst_vector_face_access,
  const DataAccessOnFaces src_vector_face_access) const
{
  internal::MFWorker<MatrixFree<dim, Number>, InVector, OutVector, CLASS, true>
    worker(*this,
           src,
           dst,
           zero_dst_vector,
           *owning_class,
           cell_operation,
           face_operation,
           boundary_operation,
           src_vector_face_access,
           dst_vector_face_access);
  task_info.loop(worker);
}



template <int dim, typename Number>
template <typename CLASS, typename OutVector, typename InVector>
inline void
MatrixFree<dim, Number>::cell_loop(
  void (CLASS::*function_pointer)(
    const MatrixFree<dim, Number> &,
    OutVector &,
    const InVector &,
    const std::pair<unsigned int, unsigned int> &),
  CLASS *         owning_class,
  OutVector &     dst,
  const InVector &src,
  const bool      zero_dst_vector) const
{
  internal::MFWorker<MatrixFree<dim, Number>, InVector, OutVector, CLASS, false>
    worker(*this,
           src,
           dst,
           zero_dst_vector,
           *owning_class,
           function_pointer,
           nullptr,
           nullptr);
  task_info.loop(worker);
}



template <int dim, typename Number>
template <typename CLASS, typename OutVector, typename InVector>
inline void
MatrixFree<dim, Number>::loop(
  void (CLASS::*cell_operation)(const MatrixFree<dim, Number> &,
                                OutVector &,
                                const InVector &,
                                const std::pair<unsigned int, unsigned int> &),
  void (CLASS::*face_operation)(const MatrixFree<dim, Number> &,
                                OutVector &,
                                const InVector &,
                                const std::pair<unsigned int, unsigned int> &),
  void (CLASS::*boundary_operation)(
    const MatrixFree<dim, Number> &,
    OutVector &,
    const InVector &,
    const std::pair<unsigned int, unsigned int> &),
  CLASS *                 owning_class,
  OutVector &             dst,
  const InVector &        src,
  const bool              zero_dst_vector,
  const DataAccessOnFaces dst_vector_face_access,
  const DataAccessOnFaces src_vector_face_access) const
{
  internal::MFWorker<MatrixFree<dim, Number>, InVector, OutVector, CLASS, false>
    worker(*this,
           src,
           dst,
           zero_dst_vector,
           *owning_class,
           cell_operation,
           face_operation,
           boundary_operation,
           src_vector_face_access,
           dst_vector_face_access);
  task_info.loop(worker);
}


#endif // ifndef DOXYGEN



DEAL_II_NAMESPACE_CLOSE

#endif