la_parallel_vector.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 2011 - 2021 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_la_parallel_vector_h
#define dealii_la_parallel_vector_h
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/communication_pattern_base.h>
#include <deal.II/base/memory_space.h>
#include <deal.II/base/memory_space_data.h>
#include <deal.II/base/mpi.h>
#include <deal.II/base/numbers.h>
#include <deal.II/base/partitioner.h>
#include <deal.II/base/thread_management.h>
 
#include <deal.II/lac/vector_operation.h>
#include <deal.II/lac/vector_space_vector.h>
#include <deal.II/lac/vector_type_traits.h>
 
#include <iomanip>
#include <memory>
 
DEAL_II_NAMESPACE_OPEN
 
// Forward declarations
#ifndef DOXYGEN
namespace LinearAlgebra
{
  namespace distributed
  {
    template <typename>
    class BlockVector;
  }
 
  template <typename>
  class ReadWriteVector;
} // namespace LinearAlgebra
 
#  ifdef DEAL_II_WITH_PETSC
namespace PETScWrappers
{
  namespace MPI
  {
    class Vector;
  }
} // namespace PETScWrappers
#  endif
 
#  ifdef DEAL_II_WITH_TRILINOS
namespace TrilinosWrappers
{
  namespace MPI
  {
    class Vector;
  }
} // namespace TrilinosWrappers
#  endif
#endif
 
