la_parallel_vector.h

Go to the documentation of this file.

// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 2012 - 2024 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Part of the source code is dual licensed under Apache-2.0 WITH
// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
// governing the source code and code contributions can be found in
// LICENSE.md and CONTRIBUTING.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/enable_observer_pointer.h>
#include <deal.II/base/memory_space.h>
#include <deal.II/base/memory_space_data.h>
#include <deal.II/base/mpi_stub.h>
#include <deal.II/base/numbers.h>
#include <deal.II/base/parallel.h>
#include <deal.II/base/partitioner.h>
 
#include <deal.II/lac/read_vector.h>
#include <deal.II/lac/vector_operation.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 ::ReadVector<Number>
    {
    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_v<MemorySpace, ::MemorySpace::Host> ||
          std::is_same_v<MemorySpace, ::MemorySpace::Default>,
        "MemorySpace should be Host or Default");
 
      Vector();
 
      Vector(const Vector<Number, MemorySpace> &in_vector);
 
      Vector(Vector<Number, MemorySpace> &&in_vector); // NOLINT
 
      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);
 
      ~Vector();
 
      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 std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner,
        const bool                                                make_ghosted,
        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) noexcept;
 
      Vector<Number, MemorySpace> &
      operator=(Vector<Number, MemorySpace> &&in_vector); // NOLINT
 
      Vector<Number, MemorySpace> &
      operator=(const Vector<Number, MemorySpace> &in_vector);
 
      template <typename Number2>
      Vector<Number, MemorySpace> &
      operator=(const Vector<Number2, MemorySpace> &in_vector);
 
      void
      compress(VectorOperation::values operation);
 
      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;
 
      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_elements(const Vector<Number, MemorySpace2> &src,
                      VectorOperation::values             operation);
 
      template <typename MemorySpace2>
      DEAL_II_DEPRECATED void
      import(const Vector<Number, MemorySpace2> &src,
             VectorOperation::values             operation)
      {
        import_elements(src, operation);
      }
 
      Vector<Number, MemorySpace> &
      operator*=(const Number factor);
 
      Vector<Number, MemorySpace> &
      operator/=(const Number factor);
 
      Vector<Number, MemorySpace> &
      operator+=(const Vector<Number, MemorySpace> &V);
 
      Vector<Number, MemorySpace> &
      operator-=(const Vector<Number, MemorySpace> &V);
 
      void
      import_elements(
        const LinearAlgebra::ReadWriteVector<Number> &V,
        const VectorOperation::values                 operation,
        const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
          &communication_pattern = {});
 
      DEAL_II_DEPRECATED void
      import(const LinearAlgebra::ReadWriteVector<Number> &V,
             VectorOperation::values                       operation,
             std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
               communication_pattern = {})
      {
        import_elements(V, operation, communication_pattern);
      }
 
      Number
      operator*(const Vector<Number, MemorySpace> &V) const;
 
      void
      add(const Number a);
 
      void
      add(const Number a, const Vector<Number, MemorySpace> &V);
 
      void
      add(const Number                       a,
          const Vector<Number, MemorySpace> &V,
          const Number                       b,
          const Vector<Number, MemorySpace> &W);
 
      void
      add(const std::vector<size_type> &indices,
          const std::vector<Number>    &values);
 
      void
      sadd(const Number                       s,
           const Number                       a,
           const Vector<Number, MemorySpace> &V);
 
      void
      scale(const Vector<Number, MemorySpace> &scaling_factors);
 
      void
      equ(const Number a, const Vector<Number, MemorySpace> &V);
 
      real_type
      l1_norm() const;
 
      real_type
      l2_norm() const;
 
      real_type
      norm_sqr() const;
 
      real_type
      linfty_norm() const;
 
      Number
      add_and_dot(const Number                       a,
                  const Vector<Number, MemorySpace> &V,
                  const Vector<Number, MemorySpace> &W);
 
      virtual size_type
      size() const override;
 
      ::IndexSet
      locally_owned_elements() const;
 
      void
      print(std::ostream      &out,
            const unsigned int precision  = 3,
            const bool         scientific = true,
            const bool         across     = true) const;
 
      std::size_t
      memory_consumption() const;
      Vector<Number, MemorySpace> &
      operator=(const Number s);
 
      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);
 
      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;
 
      virtual void
      extract_subvector_to(
        const ArrayView<const types::global_dof_index> &indices,
        const ArrayView<Number> &elements) const override;
 
      template <typename ForwardIterator, typename OutputIterator>
      void
      extract_subvector_to(ForwardIterator       indices_begin,
                           const ForwardIterator indices_end,
                           OutputIterator        values_begin) const;
      bool
      all_zero() const;
 
      Number
      mean_value() const;
 
      real_type
      lp_norm(const real_type p) 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);
 
      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 Vector<Number, MemorySpace> &V);
 
      void
      sadd_local(const Number                       s,
                 const Number                       a,
                 const Vector<Number, MemorySpace> &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
 
    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>::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 data.values.data();
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline typename Vector<Number, MemorySpace>::const_iterator
    Vector<Number, MemorySpace>::begin() const
    {
      return data.values.data();
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline typename Vector<Number, MemorySpace>::iterator
    Vector<Number, MemorySpace>::end()
    {
      return data.values.data() + partitioner->locally_owned_size();
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline typename Vector<Number, MemorySpace>::const_iterator
    Vector<Number, MemorySpace>::end() const
    {
      return data.values.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_v<MemorySpace, ::MemorySpace::Host>),
             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 == 0 ?
                                     0 :
                                     (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_v<MemorySpace, ::MemorySpace::Host>),
             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 == 0 ?
                                     0 :
                                     (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_v<MemorySpace, ::MemorySpace::Host>),
             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_v<MemorySpace, ::MemorySpace::Host>),
             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 data.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 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) noexcept
{
  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>
      using get_mpi_communicator_t =
        decltype(std::declval<T>().get_mpi_communicator());
 
      template <typename T>
      static constexpr bool has_get_mpi_communicator =
        is_supported_operation<get_mpi_communicator_t, T>;
 
      // A helper type-trait that leverage SFINAE to figure out if type T has
      // void T::locally_owned_domain_indices()
      template <typename T>
      using locally_owned_domain_indices_t =
        decltype(std::declval<T>().locally_owned_domain_indices());
 
      template <typename T>
      static constexpr bool has_locally_owned_domain_indices =
        is_supported_operation<locally_owned_domain_indices_t, T>;
 
      // A helper type-trait that leverage SFINAE to figure out if type T has
      // void T::locally_owned_range_indices()
      template <typename T>
      using locally_owned_range_indices_t =
        decltype(std::declval<T>().locally_owned_range_indices());
 
      template <typename T>
      static constexpr bool has_locally_owned_range_indices =
        is_supported_operation<locally_owned_range_indices_t, T>;
 
      // A helper type-trait that leverage SFINAE to figure out if type T has
      // void T::initialize_dof_vector(VectorType v)
      template <typename T>
      using initialize_dof_vector_t =
        decltype(std::declval<T>().initialize_dof_vector(
          std::declval<LinearAlgebra::distributed::Vector<Number> &>()));
 
      template <typename T>
      static constexpr bool has_initialize_dof_vector =
        is_supported_operation<initialize_dof_vector_t, T>;
 
      // Used for (Trilinos/PETSc)Wrappers::SparseMatrix
      template <typename MatrixType,
#if !defined(__INTEL_COMPILER) || __INTEL_COMPILER >= 1900
                std::enable_if_t<has_get_mpi_communicator<MatrixType> &&
                                   has_locally_owned_domain_indices<MatrixType>,
#else
                // workaround for Intel 18
                std::enable_if_t<
                  is_supported_operation<get_mpi_communicator_t, MatrixType> &&
                    is_supported_operation<locally_owned_domain_indices_t,
                                           MatrixType>,
#endif
                                 MatrixType> * = 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,
#if !defined(__INTEL_COMPILER) || __INTEL_COMPILER >= 1900
                std::enable_if_t<has_initialize_dof_vector<MatrixType>,
#else
                // workaround for Intel 18
                std::enable_if_t<
                  is_supported_operation<initialize_dof_vector_t, MatrixType>,
#endif
                                 MatrixType> * = 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,
#if !defined(__INTEL_COMPILER) || __INTEL_COMPILER >= 1900
                std::enable_if_t<has_get_mpi_communicator<MatrixType> &&
                                   has_locally_owned_range_indices<MatrixType>,
#else
                // workaround for Intel 18
                std::enable_if_t<
                  is_supported_operation<get_mpi_communicator_t, MatrixType> &&
                    is_supported_operation<locally_owned_range_indices_t,
                                           MatrixType>,
#endif
                                 MatrixType> * = 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,
#if !defined(__INTEL_COMPILER) || __INTEL_COMPILER >= 1900
                std::enable_if_t<has_initialize_dof_vector<MatrixType>,
#else
                // workaround for Intel 18
                std::enable_if_t<
                  is_supported_operation<initialize_dof_vector_t, MatrixType>,
#endif
                                 MatrixType> * = 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