la_parallel_vector.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2011 - 2018 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#ifndef dealii_la_parallel_vector_h
#define dealii_la_parallel_vector_h

#include <deal.II/base/config.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;
  }
}

namespace LinearAlgebra
{
  template <typename> class ReadWriteVector;
}

#ifdef DEAL_II_WITH_PETSC
namespace PETScWrappers
{
  namespace MPI
  {
    class Vector;
  }
}
#endif

#ifdef DEAL_II_WITH_TRILINOS
namespace TrilinosWrappers
{
  namespace MPI
  {
    class Vector;
  }
}
#endif

namespace LinearAlgebra
{
  namespace distributed
  {
    template <typename Number>
    class Vector : public ::LinearAlgebra::VectorSpaceVector<Number>,
      public Subscriptor
    {
    public:
      typedef Number                                            value_type;
      typedef value_type                                       *pointer;
      typedef const value_type                                 *const_pointer;
      typedef value_type                                       *iterator;
      typedef const value_type                                 *const_iterator;
      typedef value_type                                       &reference;
      typedef const value_type                                 &const_reference;
      typedef types::global_dof_index                           size_type;
      typedef typename numbers::NumberTraits<Number>::real_type real_type;

      Vector ();

      Vector (const Vector<Number> &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> &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> &v);

      Vector<Number> &
      operator = (const Vector<Number> &in_vector);

      template <typename Number2>
      Vector<Number> &
      operator = (const Vector<Number2> &in_vector);

#ifdef DEAL_II_WITH_PETSC

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

#ifdef DEAL_II_WITH_TRILINOS

      DEAL_II_DEPRECATED
      Vector<Number> &
      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> &src);



      virtual void reinit(const VectorSpaceVector<Number> &V,
                          const bool omit_zeroing_entries = false) override;

      virtual Vector<Number> &operator*= (const Number factor) override;

      virtual Vector<Number> &operator/= (const Number factor) override;

      virtual Vector<Number> &operator+= (const VectorSpaceVector<Number> &V) override;

      virtual Vector<Number> &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> &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> &V);

      DEAL_II_DEPRECATED
      void sadd (const Number          s,
                 const Number          a,
                 const Vector<Number> &V,
                 const Number          b,
                 const Vector<Number> &W);

      DEAL_II_DEPRECATED
      void equ (const Number a, const Vector<Number> &u,
                const Number b, const Vector<Number> &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);

      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;

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

      bool all_zero_local () const;

      template <typename Number2>
      Number inner_product_local (const Vector<Number2> &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> &V,
                                const Vector<Number> &W);

      std::shared_ptr<const Utilities::MPI::Partitioner> partitioner;

      size_type allocated_size;

      std::unique_ptr<Number[], decltype(&free)> values;

      mutable std::shared_ptr< ::parallel::internal::TBBPartitioner> thread_loop_partitioner;

      mutable std::unique_ptr<Number[]> 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> friend class Vector;

      template <typename Number2> friend class BlockVector;
    };
    /*----------------------- Inline functions ----------------------------------*/

#ifndef DOXYGEN


    template <typename Number>
    inline
    bool
    Vector<Number>::has_ghost_elements () const
    {
      return vector_is_ghosted;
    }



    template <typename Number>
    inline
    typename Vector<Number>::size_type
    Vector<Number>::size () const
    {
      return partitioner->size();
    }



    template <typename Number>
    inline
    typename Vector<Number>::size_type
    Vector<Number>::local_size () const
    {
      return partitioner->local_size();
    }



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



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



    template <typename Number>
    inline
    IndexSet
    Vector<Number>::locally_owned_elements() const
    {
      IndexSet is (size());

      is.add_range (partitioner->local_range().first,
                    partitioner->local_range().second);

      return is;
    }



    template <typename Number>
    inline
    typename Vector<Number>::size_type
    Vector<Number>::n_ghost_entries () const
    {
      return partitioner->n_ghost_indices();
    }



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



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



    template <typename Number>
    inline
    typename Vector<Number>::iterator
    Vector<Number>::begin ()
    {
      return values.get();
    }



    template <typename Number>
    inline
    typename Vector<Number>::const_iterator
    Vector<Number>::begin () const
    {
      return values.get();
    }



    template <typename Number>
    inline
    typename Vector<Number>::iterator
    Vector<Number>::end ()
    {
      return values.get()+partitioner->local_size();
    }



    template <typename Number>
    inline
    typename Vector<Number>::const_iterator
    Vector<Number>::end () const
    {
      return values.get()+partitioner->local_size();
    }



    template <typename Number>
    inline
    Number
    Vector<Number>::operator() (const size_type global_index) const
    {
      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 values[partitioner->global_to_local(global_index)];
    }



    template <typename Number>
    inline
    Number &
    Vector<Number>::operator() (const size_type global_index)
    {
      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 values[partitioner->global_to_local (global_index)];
    }



    template <typename Number>
    inline
    Number
    Vector<Number>::operator[] (const size_type global_index) const
    {
      return operator()(global_index);
    }



    template <typename Number>
    inline
    Number &
    Vector<Number>::operator[] (const size_type global_index)
    {
      return operator()(global_index);
    }



    template <typename Number>
    inline
    Number
    Vector<Number>::local_element (const size_type local_index) const
    {
      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 values[local_index];
    }



    template <typename Number>
    inline
    Number &
    Vector<Number>::local_element (const size_type local_index)
    {
      AssertIndexRange (local_index,
                        partitioner->local_size()+
                        partitioner->n_ghost_indices());
      return values[local_index];
    }



    template <typename Number>
    template <typename OtherNumber>
    inline
    void Vector<Number>::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>
    template <typename ForwardIterator, typename OutputIterator>
    inline
    void Vector<Number>::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>
    template <typename OtherNumber>
    inline
    void
    Vector<Number>::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>
    template <typename OtherNumber>
    inline
    void
    Vector<Number>::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>
    inline
    const MPI_Comm &
    Vector<Number>::get_mpi_communicator() const
    {
      return partitioner->get_mpi_communicator();
    }



    template <typename Number>
    inline
    const std::shared_ptr<const Utilities::MPI::Partitioner> &
    Vector<Number>::get_partitioner () const
    {
      return partitioner;
    }

#endif

  }
}


template <typename Number>
inline
void swap (LinearAlgebra::distributed::Vector<Number> &u,
           LinearAlgebra::distributed::Vector<Number> &v)
{
  u.swap (v);
}


template <typename Number>
struct is_serial_vector< LinearAlgebra::distributed::Vector< Number > > : 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 LinearOperator */
} /* namespace internal */


DEAL_II_NAMESPACE_CLOSE

#endif