namespace LinearAlgebra
{
  namespace distributed
  {
    template <typename Number, typename MemorySpace = MemorySpace::Host>
    class Vector : public ::LinearAlgebra::VectorSpaceVector<Number>,
                   public Subscriptor
    {
    public:
      using memory_space    = MemorySpace;
      using value_type      = Number;
      using pointer         = value_type *;
      using const_pointer   = const value_type *;
      using iterator        = value_type *;
      using const_iterator  = const value_type *;
      using reference       = value_type &;
      using const_reference = const value_type &;
      using size_type       = types::global_dof_index;
      using real_type       = typename numbers::NumberTraits<Number>::real_type;
 
      static_assert(
        std::is_same<MemorySpace, ::MemorySpace::Host>::value ||
          std::is_same<MemorySpace, ::MemorySpace::CUDA>::value,
        "MemorySpace should be Host or CUDA");
 
      Vector();
 
      Vector(const Vector<Number, MemorySpace> &in_vector);
 
      Vector(const size_type size);
 
      Vector(const IndexSet &local_range,
             const IndexSet &ghost_indices,
             const MPI_Comm &communicator);
 
      Vector(const IndexSet &local_range, const MPI_Comm &communicator);
 
      Vector(
        const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner);
 
      virtual ~Vector() override;
 
      void
      reinit(const size_type size, const bool omit_zeroing_entries = false);
 
      template <typename Number2>
      void
      reinit(const Vector<Number2, MemorySpace> &in_vector,
             const bool                          omit_zeroing_entries = false);
 
      void
      reinit(const IndexSet &local_range,
             const IndexSet &ghost_indices,
             const MPI_Comm &communicator);
 
      void
      reinit(const IndexSet &local_range, const MPI_Comm &communicator);
 
      void
      reinit(
        const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner,
        const MPI_Comm &comm_sm = MPI_COMM_SELF);
 
      void
      reinit(const types::global_dof_index local_size,
             const types::global_dof_index ghost_size,
             const MPI_Comm &              comm,
             const MPI_Comm &              comm_sm = MPI_COMM_SELF);
 
      void
      swap(Vector<Number, MemorySpace> &v);
 
      Vector<Number, MemorySpace> &
      operator=(const Vector<Number, MemorySpace> &in_vector);
 
      template <typename Number2>
      Vector<Number, MemorySpace> &
      operator=(const Vector<Number2, MemorySpace> &in_vector);
 
 
      virtual void
      compress(::VectorOperation::values operation) override;
 
      void
      update_ghost_values() const;
 
      void
      compress_start(
        const unsigned int                communication_channel = 0,
        ::VectorOperation::values operation = VectorOperation::add);
 
      void
      compress_finish(::VectorOperation::values operation);
 
      void
      update_ghost_values_start(
        const unsigned int communication_channel = 0) const;
 
 
      void
      update_ghost_values_finish() const;
 
      DEAL_II_DEPRECATED_EARLY void
      zero_out_ghosts() const;
 
      void
      zero_out_ghost_values() const;
 
      bool
      has_ghost_elements() const;
 
      template <typename Number2>
      void
      copy_locally_owned_data_from(const Vector<Number2, MemorySpace> &src);
 
      template <typename MemorySpace2>
      void
      import(const Vector<Number, MemorySpace2> &src,
             VectorOperation::values             operation);
 
 
 
      virtual void
      reinit(const VectorSpaceVector<Number> &V,
             const bool omit_zeroing_entries = false) override;
 
      virtual Vector<Number, MemorySpace> &
      operator*=(const Number factor) override;
 
      virtual Vector<Number, MemorySpace> &
      operator/=(const Number factor) override;
 
      virtual Vector<Number, MemorySpace> &
      operator+=(const VectorSpaceVector<Number> &V) override;
 
      virtual Vector<Number, MemorySpace> &
      operator-=(const VectorSpaceVector<Number> &V) override;
 
      virtual void
      import(const LinearAlgebra::ReadWriteVector<Number> &V,
             VectorOperation::values                       operation,
             std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
               communication_pattern = {}) override;
 
      virtual Number
      operator*(const VectorSpaceVector<Number> &V) const override;
 
      virtual void
      add(const Number a) override;
 
      virtual void
      add(const Number a, const VectorSpaceVector<Number> &V) override;
 
      virtual void
      add(const Number                     a,
          const VectorSpaceVector<Number> &V,
          const Number                     b,
          const VectorSpaceVector<Number> &W) override;
 
      virtual void
      add(const std::vector<size_type> &indices,
          const std::vector<Number> &   values);
 
      virtual void
      sadd(const Number                     s,
           const Number                     a,
           const VectorSpaceVector<Number> &V) override;
 
      virtual void
      scale(const VectorSpaceVector<Number> &scaling_factors) override;
 
      virtual void
      equ(const Number a, const VectorSpaceVector<Number> &V) override;
 
      virtual real_type
      l1_norm() const override;
 
      virtual real_type
      l2_norm() const override;
 
      real_type
      norm_sqr() const;
 
      virtual real_type
      linfty_norm() const override;
 
      virtual Number
      add_and_dot(const Number                     a,
                  const VectorSpaceVector<Number> &V,
                  const VectorSpaceVector<Number> &W) override;
 
      virtual size_type
      size() const override;
 
      virtual ::IndexSet
      locally_owned_elements() const override;
 
      virtual void
      print(std::ostream &     out,
            const unsigned int precision  = 3,
            const bool         scientific = true,
            const bool         across     = true) const override;
 
      virtual std::size_t
      memory_consumption() const override;
 
 
      virtual Vector<Number, MemorySpace> &
      operator=(const Number s) override;
 
      template <typename OtherNumber>
      void
      add(const std::vector<size_type> &       indices,
          const ::Vector<OtherNumber> &values);
 
      template <typename OtherNumber>
      void
      add(const size_type    n_elements,
          const size_type *  indices,
          const OtherNumber *values);
 
      void
      sadd(const Number s, const Vector<Number, MemorySpace> &V);
 
 
 
 
      DEAL_II_DEPRECATED_EARLY
      size_type
      local_size() const;
 
      size_type
      locally_owned_size() const;
 
      bool
      in_local_range(const size_type global_index) const;
 
      iterator
      begin();
 
      const_iterator
      begin() const;
 
      iterator
      end();
 
      const_iterator
      end() const;
 
      Number
      operator()(const size_type global_index) const;
 
      Number &
      operator()(const size_type global_index);
 
      Number operator[](const size_type global_index) const;
      Number &operator[](const size_type global_index);
 
      Number
      local_element(const size_type local_index) const;
 
      Number &
      local_element(const size_type local_index);
 
      Number *
      get_values() const;
 
      template <typename OtherNumber>
      void
      extract_subvector_to(const std::vector<size_type> &indices,
                           std::vector<OtherNumber> &    values) const;
 
      template <typename ForwardIterator, typename OutputIterator>
      void
      extract_subvector_to(ForwardIterator       indices_begin,
                           const ForwardIterator indices_end,
                           OutputIterator        values_begin) const;
      virtual bool
      all_zero() const override;
 
      virtual Number
      mean_value() const override;
 
      real_type
      lp_norm(const real_type p) const;
 
 
      const MPI_Comm &
      get_mpi_communicator() const;
 
      const std::shared_ptr<const Utilities::MPI::Partitioner> &
      get_partitioner() const;
 
      bool
      partitioners_are_compatible(
        const Utilities::MPI::Partitioner &part) const;
 
      bool
      partitioners_are_globally_compatible(
        const Utilities::MPI::Partitioner &part) const;
 
      void
      set_ghost_state(const bool ghosted) const;
 
      const std::vector<ArrayView<const Number>> &
      shared_vector_data() const;
 
 
      DeclException0(ExcVectorTypeNotCompatible);
 
      DeclException0(ExcNotAllowedForCuda);
 
      DeclException3(ExcNonMatchingElements,
                     Number,
                     Number,
                     unsigned int,
                     << "Called compress(VectorOperation::insert), but"
                     << " the element received from a remote processor, value "
                     << std::setprecision(16) << arg1
                     << ", does not match with the value "
                     << std::setprecision(16) << arg2
                     << " on the owner processor " << arg3);
 
      DeclException4(
        ExcAccessToNonLocalElement,
        size_type,
        size_type,
        size_type,
        size_type,
        << "You tried to access element " << arg1
        << " of a distributed vector, but this element is not "
        << "stored on the current processor. Note: The range of "
        << "locally owned elements is [" << arg2 << "," << arg3
        << "], and there are " << arg4 << " ghost elements "
        << "that this vector can access."
        << "\n\n"
        << "A common source for this kind of problem is that you "
        << "are passing a 'fully distributed' vector into a function "
        << "that needs read access to vector elements that correspond "
        << "to degrees of freedom on ghost cells (or at least to "
        << "'locally active' degrees of freedom that are not also "
        << "'locally owned'). You need to pass a vector that has these "
        << "elements as ghost entries.");
 
    private:
      void
      add_local(const Number a, const VectorSpaceVector<Number> &V);
 
      void
      sadd_local(const Number                     s,
                 const Number                     a,
                 const VectorSpaceVector<Number> &V);
 
      template <typename Number2>
      Number
      inner_product_local(const Vector<Number2, MemorySpace> &V) const;
 
      real_type
      norm_sqr_local() const;
 
      Number
      mean_value_local() const;
 
      real_type
      l1_norm_local() const;
 
      real_type
      lp_norm_local(const real_type p) const;
 
      real_type
      linfty_norm_local() const;
 
      Number
      add_and_dot_local(const Number                       a,
                        const Vector<Number, MemorySpace> &V,
                        const Vector<Number, MemorySpace> &W);
 
      std::shared_ptr<const Utilities::MPI::Partitioner> partitioner;
 
      size_type allocated_size;
 
      mutable ::MemorySpace::MemorySpaceData<Number, MemorySpace> data;
 
      mutable std::shared_ptr<::parallel::internal::TBBPartitioner>
        thread_loop_partitioner;
 
      mutable ::MemorySpace::MemorySpaceData<Number, MemorySpace>
        import_data;
 
      mutable bool vector_is_ghosted;
 
#ifdef DEAL_II_WITH_MPI
      std::vector<MPI_Request> compress_requests;
 
      mutable std::vector<MPI_Request> update_ghost_values_requests;
#endif
 
      mutable std::mutex mutex;
 
      MPI_Comm comm_sm;
 
      void
      clear_mpi_requests();
 
      void
      resize_val(const size_type new_allocated_size,
                 const MPI_Comm &comm_sm = MPI_COMM_SELF);
 
      // Make all other vector types friends.
      template <typename Number2, typename MemorySpace2>
      friend class Vector;
 
      // Make BlockVector type friends.
      template <typename Number2>
      friend class BlockVector;
    };
    /*-------------------- Inline functions ---------------------------------*/
 
#ifndef DOXYGEN
 
