matrix_free.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2011 - 2017 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 at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------


#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/parallel.h>
#include <deal.II/base/quadrature.h>
#include <deal.II/base/vectorization.h>
#include <deal.II/base/thread_local_storage.h>
#include <deal.II/base/template_constraints.h>
#include <deal.II/fe/fe.h>
#include <deal.II/fe/mapping.h>
#include <deal.II/fe/mapping_q1.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/la_parallel_vector.h>
#include <deal.II/lac/block_vector_base.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/hp/dof_handler.h>
#include <deal.II/hp/q_collection.h>
#include <deal.II/matrix_free/helper_functions.h>
#include <deal.II/matrix_free/shape_info.h>
#include <deal.II/matrix_free/dof_info.h>
#include <deal.II/matrix_free/mapping_info.h>

#ifdef DEAL_II_WITH_THREADS
#include <tbb/task.h>
#include <tbb/task_scheduler_init.h>
#include <tbb/parallel_for.h>
#include <tbb/blocked_range.h>
#endif

#include <stdlib.h>
#include <memory>
#include <limits>
#include <list>


DEAL_II_NAMESPACE_OPEN



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

  struct AdditionalData
  {
    enum TasksParallelScheme
    {
      none,
      partition_partition,
      partition_color,
      color
    };

    // avoid warning about use of deprecated variables
    DEAL_II_DISABLE_EXTRA_DIAGNOSTICS

    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 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)
      :
      mpi_communicator      (MPI_COMM_SELF),
      tasks_parallel_scheme (tasks_parallel_scheme),
      tasks_block_size      (tasks_block_size),
      mapping_update_flags  (mapping_update_flags),
      level_mg_handler      (level_mg_handler),
      store_plain_indices   (store_plain_indices),
      initialize_indices    (initialize_indices),
      initialize_mapping    (initialize_mapping)
    {};

    AdditionalData (const MPI_Comm            mpi_communicator,
                    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 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)
      :
      mpi_communicator      (mpi_communicator),
      tasks_parallel_scheme (tasks_parallel_scheme),
      tasks_block_size      (tasks_block_size),
      mapping_update_flags  (mapping_update_flags),
      level_mg_handler      (level_mg_handler),
      store_plain_indices   (store_plain_indices),
      initialize_indices    (initialize_indices),
      initialize_mapping    (initialize_mapping)
    {} DEAL_II_DEPRECATED

    MPI_Comm            mpi_communicator DEAL_II_DEPRECATED;

    TasksParallelScheme tasks_parallel_scheme;

    unsigned int        tasks_block_size;

    UpdateFlags         mapping_update_flags;

    unsigned int        level_mg_handler;

    bool                store_plain_indices;

    bool                initialize_indices;

    bool                initialize_mapping;
  };

  DEAL_II_ENABLE_EXTRA_DIAGNOSTICS

  MatrixFree ();

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

  ~MatrixFree();

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

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

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

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

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

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

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

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

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

  void clear();


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

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

  template <typename CLASS, typename OutVector, typename InVector>
  void cell_loop (void (CLASS::*function_pointer)(const MatrixFree &,
                                                  OutVector &,
                                                  const InVector &,
                                                  const std::pair<unsigned int,
                                                  unsigned int> &),
                  CLASS          *owning_class,
                  OutVector      &dst,
                  const InVector &src) 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 vector_component = 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 vector_component = 0) const;


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

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

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

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

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

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

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


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

  unsigned int n_components () const;

  unsigned int n_physical_cells () const;

  unsigned int n_macro_cells () const;

  const DoFHandler<dim> &
  get_dof_handler (const unsigned int fe_component = 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 fe_component = 0) const;

  bool
  at_irregular_cell (const unsigned int macro_cell_number) const;

  unsigned int
  n_components_filled (const unsigned int macro_cell_number) const;

  unsigned int
  get_dofs_per_cell (const unsigned int fe_component = 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;

  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;

  const internal::MatrixFreeFunctions::SizeInfo &
  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 fe_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<Number> &
  get_shape_info (const unsigned int fe_component = 0,
                  const unsigned int quad_index   = 0,
                  const unsigned int hp_active_fe_index = 0,
                  const unsigned int hp_active_quad_index = 0) const;

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

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


private:

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

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

  void
  initialize_indices (const std::vector<const ConstraintMatrix *> &constraint,
                      const std::vector<IndexSet> &locally_owned_set);

  void initialize_dof_handlers (const std::vector<const DoFHandler<dim>*> &dof_handlers,
                                const unsigned int                         level);

  void initialize_dof_handlers (const std::vector<const hp::DoFHandler<dim>*> &dof_handlers,
                                const unsigned int                             level);

  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<Number> > shape_info;

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

  internal::MatrixFreeFunctions::SizeInfo size_info;

  internal::MatrixFreeFunctions::TaskInfo task_info;

  bool                               indices_are_initialized;

  bool                               mapping_is_initialized;

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



/*----------------------- 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_cxx11::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
const internal::MatrixFreeFunctions::TaskInfo &
MatrixFree<dim,Number>::get_task_info () const
{
  return task_info;
}



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



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



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



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 (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() ? 0 :
         &constraint_pool_data[0] + 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() ? 0 :
         &constraint_pool_data[0] + 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 vector_component) const
{
  AssertIndexRange (vector_component, dof_info.size());
  if (dof_info[vector_component].cell_active_fe_index.empty())
    {
      AssertDimension (dof_info[vector_component].fe_index_conversion.size(),1);
      if (dof_info[vector_component].fe_index_conversion[0].first == degree)
        return range;
      else
        return std::pair<unsigned int,unsigned int> (range.second,range.second);
    }

  const unsigned int fe_index =
    dof_info[vector_component].fe_index_from_degree(degree);
  if (fe_index >= dof_info[vector_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, vector_component);
}



template <int dim, typename Number>
inline
std::pair<unsigned int,unsigned int>
MatrixFree<dim,Number>::create_cell_subrange_hp_by_index
(const std::pair<unsigned int,unsigned int> &range,
 const unsigned int fe_index,
 const unsigned int vector_component) const
{
  AssertIndexRange (fe_index, dof_info[vector_component].max_fe_index);
  const std::vector<unsigned int> &fe_indices =
    dof_info[vector_component].cell_active_fe_index;
  if (fe_indices.size() == 0)
    return range;
  else
    {
      // the range over which we are searching must be ordered, otherwise we
      // got a range that spans over too many cells
#ifdef DEBUG
      for (unsigned int i=range.first+1; i<range.second; ++i)
        Assert (fe_indices[i] >= fe_indices[i-1],
                ExcMessage ("Cell range must be over sorted range of fe indices in hp case!"));
      AssertIndexRange(range.first,fe_indices.size()+1);
      AssertIndexRange(range.second,fe_indices.size()+1);
#endif
      std::pair<unsigned int,unsigned int> return_range;
      return_range.first =
        std::lower_bound(&fe_indices[0] + range.first,
                         &fe_indices[0] + range.second, fe_index)
        -&fe_indices[0] ;
      return_range.second =
        std::lower_bound(&fe_indices[0] + return_range.first,
                         &fe_indices[0] + range.second,
                         fe_index + 1)-&fe_indices[0];
      Assert(return_range.first >= range.first &&
             return_range.second <= range.second, ExcInternalError());
      return return_range;
    }
}



template <int dim, typename Number>
inline
void
MatrixFree<dim,Number>::renumber_dofs (std::vector<types::global_dof_index> &renumbering,
                                       const unsigned int vector_component)
{
  AssertIndexRange(vector_component, dof_info.size());
  dof_info[vector_component].renumber_dofs (renumbering);
}



template <int dim, typename Number>
inline
const DoFHandler<dim> &
MatrixFree<dim,Number>::get_dof_handler (const unsigned int dof_index) const
{
  AssertIndexRange (dof_index, n_components());
  if (dof_handlers.active_dof_handler == DoFHandlers::usual)
    {
      AssertDimension (dof_handlers.dof_handler.size(),
                       dof_handlers.n_dof_handlers);
      return *dof_handlers.dof_handler[dof_index];
    }
  else
    {
      Assert (false, ExcNotImplemented());
      // put pseudo return argument to avoid compiler error, but trigger a
      // segfault in case this is only run in optimized mode
      return *dof_handlers.dof_handler[numbers::invalid_unsigned_int];
    }
}



template <int dim, typename Number>
inline
typename DoFHandler<dim>::cell_iterator
MatrixFree<dim,Number>::get_cell_iterator(const unsigned int macro_cell_number,
                                          const unsigned int vector_number,
                                          const unsigned int dof_index) const
{
  const unsigned int vectorization_length=VectorizedArray<Number>::n_array_elements;
#ifdef DEBUG
  AssertIndexRange (dof_index, dof_handlers.n_dof_handlers);
  AssertIndexRange (macro_cell_number, size_info.n_macro_cells);
  AssertIndexRange (vector_number, vectorization_length);
  const unsigned int irreg_filled = dof_info[dof_index].row_starts[macro_cell_number][2];
  if (irreg_filled > 0)
    AssertIndexRange (vector_number, irreg_filled);
#endif

  const DoFHandler<dim> *dofh = 0;
  if (dof_handlers.active_dof_handler == DoFHandlers::usual)
    {
      AssertDimension (dof_handlers.dof_handler.size(),
                       dof_handlers.n_dof_handlers);
      dofh = dof_handlers.dof_handler[dof_index];
    }
  else
    {
      Assert (false, ExcMessage ("Cannot return DoFHandler<dim>::cell_iterator "
                                 "for underlying DoFHandler!"));
    }

  std::pair<unsigned int,unsigned int> index =
    cell_level_index[macro_cell_number*vectorization_length+vector_number];
  return typename DoFHandler<dim>::cell_iterator
         (&dofh->get_triangulation(), index.first, index.second, dofh);
}



template <int dim, typename Number>
inline
typename hp::DoFHandler<dim>::active_cell_iterator
MatrixFree<dim,Number>::get_hp_cell_iterator(const unsigned int macro_cell_number,
                                             const unsigned int vector_number,
                                             const unsigned int dof_index) const
{
  const unsigned int vectorization_length=VectorizedArray<Number>::n_array_elements;
#ifdef DEBUG
  AssertIndexRange (dof_index, dof_handlers.n_dof_handlers);
  AssertIndexRange (macro_cell_number, size_info.n_macro_cells);
  AssertIndexRange (vector_number, vectorization_length);
  const unsigned int irreg_filled = dof_info[dof_index].row_starts[macro_cell_number][2];
  if (irreg_filled > 0)
    AssertIndexRange (vector_number, irreg_filled);
#endif

  Assert (dof_handlers.active_dof_handler == DoFHandlers::hp,
          ExcNotImplemented());
  const hp::DoFHandler<dim> *dofh = dof_handlers.hp_dof_handler[dof_index];
  std::pair<unsigned int,unsigned int> index =
    cell_level_index[macro_cell_number*vectorization_length+vector_number];
  return typename hp::DoFHandler<dim>::cell_iterator
         (&dofh->get_triangulation(), index.first, index.second, dofh);
}



template <int dim, typename Number>
inline
bool
MatrixFree<dim,Number>::at_irregular_cell (const unsigned int macro_cell) const
{
  AssertIndexRange (macro_cell, size_info.n_macro_cells);
  return dof_info[0].row_starts[macro_cell][2] > 0;
}



template <int dim, typename Number>
inline
unsigned int
MatrixFree<dim,Number>::n_components_filled (const unsigned int macro_cell) const
{
  AssertIndexRange (macro_cell, size_info.n_macro_cells);
  const unsigned int n_filled = dof_info[0].row_starts[macro_cell][2];
  if (n_filled == 0)
    return VectorizedArray<Number>::n_array_elements;
  else
    return n_filled;
}



template <int dim, typename Number>
inline
unsigned int
MatrixFree<dim,Number>::get_dofs_per_cell(const unsigned int dof_index,
                                          const unsigned int active_fe_index) const
{
  AssertIndexRange (dof_index, dof_info.size());
  return dof_info[dof_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.mapping_data_gen.size());
  return mapping_info.mapping_data_gen[quad_index].n_q_points[active_fe_index];
}



template <int dim, typename Number>
inline
unsigned int
MatrixFree<dim,Number>::get_dofs_per_face(const unsigned int dof_index,
                                          const unsigned int active_fe_index) const
{
  AssertIndexRange (dof_index, dof_info.size());
  return dof_info[dof_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.mapping_data_gen.size());
  return mapping_info.mapping_data_gen[quad_index].n_q_points_face[active_fe_index];
}



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



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



template <int dim, typename Number>
inline
const internal::MatrixFreeFunctions::ShapeInfo<Number> &
MatrixFree<dim,Number>::get_shape_info (const unsigned int index_fe,
                                        const unsigned int index_quad,
                                        const unsigned int active_fe_index,
                                        const unsigned int active_quad_index) const
{
  AssertIndexRange (index_fe, 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(index_fe, index_quad,
                    active_fe_index, active_quad_index);
}



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.mapping_data_gen.size());
  return mapping_info.mapping_data_gen[quad_index].
         quadrature[active_fe_index];
}



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.mapping_data_gen.size());
  return mapping_info.mapping_data_gen[quad_index].
         face_quadrature[active_fe_index];
}



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
{
  typedef std::list<std::pair<bool, AlignedVector<VectorizedArray<Number> > > > list_type;
  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
{
  typedef std::list<std::pair<bool, AlignedVector<VectorizedArray<Number> > > > list_type;
  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"));
}



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

namespace internal
{
  namespace MatrixFree
  {
    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;
    }
  }
}



template <int dim, typename Number>
template <typename DoFHandlerType, typename QuadratureType>
void MatrixFree<dim,Number>::
reinit(const DoFHandlerType                                  &dof_handler,
       const ConstraintMatrix                                &constraints_in,
       const QuadratureType                                  &quad,
       const typename MatrixFree<dim,Number>::AdditionalData additional_data)
{
  std::vector<const DoFHandlerType *>   dof_handlers;
  std::vector<const ConstraintMatrix *> 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::MatrixFree::extract_locally_owned_index_sets
    (dof_handlers, additional_data.level_mg_handler);

  std::vector<hp::QCollection<1> > quad_hp;
  quad_hp.push_back (hp::QCollection<1>(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>
void MatrixFree<dim,Number>::
reinit(const Mapping<dim>                                    &mapping,
       const DoFHandlerType                                  &dof_handler,
       const ConstraintMatrix                                &constraints_in,
       const QuadratureType                                  &quad,
       const typename MatrixFree<dim,Number>::AdditionalData additional_data)
{
  std::vector<const DoFHandlerType *>   dof_handlers;
  std::vector<const ConstraintMatrix *> constraints;

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

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

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

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



template <int dim, typename Number>
template <typename DoFHandlerType, typename QuadratureType>
void MatrixFree<dim,Number>::
reinit(const std::vector<const DoFHandlerType *>   &dof_handler,
       const std::vector<const ConstraintMatrix *> &constraint,
       const std::vector<QuadratureType>           &quad,
       const typename MatrixFree<dim,Number>::AdditionalData additional_data)
{
  std::vector<IndexSet> locally_owned_set =
    internal::MatrixFree::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.push_back (hp::QCollection<1>(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>
void MatrixFree<dim,Number>::
reinit(const std::vector<const DoFHandlerType *>             &dof_handler,
       const std::vector<const ConstraintMatrix *>           &constraint,
       const QuadratureType                                  &quad,
       const typename MatrixFree<dim,Number>::AdditionalData additional_data)
{
  std::vector<IndexSet> locally_owned_set =
    internal::MatrixFree::extract_locally_owned_index_sets
    (dof_handler, additional_data.level_mg_handler);
  std::vector<hp::QCollection<1> > quad_hp;
  quad_hp.push_back (hp::QCollection<1>(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>
void MatrixFree<dim,Number>::
reinit(const Mapping<dim>                                    &mapping,
       const std::vector<const DoFHandlerType *>             &dof_handler,
       const std::vector<const ConstraintMatrix *>           &constraint,
       const QuadratureType                                  &quad,
       const typename MatrixFree<dim,Number>::AdditionalData additional_data)
{
  std::vector<IndexSet> locally_owned_set =
    internal::MatrixFree::extract_locally_owned_index_sets
    (dof_handler, additional_data.level_mg_handler);
  std::vector<hp::QCollection<1> > quad_hp;
  quad_hp.push_back (hp::QCollection<1>(quad));
  internal_reinit(mapping, dof_handler,constraint,
                  locally_owned_set, quad_hp, additional_data);
}



template <int dim, typename Number>
template <typename DoFHandlerType, typename QuadratureType>
void MatrixFree<dim,Number>::
reinit(const Mapping<dim>                                   &mapping,
       const std::vector<const DoFHandlerType *>            &dof_handler,
       const std::vector<const ConstraintMatrix *>          &constraint,
       const std::vector<QuadratureType>                    &quad,
       const typename MatrixFree<dim,Number>::AdditionalData additional_data)
{
  std::vector<IndexSet> locally_owned_set =
    internal::MatrixFree::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.push_back (hp::QCollection<1>(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>
void MatrixFree<dim,Number>::
reinit(const Mapping<dim>                                    &mapping,
       const std::vector<const DoFHandlerType *>             &dof_handler,
       const std::vector<const ConstraintMatrix *>           &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.push_back (hp::QCollection<1>(quad[q]));
  internal_reinit (mapping,
                   dof_handler,
                   constraint, locally_owned_set, quad_hp, additional_data);
}



// ------------------------------ implementation of cell_loop ---------------

// 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). 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<typename VectorStruct>
  bool update_ghost_values_start_block (const VectorStruct &vec,
                                        const unsigned int channel,
                                        internal::bool2type<true>);
  template<typename VectorStruct>
  void reset_ghost_values_block (const VectorStruct &vec,
                                 const bool          zero_out_ghosts,
                                 internal::bool2type<true>);
  template<typename VectorStruct>
  void update_ghost_values_finish_block (const VectorStruct &vec,
                                         internal::bool2type<true>);
  template<typename VectorStruct>
  void compress_start_block (const VectorStruct &vec,
                             const unsigned int channel,
                             internal::bool2type<true>);
  template<typename VectorStruct>
  void compress_finish_block (const VectorStruct &vec,
                              internal::bool2type<true>);

  template<typename VectorStruct>
  bool update_ghost_values_start_block (const VectorStruct &,
                                        const unsigned int,
                                        internal::bool2type<false>)
  {
    return false;
  }
  template<typename VectorStruct>
  void reset_ghost_values_block (const VectorStruct &,
                                 const bool,
                                 internal::bool2type<false>)
  {}
  template<typename VectorStruct>
  void update_ghost_values_finish_block (const VectorStruct &,
                                         internal::bool2type<false>)
  {}
  template<typename VectorStruct>
  void compress_start_block (const VectorStruct &,
                             const unsigned int,
                             internal::bool2type<false>)
  {}
  template<typename VectorStruct>
  void compress_finish_block (const VectorStruct &,
                              internal::bool2type<false>)
  {}



  // returns true if the vector was in a state without ghost values before,
  // i.e., we need to zero out ghosts in the very end
  template<typename VectorStruct>
  inline
  bool update_ghost_values_start (const VectorStruct &vec,
                                  const unsigned int channel = 0)
  {
    return
      update_ghost_values_start_block(vec, channel,
                                      internal::bool2type<IsBlockVector<VectorStruct>::value>());
  }



  template<typename Number>
  inline
  bool update_ghost_values_start (const LinearAlgebra::distributed::Vector<Number> &vec,
                                  const unsigned int                  channel = 0)
  {
    bool return_value = !vec.has_ghost_elements();
    vec.update_ghost_values_start(channel);
    return return_value;
  }



  template <typename VectorStruct>
  inline
  bool update_ghost_values_start (const std::vector<VectorStruct> &vec)
  {
    bool return_value = false;
    for (unsigned int comp=0; comp<vec.size(); comp++)
      return_value = update_ghost_values_start(vec[comp], comp);
    return return_value;
  }



  template <typename VectorStruct>
  inline
  bool update_ghost_values_start (const std::vector<VectorStruct *> &vec)
  {
    bool return_value = false;
    for (unsigned int comp=0; comp<vec.size(); comp++)
      return_value = update_ghost_values_start(*vec[comp], comp);
    return return_value;
  }



  template<typename VectorStruct>
  inline
  bool update_ghost_values_start_block (const VectorStruct &vec,
                                        const unsigned int channel,
                                        internal::bool2type<true>)
  {
    bool return_value = false;
    for (unsigned int i=0; i<vec.n_blocks(); ++i)
      return_value = update_ghost_values_start(vec.block(i), channel+509*i);
    return return_value;
  }



  // 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<typename VectorStruct>
  inline
  void reset_ghost_values (const VectorStruct &vec,
                           const bool          zero_out_ghosts)
  {
    reset_ghost_values_block(vec, zero_out_ghosts,
                             internal::bool2type<IsBlockVector<VectorStruct>::value>());
  }



  template<typename Number>
  inline
  void reset_ghost_values (const LinearAlgebra::distributed::Vector<Number> &vec,
                           const bool zero_out_ghosts)
  {
    if (zero_out_ghosts)
      const_cast<LinearAlgebra::distributed::Vector<Number>&>(vec).zero_out_ghosts();
  }



  template <typename VectorStruct>
  inline
  void reset_ghost_values (const std::vector<VectorStruct> &vec,
                           const bool zero_out_ghosts)
  {
    for (unsigned int comp=0; comp<vec.size(); comp++)
      reset_ghost_values(vec[comp], zero_out_ghosts);
  }



  template <typename VectorStruct>
  inline
  void reset_ghost_values (const std::vector<VectorStruct *> &vec,
                           const bool zero_out_ghosts)
  {
    for (unsigned int comp=0; comp<vec.size(); comp++)
      reset_ghost_values(*vec[comp], zero_out_ghosts);
  }



  template<typename VectorStruct>
  inline
  void reset_ghost_values_block (const VectorStruct &vec,
                                 const bool          zero_out_ghosts,
                                 internal::bool2type<true>)
  {
    for (unsigned int i=0; i<vec.n_blocks(); ++i)
      reset_ghost_values(vec.block(i), zero_out_ghosts);
  }



  template <typename VectorStruct>
  inline
  void update_ghost_values_finish (const VectorStruct &vec)
  {
    update_ghost_values_finish_block(vec,
                                     internal::bool2type<IsBlockVector<VectorStruct>::value>());
  }



  template <typename Number>
  inline
  void update_ghost_values_finish (const LinearAlgebra::distributed::Vector<Number> &vec)
  {
    vec.update_ghost_values_finish();
  }



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



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



  template <typename VectorStruct>
  inline
  void update_ghost_values_finish_block (const VectorStruct &vec,
                                         internal::bool2type<true>)
  {
    for (unsigned int i=0; i<vec.n_blocks(); ++i)
      update_ghost_values_finish(vec.block(i));
  }



  template <typename VectorStruct>
  inline
  void compress_start (VectorStruct &vec,
                       const unsigned int channel = 0)
  {
    compress_start_block (vec, channel,
                          internal::bool2type<IsBlockVector<VectorStruct>::value>());
  }



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



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



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



  template <typename VectorStruct>
  inline
  void compress_start_block (VectorStruct      &vec,
                             const unsigned int channel,
                             internal::bool2type<true>)
  {
    for (unsigned int i=0; i<vec.n_blocks(); ++i)
      compress_start(vec.block(i), channel + 500*i);
  }



  template <typename VectorStruct>
  inline
  void compress_finish (VectorStruct &vec)
  {
    compress_finish_block(vec,
                          internal::bool2type<IsBlockVector<VectorStruct>::value>());
  }



  template <typename Number>
  inline
  void compress_finish (LinearAlgebra::distributed::Vector<Number> &vec)
  {
    vec.compress_finish(::VectorOperation::add);
  }



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



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



  template <typename VectorStruct>
  inline
  void compress_finish_block (VectorStruct &vec,
                              internal::bool2type<true>)
  {
    for (unsigned int i=0; i<vec.n_blocks(); ++i)
      compress_finish(vec.block(i));
  }



#ifdef DEAL_II_WITH_THREADS

  // This defines the TBB data structures that are needed to schedule the
  // partition-partition variant

  namespace partition
  {
    template<typename Worker>
    class CellWork : public tbb::task
    {
    public:
      CellWork (const Worker &worker_in,
                const unsigned int partition_in,
                const internal::MatrixFreeFunctions::TaskInfo &task_info_in,
                const bool is_blocked_in)
        :
        dummy (NULL),
        worker (worker_in),
        partition (partition_in),
        task_info (task_info_in),
        is_blocked (is_blocked_in)
      {};
      tbb::task *execute ()
      {
        std::pair<unsigned int, unsigned int> cell_range
        (task_info.partition_color_blocks_data[partition],
         task_info.partition_color_blocks_data[partition+1]);
        worker(cell_range);
        if (is_blocked==true)
          dummy->spawn (*dummy);
        return NULL;
      }

      tbb::empty_task *dummy;

    private:
      const Worker      &worker;
      const unsigned int partition;
      const internal::MatrixFreeFunctions::TaskInfo &task_info;
      const bool         is_blocked;
    };



    template<typename Worker>
    class PartitionWork : public tbb::task
    {
    public:
      PartitionWork (const Worker &function_in,
                     const unsigned int partition_in,
                     const internal::MatrixFreeFunctions::TaskInfo &task_info_in,
                     const bool    is_blocked_in = false)
        :
        dummy (NULL),
        function (function_in),
        partition (partition_in),
        task_info (task_info_in),
        is_blocked (is_blocked_in)
      {};
      tbb::task *execute ()
      {
        tbb::empty_task *root = new( tbb::task::allocate_root() )
        tbb::empty_task;
        unsigned int evens = task_info.partition_evens[partition];
        unsigned int odds  = task_info.partition_odds[partition];
        unsigned int n_blocked_workers =
          task_info.partition_n_blocked_workers[partition];
        unsigned int n_workers = task_info.partition_n_workers[partition];
        std::vector<CellWork<Worker>*> worker(n_workers);
        std::vector<CellWork<Worker>*> blocked_worker(n_blocked_workers);

        root->set_ref_count(evens+1);
        for (unsigned int j=0; j<evens; j++)
          {
            worker[j] = new(root->allocate_child())
            CellWork<Worker>(function, task_info.
                             partition_color_blocks_row_index[partition]+2*j,
                             task_info, false);
            if (j>0)
              {
                worker[j]->set_ref_count(2);
                blocked_worker[j-1]->dummy = new(worker[j]->allocate_child())
                tbb::empty_task;
                worker[j-1]->spawn(*blocked_worker[j-1]);
              }
            else
              worker[j]->set_ref_count(1);
            if (j<evens-1)
              {
                blocked_worker[j] = new(worker[j]->allocate_child())
                CellWork<Worker>(function, task_info.
                                 partition_color_blocks_row_index
                                 [partition] + 2*j+1, task_info, true);
              }
            else
              {
                if (odds==evens)
                  {
                    worker[evens] = new(worker[j]->allocate_child())
                    CellWork<Worker>(function, task_info.
                                     partition_color_blocks_row_index[partition]+2*j+1,
                                     task_info, false);
                    worker[j]->spawn(*worker[evens]);
                  }
                else
                  {
                    tbb::empty_task *child = new(worker[j]->allocate_child())
                    tbb::empty_task();
                    worker[j]->spawn(*child);
                  }
              }
          }

        root->wait_for_all();
        root->destroy(*root);
        if (is_blocked==true)
          dummy->spawn (*dummy);
        return NULL;
      }

      tbb::empty_task *dummy;

    private:
      const Worker  &function;
      const unsigned int partition;
      const internal::MatrixFreeFunctions::TaskInfo &task_info;
      const bool     is_blocked;
    };

  } // end of namespace partition



  namespace color
  {
    template <typename Worker>
    class CellWork
    {
    public:
      CellWork (const Worker                   &worker_in,
                const internal::MatrixFreeFunctions::TaskInfo &task_info_in)
        :
        worker (worker_in),
        task_info (task_info_in)
      {};
      void operator()(const tbb::blocked_range<unsigned int> &r) const
      {
        for (unsigned int block=r.begin(); block<r.end(); block++)
          {
            std::pair<unsigned int,unsigned int> cell_range;
            if (task_info.position_short_block<block)
              {
                cell_range.first = (block-1)*task_info.block_size+
                                   task_info.block_size_last;
                cell_range.second = cell_range.first + task_info.block_size;
              }
            else
              {
                cell_range.first = block*task_info.block_size;
                cell_range.second = cell_range.first +
                                    ((block == task_info.position_short_block)?
                                     (task_info.block_size_last):(task_info.block_size));
              }
            worker (cell_range);
          }
      }
    private:
      const Worker   &worker;
      const internal::MatrixFreeFunctions::TaskInfo &task_info;
    };


    template<typename Worker>
    class PartitionWork : public tbb::task
    {
    public:
      PartitionWork (const Worker &worker_in,
                     const unsigned int partition_in,
                     const internal::MatrixFreeFunctions::TaskInfo &task_info_in,
                     const bool    is_blocked_in)
        :
        dummy (NULL),
        worker (worker_in),
        partition (partition_in),
        task_info (task_info_in),
        is_blocked (is_blocked_in)
      {};
      tbb::task *execute ()
      {
        unsigned int lower = task_info.partition_color_blocks_data[partition],
                     upper = task_info.partition_color_blocks_data[partition+1];
        parallel_for(tbb::blocked_range<unsigned int>(lower,upper,1),
                     CellWork<Worker> (worker,task_info));
        if (is_blocked==true)
          dummy->spawn (*dummy);
        return NULL;
      }

      tbb::empty_task *dummy;

    private:
      const Worker &worker;
      const unsigned int partition;
      const internal::MatrixFreeFunctions::TaskInfo &task_info;
      const bool is_blocked;
    };

  } // end of namespace color


  template<typename VectorStruct>
  class MPIComDistribute : public tbb::task
  {
  public:
    MPIComDistribute (const VectorStruct  &src_in)
      :
      src(src_in)
    {};

    tbb::task *execute ()
    {
      internal::update_ghost_values_finish(src);
      return 0;
    }

  private:
    const VectorStruct &src;
  };



  template<typename VectorStruct>
  class MPIComCompress : public tbb::task
  {
  public:
    MPIComCompress (VectorStruct        &dst_in)
      :
      dst(dst_in)
    {};

    tbb::task *execute ()
    {
      internal::compress_start(dst);
      return 0;
    }

  private:
    VectorStruct &dst;
  };

#endif // DEAL_II_WITH_THREADS

} // end of namespace internal



template <int dim, typename Number>
template <typename OutVector, typename InVector>
inline
void
MatrixFree<dim, Number>::cell_loop
(const std_cxx11::function<void (const MatrixFree<dim,Number> &,
                                 OutVector &,
                                 const InVector &,
                                 const std::pair<unsigned int,
                                 unsigned int> &)> &cell_operation,
 OutVector       &dst,
 const InVector  &src) const
{
  // in any case, need to start the ghost import at the beginning
  bool ghosts_were_not_set = internal::update_ghost_values_start (src);

#ifdef DEAL_II_WITH_THREADS

  // Use multithreading if so requested and if there is enough work to do in
  // parallel (the code might hang if there are less than two chunks!)
  if (task_info.use_multithreading == true && task_info.n_blocks > 3)
    {
      // to simplify the function calls, bind away all arguments except the
      // cell range
      typedef
      std_cxx11::function<void (const std::pair<unsigned int,unsigned int> &range)>
      Worker;

      const Worker func = std_cxx11::bind (std_cxx11::ref(cell_operation),
                                           std_cxx11::cref(*this),
                                           std_cxx11::ref(dst),
                                           std_cxx11::cref(src),
                                           std_cxx11::_1);

      if (task_info.use_partition_partition == true)
        {
          tbb::empty_task *root = new( tbb::task::allocate_root() )
          tbb::empty_task;
          unsigned int evens = task_info.evens;
          unsigned int odds  = task_info.odds;
          root->set_ref_count(evens+1);
          unsigned int n_blocked_workers = task_info.n_blocked_workers;
          unsigned int n_workers = task_info.n_workers;
          std::vector<internal::partition::PartitionWork<Worker>*>
          worker(n_workers);
          std::vector<internal::partition::PartitionWork<Worker>*>
          blocked_worker(n_blocked_workers);
          internal::MPIComCompress<OutVector> *worker_compr =
            new(root->allocate_child())
          internal::MPIComCompress<OutVector>(dst);
          worker_compr->set_ref_count(1);
          for (unsigned int j=0; j<evens; j++)
            {
              if (j>0)
                {
                  worker[j] = new(root->allocate_child())
                  internal::partition::PartitionWork<Worker>
                  (func,2*j,task_info,false);
                  worker[j]->set_ref_count(2);
                  blocked_worker[j-1]->dummy = new(worker[j]->allocate_child())
                  tbb::empty_task;
                  if (j>1)
                    worker[j-1]->spawn(*blocked_worker[j-1]);
                  else
                    worker_compr->spawn(*blocked_worker[j-1]);
                }
              else
                {
                  worker[j] = new(worker_compr->allocate_child())
                  internal::partition::PartitionWork<Worker>
                  (func,2*j,task_info,false);
                  worker[j]->set_ref_count(2);
                  internal::MPIComDistribute<InVector> *worker_dist =
                    new (worker[j]->allocate_child())
                  internal::MPIComDistribute<InVector>(src);
                  worker_dist->spawn(*worker_dist);
                }
              if (j<evens-1)
                {
                  blocked_worker[j] = new(worker[j]->allocate_child())
                  internal::partition::PartitionWork<Worker>
                  (func,2*j+1,task_info,true);
                }
              else
                {
                  if (odds==evens)
                    {
                      worker[evens] = new(worker[j]->allocate_child())
                      internal::partition::PartitionWork<Worker>
                      (func,2*j+1,task_info,false);
                      worker[j]->spawn(*worker[evens]);
                    }
                  else
                    {
                      tbb::empty_task *child = new(worker[j]->allocate_child())
                      tbb::empty_task();
                      worker[j]->spawn(*child);
                    }
                }
            }

          root->wait_for_all();
          root->destroy(*root);
        }
      else // end of partition-partition, start of partition-color
        {
          unsigned int evens = task_info.evens;
          unsigned int odds  = task_info.odds;

          // check whether there is only one partition. if not, build up the
          // tree of partitions
          if (odds > 0)
            {
              tbb::empty_task *root = new( tbb::task::allocate_root() ) tbb::empty_task;
              root->set_ref_count(evens+1);
              unsigned int n_blocked_workers = odds-(odds+evens+1)%2;
              unsigned int n_workers = task_info.partition_color_blocks_data.size()-1-
                                       n_blocked_workers;
              std::vector<internal::color::PartitionWork<Worker>*> worker(n_workers);
              std::vector<internal::color::PartitionWork<Worker>*> blocked_worker(n_blocked_workers);
              unsigned int worker_index = 0, slice_index = 0;
              unsigned int spawn_index =  0;
              int spawn_index_child = -2;
              internal::MPIComCompress<OutVector> *worker_compr = new(root->allocate_child())
              internal::MPIComCompress<OutVector>(dst);
              worker_compr->set_ref_count(1);
              for (unsigned int part=0;
                   part<task_info.partition_color_blocks_row_index.size()-1; part++)
                {
                  const unsigned int spawn_index_new = worker_index;
                  if (part == 0)
                    worker[worker_index] = new(worker_compr->allocate_child())
                    internal::color::PartitionWork<Worker>(func,slice_index,task_info,false);
                  else
                    worker[worker_index] = new(root->allocate_child())
                    internal::color::PartitionWork<Worker>(func,slice_index,task_info,false);
                  slice_index++;
                  for (; slice_index<task_info.partition_color_blocks_row_index[part+1];
                       slice_index++)
                    {
                      worker[worker_index]->set_ref_count(1);
                      worker_index++;
                      worker[worker_index] = new (worker[worker_index-1]->allocate_child())
                      internal::color::PartitionWork<Worker>(func,slice_index,task_info,false);
                    }
                  worker[worker_index]->set_ref_count(2);
                  if (part>0)
                    {
                      blocked_worker[(part-1)/2]->dummy =
                        new (worker[worker_index]->allocate_child()) tbb::empty_task;
                      worker_index++;
                      if (spawn_index_child == -1)
                        worker[spawn_index]->spawn(*blocked_worker[(part-1)/2]);
                      else
                        {
                          Assert(spawn_index_child>=0, ExcInternalError());
                          worker[spawn_index]->spawn(*worker[spawn_index_child]);
                        }
                      spawn_index = spawn_index_new;
                    }
                  else
                    {
                      internal::MPIComDistribute<InVector> *worker_dist =
                        new (worker[worker_index]->allocate_child())
                      internal::MPIComDistribute<InVector>(src);
                      worker_dist->spawn(*worker_dist);
                      worker_index++;
                    }
                  part += 1;
                  if (part<task_info.partition_color_blocks_row_index.size()-1)
                    {
                      if (part<task_info.partition_color_blocks_row_index.size()-2)
                        {
                          blocked_worker[part/2] = new(worker[worker_index-1]->allocate_child())
                          internal::color::PartitionWork<Worker>(func,slice_index,task_info,true);
                          slice_index++;
                          if (slice_index<
                              task_info.partition_color_blocks_row_index[part+1])
                            {
                              blocked_worker[part/2]->set_ref_count(1);
                              worker[worker_index] = new(blocked_worker[part/2]->allocate_child())
                              internal::color::PartitionWork<Worker>(func,slice_index,task_info,false);
                              slice_index++;
                            }
                          else
                            {
                              spawn_index_child = -1;
                              continue;
                            }
                        }
                      for (; slice_index<task_info.partition_color_blocks_row_index[part+1];
                           slice_index++)
                        {
                          if (slice_index>
                              task_info.partition_color_blocks_row_index[part])
                            {
                              worker[worker_index]->set_ref_count(1);
                              worker_index++;
                            }
                          worker[worker_index] = new (worker[worker_index-1]->allocate_child())
                          internal::color::PartitionWork<Worker>(func,slice_index,task_info,false);
                        }
                      spawn_index_child = worker_index;
                      worker_index++;
                    }
                  else
                    {
                      tbb::empty_task *final = new (worker[worker_index-1]->allocate_child())
                      tbb::empty_task;
                      worker[spawn_index]->spawn(*final);
                      spawn_index_child = worker_index-1;
                    }
                }
              if (evens==odds)
                {
                  Assert(spawn_index_child>=0, ExcInternalError());
                  worker[spawn_index]->spawn(*worker[spawn_index_child]);
                }
              root->wait_for_all();
              root->destroy(*root);
            }
          // case when we only have one partition: this is the usual coloring
          // scheme, and we just schedule a parallel for loop for each color
          else
            {
              Assert(evens==1,ExcInternalError());
              internal::update_ghost_values_finish(src);

              for (unsigned int color=0;
                   color < task_info.partition_color_blocks_row_index[1];
                   ++color)
                {
                  unsigned int lower = task_info.partition_color_blocks_data[color],
                               upper = task_info.partition_color_blocks_data[color+1];
                  parallel_for(tbb::blocked_range<unsigned int>(lower,upper,1),
                               internal::color::CellWork<Worker>
                               (func,task_info));
                }

              internal::compress_start(dst);
            }
        }
    }
  else
#endif
    // serial loop
    {
      std::pair<unsigned int,unsigned int> cell_range;

      // First operate on cells where no ghost data is needed (inner cells)
      {
        cell_range.first = 0;
        cell_range.second = size_info.boundary_cells_start;
        cell_operation (*this, dst, src, cell_range);
      }

      // before starting operations on cells that contain ghost nodes (outer
      // cells), wait for the MPI commands to finish
      internal::update_ghost_values_finish(src);

      // For the outer cells, do the same procedure as for inner cells.
      if (size_info.boundary_cells_end > size_info.boundary_cells_start)
        {
          cell_range.first = size_info.boundary_cells_start;
          cell_range.second = size_info.boundary_cells_end;
          cell_operation (*this, dst, src, cell_range);
        }

      internal::compress_start(dst);

      // Finally operate on cells where no ghost data is needed (inner cells)
      if (size_info.n_macro_cells > size_info.boundary_cells_end)
        {
          cell_range.first = size_info.boundary_cells_end;
          cell_range.second = size_info.n_macro_cells;
          cell_operation (*this, dst, src, cell_range);
        }
    }

  // In every case, we need to finish transfers at the very end
  internal::compress_finish(dst);
  internal::reset_ghost_values(src, ghosts_were_not_set);
}



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
{
  // here, use std_cxx11::bind to hand a function handler with the appropriate
  // argument to the other loop function
  std_cxx11::function<void (const MatrixFree<dim,Number> &,
                            OutVector &,
                            const InVector &,
                            const std::pair<unsigned int,
                            unsigned int> &)>
  function = std_cxx11::bind<void>(function_pointer,
                                   owning_class,
                                   std_cxx11::_1,
                                   std_cxx11::_2,
                                   std_cxx11::_3,
                                   std_cxx11::_4);
  cell_loop (function, dst, src);
}



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
{
  // here, use std_cxx11::bind to hand a function handler with the appropriate
  // argument to the other loop function
  std_cxx11::function<void (const MatrixFree<dim,Number> &,
                            OutVector &,
                            const InVector &,
                            const std::pair<unsigned int,
                            unsigned int> &)>
  function = std_cxx11::bind<void>(function_pointer,
                                   owning_class,
                                   std_cxx11::_1,
                                   std_cxx11::_2,
                                   std_cxx11::_3,
                                   std_cxx11::_4);
  cell_loop (function, dst, src);
}


#endif  // ifndef DOXYGEN



DEAL_II_NAMESPACE_CLOSE

#endif