la_parallel_vector.h

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

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

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);

      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);

#ifdef DEAL_II_WITH_PETSC

      DEAL_II_DEPRECATED
      Vector<Number, MemorySpace> &
      operator=(const PETScWrappers::MPI::Vector &petsc_vec);
#endif

#ifdef DEAL_II_WITH_TRILINOS

      DEAL_II_DEPRECATED
      Vector<Number, MemorySpace> &
      operator=(const TrilinosWrappers::MPI::Vector &trilinos_vec);
#endif


      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;

      void
      zero_out_ghosts() const;

      bool
      has_ghost_elements() const;

      template <typename Number2>
      void
      copy_locally_owned_data_from(const Vector<Number2, MemorySpace> &src);



      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 CommunicationPatternBase> communication_pattern =
          std::shared_ptr<const CommunicationPatternBase>()) 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
      void
      sadd(const Number                       s,
           const Number                       a,
           const Vector<Number, MemorySpace> &V,
           const Number                       b,
           const Vector<Number, MemorySpace> &W);

      DEAL_II_DEPRECATED
      void
      equ(const Number                       a,
          const Vector<Number, MemorySpace> &u,
          const Number                       b,
          const Vector<Number, MemorySpace> &v);




      size_type
      local_size() const;

      DEAL_II_DEPRECATED
      std::pair<size_type, size_type>
      local_range() const;

      DEAL_II_DEPRECATED
      bool
      in_local_range(const size_type global_index) const;

      DEAL_II_DEPRECATED
      size_type
      n_ghost_entries() const;

      DEAL_II_DEPRECATED
      const IndexSet &
      ghost_elements() const;

      DEAL_II_DEPRECATED
      bool
      is_ghost_entry(const types::global_dof_index 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;


      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 << " to " << arg3
                     << ", and there are " << arg4 << " ghost elements "
                     << "that this vector can access.");

    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 Threads::Mutex mutex;

      void
      clear_mpi_requests();

      void
      resize_val(const size_type new_allocated_size);

      /*
       * Make all other vector types friends.
       */
      template <typename Number2, typename MemorySpace2>
      friend class Vector;

      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 partitioner->local_size();
    }



    template <typename Number, typename MemorySpace>
    inline std::pair<typename Vector<Number, MemorySpace>::size_type,
                     typename Vector<Number, MemorySpace>::size_type>
    Vector<Number, MemorySpace>::local_range() const
    {
      return partitioner->local_range();
    }



    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>::size_type
    Vector<Number, MemorySpace>::n_ghost_entries() const
    {
      return partitioner->n_ghost_indices();
    }



    template <typename Number, typename MemorySpace>
    inline const IndexSet &
    Vector<Number, MemorySpace>::ghost_elements() const
    {
      return partitioner->ghost_indices();
    }



    template <typename Number, typename MemorySpace>
    inline bool
    Vector<Number, MemorySpace>::is_ghost_entry(
      const size_type global_index) const
    {
      return partitioner->is_ghost_entry(global_index);
    }



    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->local_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->local_size();
    }



    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,
                                   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,
                                   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->local_size() +
                         partitioner->n_ghost_indices());
      // do not allow reading a vector which is not in ghost mode
      Assert(local_index < local_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->local_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:
      template <typename Matrix>
      static void
      reinit_range_vector(const Matrix &                              matrix,
                          LinearAlgebra::distributed::Vector<Number> &v,
                          bool omit_zeroing_entries)
      {
        matrix.initialize_dof_vector(v);
        if (!omit_zeroing_entries)
          v = Number();
      }

      template <typename Matrix>
      static void
      reinit_domain_vector(const Matrix &                              matrix,
                           LinearAlgebra::distributed::Vector<Number> &v,
                           bool omit_zeroing_entries)
      {
        matrix.initialize_dof_vector(v);
        if (!omit_zeroing_entries)
          v = Number();
      }
    };

  } // namespace LinearOperatorImplementation
} /* namespace internal */


DEAL_II_NAMESPACE_CLOSE

#endif