    namespace internal
    {
      template <typename Number, typename MemorySpace>
      struct Policy
      {
        static inline typename Vector<Number, MemorySpace>::iterator
        begin(::MemorySpace::MemorySpaceData<Number, MemorySpace> &)
        {
          return nullptr;
        }
 
        static inline typename Vector<Number, MemorySpace>::const_iterator
        begin(
          const ::MemorySpace::MemorySpaceData<Number, MemorySpace> &)
        {
          return nullptr;
        }
 
        static inline Number *
        get_values(
          ::MemorySpace::MemorySpaceData<Number, MemorySpace> &)
        {
          return nullptr;
        }
      };
 
 
 
      template <typename Number>
      struct Policy<Number, ::MemorySpace::Host>
      {
        static inline
          typename Vector<Number, ::MemorySpace::Host>::iterator
          begin(::MemorySpace::
                  MemorySpaceData<Number, ::MemorySpace::Host> &data)
        {
          return data.values.get();
        }
 
        static inline
          typename Vector<Number, ::MemorySpace::Host>::const_iterator
          begin(const ::MemorySpace::
                  MemorySpaceData<Number, ::MemorySpace::Host> &data)
        {
          return data.values.get();
        }
 
        static inline Number *
        get_values(::MemorySpace::
                     MemorySpaceData<Number, ::MemorySpace::Host> &data)
        {
          return data.values.get();
        }
      };
 
 
 
      template <typename Number>
      struct Policy<Number, ::MemorySpace::CUDA>
      {
        static inline
          typename Vector<Number, ::MemorySpace::CUDA>::iterator
          begin(::MemorySpace::
                  MemorySpaceData<Number, ::MemorySpace::CUDA> &data)
        {
          return data.values_dev.get();
        }
 
        static inline
          typename Vector<Number, ::MemorySpace::CUDA>::const_iterator
          begin(const ::MemorySpace::
                  MemorySpaceData<Number, ::MemorySpace::CUDA> &data)
        {
          return data.values_dev.get();
        }
 
        static inline Number *
        get_values(::MemorySpace::
                     MemorySpaceData<Number, ::MemorySpace::CUDA> &data)
        {
          return data.values_dev.get();
        }
      };
    } // namespace internal
 
 
    template <typename Number, typename MemorySpace>
    inline bool
    Vector<Number, MemorySpace>::has_ghost_elements() const
    {
      return vector_is_ghosted;
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline typename Vector<Number, MemorySpace>::size_type
    Vector<Number, MemorySpace>::size() const
    {
      return partitioner->size();
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline typename Vector<Number, MemorySpace>::size_type
    Vector<Number, MemorySpace>::local_size() const
    {
      return locally_owned_size();
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline typename Vector<Number, MemorySpace>::size_type
    Vector<Number, MemorySpace>::locally_owned_size() const
    {
      return partitioner->locally_owned_size();
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline bool
    Vector<Number, MemorySpace>::in_local_range(
      const size_type global_index) const
    {
      return partitioner->in_local_range(global_index);
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline IndexSet
    Vector<Number, MemorySpace>::locally_owned_elements() const
    {
      IndexSet is(size());
 
      is.add_range(partitioner->local_range().first,
                   partitioner->local_range().second);
 
      return is;
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline typename Vector<Number, MemorySpace>::iterator
    Vector<Number, MemorySpace>::begin()
    {
      return internal::Policy<Number, MemorySpace>::begin(data);
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline typename Vector<Number, MemorySpace>::const_iterator
    Vector<Number, MemorySpace>::begin() const
    {
      return internal::Policy<Number, MemorySpace>::begin(data);
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline typename Vector<Number, MemorySpace>::iterator
    Vector<Number, MemorySpace>::end()
    {
      return internal::Policy<Number, MemorySpace>::begin(data) +
             partitioner->locally_owned_size();
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline typename Vector<Number, MemorySpace>::const_iterator
    Vector<Number, MemorySpace>::end() const
    {
      return internal::Policy<Number, MemorySpace>::begin(data) +
             partitioner->locally_owned_size();
    }
 
 
 
    template <typename Number, typename MemorySpace>
    const std::vector<ArrayView<const Number>> &
    Vector<Number, MemorySpace>::shared_vector_data() const
    {
      return data.values_sm;
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline Number
    Vector<Number, MemorySpace>::operator()(const size_type global_index) const
    {
      Assert((std::is_same<MemorySpace, ::MemorySpace::Host>::value),
             ExcMessage(
               "This function is only implemented for the Host memory space"));
      Assert(
        partitioner->in_local_range(global_index) ||
          partitioner->ghost_indices().is_element(global_index),
        ExcAccessToNonLocalElement(global_index,
                                   partitioner->local_range().first,
                                   partitioner->local_range().second - 1,
                                   partitioner->ghost_indices().n_elements()));
      // do not allow reading a vector which is not in ghost mode
      Assert(partitioner->in_local_range(global_index) ||
               vector_is_ghosted == true,
             ExcMessage("You tried to read a ghost element of this vector, "
                        "but it has not imported its ghost values."));
      return data.values[partitioner->global_to_local(global_index)];
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline Number &
    Vector<Number, MemorySpace>::operator()(const size_type global_index)
    {
      Assert((std::is_same<MemorySpace, ::MemorySpace::Host>::value),
             ExcMessage(
               "This function is only implemented for the Host memory space"));
      Assert(
        partitioner->in_local_range(global_index) ||
          partitioner->ghost_indices().is_element(global_index),
        ExcAccessToNonLocalElement(global_index,
                                   partitioner->local_range().first,
                                   partitioner->local_range().second - 1,
                                   partitioner->ghost_indices().n_elements()));
      // we would like to prevent reading ghosts from a vector that does not
      // have them imported, but this is not possible because we might be in a
      // part of the code where the vector has enabled ghosts but is non-const
      // (then, the compiler picks this method according to the C++ rule book
      // even if a human would pick the const method when this subsequent use
      // is just a read)
      return data.values[partitioner->global_to_local(global_index)];
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline Number Vector<Number, MemorySpace>::
                  operator[](const size_type global_index) const
    {
      return operator()(global_index);
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline Number &Vector<Number, MemorySpace>::
                   operator[](const size_type global_index)
    {
      return operator()(global_index);
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline Number
    Vector<Number, MemorySpace>::local_element(
      const size_type local_index) const
    {
      Assert((std::is_same<MemorySpace, ::MemorySpace::Host>::value),
             ExcMessage(
               "This function is only implemented for the Host memory space"));
      AssertIndexRange(local_index,
                       partitioner->locally_owned_size() +
                         partitioner->n_ghost_indices());
      // do not allow reading a vector which is not in ghost mode
      Assert(local_index < locally_owned_size() || vector_is_ghosted == true,
             ExcMessage("You tried to read a ghost element of this vector, "
                        "but it has not imported its ghost values."));
 
      return data.values[local_index];
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline Number &
    Vector<Number, MemorySpace>::local_element(const size_type local_index)
    {
      Assert((std::is_same<MemorySpace, ::MemorySpace::Host>::value),
             ExcMessage(
               "This function is only implemented for the Host memory space"));
 
      AssertIndexRange(local_index,
                       partitioner->locally_owned_size() +
                         partitioner->n_ghost_indices());
 
      return data.values[local_index];
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline Number *
    Vector<Number, MemorySpace>::get_values() const
    {
      return internal::Policy<Number, MemorySpace>::get_values(data);
    }
 
 
 
    template <typename Number, typename MemorySpace>
    template <typename OtherNumber>
    inline void
    Vector<Number, MemorySpace>::extract_subvector_to(
      const std::vector<size_type> &indices,
      std::vector<OtherNumber> &    values) const
    {
      for (size_type i = 0; i < indices.size(); ++i)
        values[i] = operator()(indices[i]);
    }
 
 
 
    template <typename Number, typename MemorySpace>
    template <typename ForwardIterator, typename OutputIterator>
    inline void
    Vector<Number, MemorySpace>::extract_subvector_to(
      ForwardIterator       indices_begin,
      const ForwardIterator indices_end,
      OutputIterator        values_begin) const
    {
      while (indices_begin != indices_end)
        {
          *values_begin = operator()(*indices_begin);
          indices_begin++;
          values_begin++;
        }
    }
 
 
 
    template <typename Number, typename MemorySpace>
    template <typename OtherNumber>
    inline void
    Vector<Number, MemorySpace>::add(
      const std::vector<size_type> &       indices,
      const ::Vector<OtherNumber> &values)
    {
      AssertDimension(indices.size(), values.size());
      for (size_type i = 0; i < indices.size(); ++i)
        {
          Assert(
            numbers::is_finite(values[i]),
            ExcMessage(
              "The given value is not finite but either infinite or Not A Number (NaN)"));
          this->operator()(indices[i]) += values(i);
        }
    }
 
 
 
    template <typename Number, typename MemorySpace>
    template <typename OtherNumber>
    inline void
    Vector<Number, MemorySpace>::add(const size_type    n_elements,
                                     const size_type *  indices,
                                     const OtherNumber *values)
    {
      for (size_type i = 0; i < n_elements; ++i, ++indices, ++values)
        {
          Assert(
            numbers::is_finite(*values),
            ExcMessage(
              "The given value is not finite but either infinite or Not A Number (NaN)"));
          this->operator()(*indices) += *values;
        }
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline const MPI_Comm &
    Vector<Number, MemorySpace>::get_mpi_communicator() const
    {
      return partitioner->get_mpi_communicator();
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline const std::shared_ptr<const Utilities::MPI::Partitioner> &
    Vector<Number, MemorySpace>::get_partitioner() const
    {
      return partitioner;
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline void
    Vector<Number, MemorySpace>::set_ghost_state(const bool ghosted) const
    {
      vector_is_ghosted = ghosted;
    }
 
#endif
 
  } // namespace distributed
} // namespace LinearAlgebra
 
 
template <typename Number, typename MemorySpace>
inline void
swap(LinearAlgebra::distributed::Vector<Number, MemorySpace> &u,
     LinearAlgebra::distributed::Vector<Number, MemorySpace> &v)
{
  u.swap(v);
}
 
 
template <typename Number, typename MemorySpace>
struct is_serial_vector<LinearAlgebra::distributed::Vector<Number, MemorySpace>>
  : std::false_type
{};
 
 
 
namespace internal
{
  namespace LinearOperatorImplementation
  {
    template <typename>
    class ReinitHelper;
 
    template <typename Number>
    class ReinitHelper<LinearAlgebra::distributed::Vector<Number>>
    {
    public:
      // A helper type-trait that leverage SFINAE to figure out if type T has
      // void T::get_mpi_communicator()
      template <typename T>
      struct has_get_mpi_communicator
      {
      private:
        static bool
        detect(...);
 
        template <typename U>
        static decltype(std::declval<U>().get_mpi_communicator())
        detect(const U &);
 
      public:
        static const bool value =
          !std::is_same<bool, decltype(detect(std::declval<T>()))>::value;
      };
 
      // A helper type-trait that leverage SFINAE to figure out if type T has
      // void T::locally_owned_domain_indices()
      template <typename T>
      struct has_locally_owned_domain_indices
      {
      private:
        static bool
        detect(...);
 
        template <typename U>
        static decltype(std::declval<U>().locally_owned_domain_indices())
        detect(const U &);
 
      public:
        static const bool value =
          !std::is_same<bool, decltype(detect(std::declval<T>()))>::value;
      };
 
      // A helper type-trait that leverage SFINAE to figure out if type T has
      // void T::locally_owned_range_indices()
      template <typename T>
      struct has_locally_owned_range_indices
      {
      private:
        static bool
        detect(...);
 
        template <typename U>
        static decltype(std::declval<U>().locally_owned_range_indices())
        detect(const U &);
 
      public:
        static const bool value =
          !std::is_same<bool, decltype(detect(std::declval<T>()))>::value;
      };
 
      // A helper type-trait that leverage SFINAE to figure out if type T has
      // void T::initialize_dof_vector(VectorType v)
      template <typename T>
      struct has_initialize_dof_vector
      {
      private:
        static bool
        detect(...);
 
        template <typename U>
        static decltype(std::declval<U>().initialize_dof_vector(
          std::declval<LinearAlgebra::distributed::Vector<Number> &>()))
        detect(const U &);
 
      public:
        static const bool value =
          !std::is_same<bool, decltype(detect(std::declval<T>()))>::value;
      };
 
      // Used for (Trilinos/PETSc)Wrappers::SparseMatrix
      template <typename MatrixType,
                typename std::enable_if<
                  has_get_mpi_communicator<MatrixType>::value &&
                    has_locally_owned_domain_indices<MatrixType>::value,
                  MatrixType>::type * = nullptr>
      static void
      reinit_domain_vector(MatrixType &                                mat,
                           LinearAlgebra::distributed::Vector<Number> &vec,
                           bool /*omit_zeroing_entries*/)
      {
        vec.reinit(mat.locally_owned_domain_indices(),
                   mat.get_mpi_communicator());
      }
 
      // Used for MatrixFree and DiagonalMatrix
      template <
        typename MatrixType,
        typename std::enable_if<has_initialize_dof_vector<MatrixType>::value,
                                MatrixType>::type * = nullptr>
      static void
      reinit_domain_vector(MatrixType &                                mat,
                           LinearAlgebra::distributed::Vector<Number> &vec,
                           bool omit_zeroing_entries)
      {
        mat.initialize_dof_vector(vec);
        if (!omit_zeroing_entries)
          vec = Number();
      }
 
      // Used for (Trilinos/PETSc)Wrappers::SparseMatrix
      template <typename MatrixType,
                typename std::enable_if<
                  has_get_mpi_communicator<MatrixType>::value &&
                    has_locally_owned_range_indices<MatrixType>::value,
                  MatrixType>::type * = nullptr>
      static void
      reinit_range_vector(MatrixType &                                mat,
                          LinearAlgebra::distributed::Vector<Number> &vec,
                          bool /*omit_zeroing_entries*/)
      {
        vec.reinit(mat.locally_owned_range_indices(),
                   mat.get_mpi_communicator());
      }
 
      // Used for MatrixFree and DiagonalMatrix
      template <
        typename MatrixType,
        typename std::enable_if<has_initialize_dof_vector<MatrixType>::value,
                                MatrixType>::type * = nullptr>
      static void
      reinit_range_vector(MatrixType &                                mat,
                          LinearAlgebra::distributed::Vector<Number> &vec,
                          bool omit_zeroing_entries)
      {
        mat.initialize_dof_vector(vec);
        if (!omit_zeroing_entries)
          vec = Number();
      }
    };
 
  } // namespace LinearOperatorImplementation
} /* namespace internal */
 
 
DEAL_II_NAMESPACE_CLOSE
 
#endif