parallel_vector.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2011 - 2016 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__parallel_vector_h
#define dealii__parallel_vector_h

#include <deal.II/base/config.h>
#include <deal.II/base/index_set.h>
#include <deal.II/base/mpi.h>
#include <deal.II/base/template_constraints.h>
#include <deal.II/base/types.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/memory_consumption.h>
#include <deal.II/base/partitioner.h>
#include <deal.II/base/thread_management.h>
#include <deal.II/lac/vector_view.h>

#include <cstring>
#include <iomanip>


DEAL_II_NAMESPACE_OPEN

#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 parallel
{
  namespace distributed
  {
    template <typename Number> class BlockVector;

    template <typename Number>
    class Vector : 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;

      static const bool supports_distributed_data = true;

      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_cxx11::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> &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_cxx11::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

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

#ifdef DEAL_II_WITH_TRILINOS

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

      void copy_from (const Vector<Number> &in_vector,
                      const bool            call_update_ghost_values = false) DEAL_II_DEPRECATED;

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

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

      bool has_ghost_elements() const;

      bool all_zero () const;

      bool is_non_negative () const;

      template <typename Number2>
      bool operator == (const Vector<Number2> &v) const;

      template <typename Number2>
      bool operator != (const Vector<Number2> &v) const;

      template <typename Number2>
      Number operator * (const Vector<Number2> &V) const;

      real_type norm_sqr () const;

      Number mean_value () const;

      real_type l1_norm () const;

      real_type l2_norm () const;

      real_type lp_norm (const real_type p) const;

      real_type linfty_norm () const;

      Number add_and_dot (const Number          a,
                          const Vector<Number> &V,
                          const Vector<Number> &W);

      size_type size () const;

      size_type local_size() const;

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

      bool in_local_range (const size_type global_index) const;

      IndexSet locally_owned_elements () const;

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

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

      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;

      Number local_element (const size_type local_index) const;

      Number &local_element (const size_type local_index);



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

      Vector<Number> &operator -= (const Vector<Number> &V);

      template <typename OtherNumber>
      void add (const std::vector<size_type>   &indices,
                const std::vector<OtherNumber>  &values);

      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 add (const Number s);

      void add (const Vector<Number> &V) DEAL_II_DEPRECATED;

      void add (const Number a, const Vector<Number> &V);

      void add (const Number a, const Vector<Number> &V,
                const Number b, const Vector<Number> &W);

      void sadd (const Number          s,
                 const Vector<Number> &V);

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

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

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

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

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

      void scale (const Vector<Number> &scaling_factors);

      template <typename Number2>
      void scale (const Vector<Number2> &scaling_factors);

      void equ (const Number a, const Vector<Number> &u);

      template <typename Number2>
      void equ (const Number a, const Vector<Number2> &u);

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

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

      void ratio (const Vector<Number> &a,
                  const Vector<Number> &b) DEAL_II_DEPRECATED;


      const MPI_Comm &get_mpi_communicator () const;

      std_cxx11::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 print (std::ostream       &out,
                  const unsigned int  precision  = 3,
                  const bool          scientific = true,
                  const bool          across     = true) const;

      std::size_t memory_consumption () const;

      DeclException3 (ExcNonMatchingElements,
                      double, double, 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:
      bool all_zero_local () const;

      bool is_non_negative_local () const;

      template <typename Number2>
      bool vectors_equal_local (const Vector<Number2> &v) 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_cxx11::shared_ptr<const Utilities::MPI::Partitioner> partitioner;

      size_type allocated_size;

      Number         *val;

      mutable Number *import_data;

      mutable bool vector_is_ghosted;

      VectorView<Number> vector_view;

#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
    Vector<Number>::Vector ()
      :
      partitioner (new Utilities::MPI::Partitioner()),
      allocated_size (0),
      val (0),
      import_data (0),
      vector_is_ghosted (false),
      vector_view (0, static_cast<Number *>(0))
    {}



    template <typename Number>
    inline
    Vector<Number>::Vector (const Vector<Number> &v)
      :
      Subscriptor(),
      allocated_size (0),
      val (0),
      import_data (0),
      vector_is_ghosted (false),
      vector_view (0, static_cast<Number *>(0))
    {
      reinit (v, true);
      vector_view = v.vector_view;
      zero_out_ghosts();
    }



    template <typename Number>
    inline
    Vector<Number>::Vector (const IndexSet &local_range,
                            const IndexSet &ghost_indices,
                            const MPI_Comm  communicator)
      :
      allocated_size (0),
      val (0),
      import_data (0),
      vector_is_ghosted (false),
      vector_view (0, static_cast<Number *>(0))
    {
      reinit (local_range, ghost_indices, communicator);
    }



    template <typename Number>
    inline
    Vector<Number>::Vector (const IndexSet &local_range,
                            const MPI_Comm  communicator)
      :
      allocated_size (0),
      val (0),
      import_data (0),
      vector_is_ghosted (false),
      vector_view (0, static_cast<Number *>(0))
    {
      reinit (local_range, communicator);
    }



    template <typename Number>
    inline
    Vector<Number>::Vector (const size_type size)
      :
      allocated_size (0),
      val (0),
      import_data (0),
      vector_is_ghosted (false),
      vector_view (0, static_cast<Number *>(0))
    {
      reinit (size, false);
    }



    template <typename Number>
    inline
    Vector<Number>::
    Vector (const std_cxx11::shared_ptr<const Utilities::MPI::Partitioner> &partitioner)
      :
      allocated_size (0),
      val (0),
      import_data (0),
      vector_is_ghosted (false),
      vector_view (0, static_cast<Number *>(0))
    {
      reinit (partitioner);
    }



    template <typename Number>
    inline
    Vector<Number>::~Vector ()
    {
      resize_val(0);

      if (import_data != 0)
        delete[] import_data;
      import_data = 0;

      clear_mpi_requests();
    }



    template <typename Number>
    inline
    Vector<Number> &
    Vector<Number>::operator = (const Vector<Number> &c)
    {
#ifdef _MSC_VER
      return this->operator=<Number>(c);
#else
      return this->template operator=<Number>(c);
#endif
    }



    template <typename Number>
    template <typename Number2>
    inline
    Vector<Number> &
    Vector<Number>::operator = (const Vector<Number2> &c)
    {
      Assert (c.partitioner.get() != 0, ExcNotInitialized());

      // we update ghost values whenever one of the input or output vector
      // already held ghost values or when we import data from a vector with
      // the same local range but different ghost layout
      bool must_update_ghost_values = c.vector_is_ghosted;

      // check whether the two vectors use the same parallel partitioner. if
      // not, check if all local ranges are the same (that way, we can
      // exchange data between different parallel layouts). One variant which
      // is included here and necessary for compatibility with the other
      // distributed vector classes (Trilinos, PETSc) is the case when vector
      // c does not have any ghosts (constructed without ghost elements given)
      // but the current vector does: In that case, we need to exchange data
      // also when none of the two vector had updated its ghost values before.
      if (partitioner.get() == 0)
        reinit (c, true);
      else if (partitioner.get() != c.partitioner.get())
        {
          // local ranges are also the same if both partitioners are empty
          // (even if they happen to define the empty range as [0,0) or [c,c)
          // for some c!=0 in a different way).
          int local_ranges_are_identical =
            (local_range() == c.local_range() ||
             (local_range().second == local_range().first &&
              c.local_range().second == c.local_range().first));
          if ((c.partitioner->n_mpi_processes() > 1 &&
               Utilities::MPI::min(local_ranges_are_identical,
                                   c.partitioner->get_communicator()) == 0)
              ||
              !local_ranges_are_identical)
            reinit (c, true);
          else
            must_update_ghost_values |= vector_is_ghosted;

          must_update_ghost_values |=
            (c.partitioner->ghost_indices_initialized() == false &&
             partitioner->ghost_indices_initialized() == true);
        }
      else
        must_update_ghost_values |= vector_is_ghosted;

      // Need to explicitly downcast to ::Vector to make templated
      // operator= available.
      AssertDimension(vector_view.size(), c.vector_view.size());
      static_cast<::Vector<Number> &>(vector_view) = c.vector_view;

      if (must_update_ghost_values)
        update_ghost_values();
      else
        zero_out_ghosts();
      return *this;
    }



    template <typename Number>
    inline
    void
    Vector<Number>::compress (::VectorOperation::values operation)
    {
      compress_start (0, operation);
      compress_finish(operation);
    }



    template <typename Number>
    inline
    void
    Vector<Number>::update_ghost_values () const
    {
      update_ghost_values_start ();
      update_ghost_values_finish ();
    }



    template <typename Number>
    inline
    void
    Vector<Number>::zero_out_ghosts ()
    {
      std::fill_n (&val[partitioner->local_size()],
                   partitioner->n_ghost_indices(),
                   Number());
      vector_is_ghosted = false;
    }



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



    template <typename Number>
    inline
    bool
    Vector<Number>::all_zero_local () const
    {
      return partitioner->local_size()>0 ? vector_view.all_zero () : true;
    }



    template <typename Number>
    inline
    bool
    Vector<Number>::all_zero () const
    {
      // use int instead of bool. in order to make global reduction operations
      // work also when MPI_Init was not called, only call MPI_Allreduce
      // commands when there is more than one processor (note that reinit()
      // functions handle this case correctly through the job_supports_mpi()
      // query). this is the same in all the functions below
      int local_result = -static_cast<int>(all_zero_local());
      if (partitioner->n_mpi_processes() > 1)
        return -Utilities::MPI::max(local_result,
                                    partitioner->get_communicator());
      else
        return -local_result;
    }



    template <typename Number>
    inline
    bool
    Vector<Number>::is_non_negative_local () const
    {
      return partitioner->local_size()>0 ? vector_view.is_non_negative () : true;
    }



    template <typename Number>
    inline
    bool
    Vector<Number>::is_non_negative () const
    {
      int local_result = -static_cast<int>(is_non_negative_local());
      if (partitioner->n_mpi_processes() > 1)
        return -Utilities::MPI::max(local_result,
                                    partitioner->get_communicator());
      else
        return -local_result;
    }



    template <typename Number>
    template <typename Number2>
    inline
    bool
    Vector<Number>::vectors_equal_local (const Vector<Number2> &v) const
    {
      return partitioner->local_size()>0 ?
             vector_view.template operator == <Number2>(v.vector_view)
             : true;
    }



    template <typename Number>
    template <typename Number2>
    inline
    bool
    Vector<Number>::operator == (const Vector<Number2> &v) const
    {
      // MPI does not support bools, so use unsigned int instead. Two vectors
      // are equal if the check for non-equal fails on all processors
      unsigned int local_result = static_cast<int>(!vectors_equal_local(v));
      unsigned int result =
        partitioner->n_mpi_processes() > 1
        ?
        Utilities::MPI::max(local_result, partitioner->get_communicator())
        :
        local_result;
      return result==0;
    }



    template <typename Number>
    template <typename Number2>
    inline
    bool
    Vector<Number>::operator != (const Vector<Number2> &v) const
    {
      return !(operator == (v));
    }



    template <typename Number>
    template <typename Number2>
    inline
    Number
    Vector<Number>::inner_product_local(const Vector<Number2> &V) const
    {
      // on some processors, the size might be zero, which is not allowed by
      // the ::Vector class. Therefore, insert a check here
      return (partitioner->local_size()>0 ?
              vector_view.operator* (V.vector_view)
              : Number());
    }



    template <typename Number>
    template <typename Number2>
    inline
    Number
    Vector<Number>::operator * (const Vector<Number2> &V) const
    {
      Number local_result = inner_product_local(V);
      if (partitioner->n_mpi_processes() > 1)
        return Utilities::MPI::sum (local_result,
                                    partitioner->get_communicator());
      else
        return local_result;
    }



    template <typename Number>
    inline
    typename Vector<Number>::real_type
    Vector<Number>::norm_sqr_local () const
    {
      return partitioner->local_size()>0 ? vector_view.norm_sqr() : real_type();
    }



    template <typename Number>
    inline
    typename Vector<Number>::real_type
    Vector<Number>::norm_sqr () const
    {
      real_type local_result = norm_sqr_local();
      if (partitioner->n_mpi_processes() > 1)
        return Utilities::MPI::sum(local_result,
                                   partitioner->get_communicator());
      else
        return local_result;
    }



    template <typename Number>
    inline
    Number
    Vector<Number>::mean_value_local () const
    {
      Assert (partitioner->size()!=0, ExcEmptyObject());
      return (partitioner->local_size() ?
              vector_view.mean_value()
              : Number());
    }



    template <typename Number>
    inline
    Number
    Vector<Number>::mean_value () const
    {
      Number local_result = mean_value_local();
      if (partitioner->n_mpi_processes() > 1)
        return Utilities::MPI::sum (local_result *
                                    (real_type)partitioner->local_size(),
                                    partitioner->get_communicator())
               /(real_type)partitioner->size();
      else
        return local_result;
    }



    template <typename Number>
    inline
    typename Vector<Number>::real_type
    Vector<Number>::l1_norm_local () const
    {
      return partitioner->local_size() ? vector_view.l1_norm() : real_type();
    }



    template <typename Number>
    inline
    typename Vector<Number>::real_type
    Vector<Number>::l1_norm () const
    {
      real_type local_result = l1_norm_local();
      if (partitioner->n_mpi_processes() > 1)
        return Utilities::MPI::sum(local_result,
                                   partitioner->get_communicator());
      else
        return local_result;
    }



    template <typename Number>
    inline
    typename Vector<Number>::real_type
    Vector<Number>::l2_norm () const
    {
      return std::sqrt(norm_sqr());
    }



    template <typename Number>
    inline
    typename Vector<Number>::real_type
    Vector<Number>::lp_norm_local (const real_type p) const
    {
      return partitioner->local_size() ? vector_view.lp_norm(p) : real_type();
    }



    template <typename Number>
    inline
    typename Vector<Number>::real_type
    Vector<Number>::lp_norm (const real_type p) const
    {
      const real_type local_result = lp_norm_local(p);
      if (partitioner->n_mpi_processes() > 1)
        return std::pow (Utilities::MPI::sum(std::pow(local_result,p),
                                             partitioner->get_communicator()),
                         static_cast<real_type>(1.0/p));
      else
        return local_result;
    }



    template <typename Number>
    inline
    typename Vector<Number>::real_type
    Vector<Number>::linfty_norm_local () const
    {
      return partitioner->local_size() ? vector_view.linfty_norm() : real_type();
    }



    template <typename Number>
    inline
    typename Vector<Number>::real_type
    Vector<Number>::linfty_norm () const
    {
      const real_type local_result = linfty_norm_local();
      if (partitioner->n_mpi_processes() > 1)
        return Utilities::MPI::max (local_result,
                                    partitioner->get_communicator());
      else
        return local_result;
    }



    template <typename Number>
    inline
    Number
    Vector<Number>::add_and_dot_local(const Number          a,
                                      const Vector<Number> &V,
                                      const Vector<Number> &W)
    {
      // on some processors, the size might be zero, which is not allowed by
      // the ::Vector class. Therefore, insert a check here
      return (partitioner->local_size()>0 ?
              vector_view.add_and_dot(a, V.vector_view, W.vector_view)
              : Number());
    }



    template <typename Number>
    inline
    Number
    Vector<Number>::add_and_dot (const Number          a,
                                 const Vector<Number> &V,
                                 const Vector<Number> &W)
    {
      Number local_result = add_and_dot_local(a, V, W);
      if (partitioner->n_mpi_processes() > 1)
        return Utilities::MPI::sum (local_result,
                                    partitioner->get_communicator());
      else
        return local_result;
    }



    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 (local_range().first, 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 vector_view.begin();
    }



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



    template <typename Number>
    inline
    typename Vector<Number>::iterator
    Vector<Number>::end ()
    {
      return vector_view.end();
    }



    template <typename Number>
    inline
    typename Vector<Number>::const_iterator
    Vector<Number>::end () const
    {
      return vector_view.end();
    }



    template <typename Number>
    inline
    Number
    Vector<Number>::operator() (const size_type global_index) const
    {
      Assert (in_local_range (global_index) ||
              partitioner->ghost_indices().is_element(global_index),
              ExcAccessToNonLocalElement(global_index, local_range().first,
                                         local_range().second,
                                         partitioner->ghost_indices().n_elements()));
      // do not allow reading a vector which is not in ghost mode
      Assert (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 val[partitioner->global_to_local(global_index)];
    }



    template <typename Number>
    inline
    Number &
    Vector<Number>::operator() (const size_type global_index)
    {
      Assert (in_local_range (global_index) ||
              partitioner->ghost_indices().is_element(global_index),
              ExcAccessToNonLocalElement(global_index, local_range().first,
                                         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 val[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>
    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>
    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 val[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 val[local_index];
    }



    template <typename Number>
    inline
    Vector<Number> &
    Vector<Number>::operator = (const Number s)
    {
      // if we call Vector::operator=0, we want to zero out all the entries
      // plus ghosts.
      if (partitioner->local_size() > 0)
        vector_view.::template Vector<Number>::operator= (s);
      if (s==Number())
        zero_out_ghosts();

      return *this;
    }



    template <typename Number>
    inline
    Vector<Number> &
    Vector<Number>::operator += (const Vector<Number> &v)
    {
      AssertDimension (local_size(), v.local_size());

      // ::Vector does not allow empty fields but this might happen on
      // some processors for parallel implementation
      if (local_size()>0)
        vector_view += v.vector_view;

      if (vector_is_ghosted)
        update_ghost_values();

      return *this;
    }



    template <typename Number>
    inline
    Vector<Number> &
    Vector<Number>::operator -= (const Vector<Number> &v)
    {
      AssertDimension (local_size(), v.local_size());

      // ::Vector does not allow empty fields but this might happen on
      // some processors for parallel implementation
      if (local_size()>0)
        vector_view -= v.vector_view;

      if (vector_is_ghosted)
        update_ghost_values();

      return *this;
    }



    template <typename Number>
    template <typename OtherNumber>
    inline
    void
    Vector<Number>::add (const std::vector<size_type> &indices,
                         const std::vector<OtherNumber>  &values)
    {
      AssertDimension (indices.size(), values.size());
      add (indices.size(), &indices[0], &values[0]);
    }



    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());
      add (indices.size(), &indices[0], values.begin());
    }



    template <typename Number>
    template <typename OtherNumber>
    inline
    void
    Vector<Number>::add (const size_type    n_indices,
                         const size_type   *indices,
                         const OtherNumber *values)
    {
      for (size_type i=0; i<n_indices; ++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>
    inline
    void
    Vector<Number>::add (const Number a)
    {
      // ::Vector does not allow empty fields but this might happen on
      // some processors for parallel implementation
      if (local_size())
        vector_view.add (a);

      if (vector_is_ghosted)
        update_ghost_values();
    }



    template <typename Number>
    inline
    void
    Vector<Number>::add (const Vector<Number> &v)
    {
      // ::Vector does not allow empty fields but this might happen on
      // some processors for parallel implementation
      if (local_size())
        vector_view.add (v.vector_view);

      if (vector_is_ghosted)
        update_ghost_values();
    }



    template <typename Number>
    inline
    void
    Vector<Number>::add (const Number a,
                         const Vector<Number> &v)
    {
      // ::Vector does not allow empty fields but this might happen on
      // some processors for parallel implementation
      if (local_size())
        vector_view.add (a, v.vector_view);

      if (vector_is_ghosted)
        update_ghost_values();
    }



    template <typename Number>
    inline
    void
    Vector<Number>::add (const Number a,
                         const Vector<Number> &v,
                         const Number b,
                         const Vector<Number> &w)
    {
      // ::Vector does not allow empty fields but this might happen on
      // some processors for parallel implementation
      if (local_size())
        vector_view.add (a, v.vector_view, b, w.vector_view);

      if (vector_is_ghosted)
        update_ghost_values();
    }



    template <typename Number>
    inline
    void
    Vector<Number>::sadd (const Number x,
                          const Vector<Number> &v)
    {
      // ::Vector does not allow empty fields but this might happen on
      // some processors for parallel implementation
      if (local_size())
        vector_view.sadd (x, v.vector_view);

      if (vector_is_ghosted)
        update_ghost_values();
    }



    template <typename Number>
    inline
    void
    Vector<Number>::sadd (const Number x,
                          const Number a,
                          const Vector<Number> &v)
    {
      // ::Vector does not allow empty fields but this might happen on
      // some processors for parallel implementation
      if (local_size())
        vector_view.sadd (x, a, v.vector_view);

      if (vector_is_ghosted)
        update_ghost_values();
    }



    template <typename Number>
    inline
    void
    Vector<Number>::sadd (const Number x,
                          const Number a,
                          const Vector<Number> &v,
                          const Number b,
                          const Vector<Number> &w)
    {
      // ::Vector does not allow empty fields but this might happen on
      // some processors for parallel implementation
      if (local_size())
        vector_view.sadd (x, a, v.vector_view, b, w.vector_view);

      if (vector_is_ghosted)
        update_ghost_values();
    }



    template <typename Number>
    inline
    void
    Vector<Number>::sadd (const Number s,
                          const Number a,
                          const Vector<Number> &v,
                          const Number b,
                          const Vector<Number> &w,
                          const Number c,
                          const Vector<Number> &x)
    {
      // ::Vector does not allow empty fields but this might happen on
      // some processors for parallel implementation
      if (local_size())
        vector_view.sadd (s, a, v.vector_view, b, w.vector_view,
                          c, x.vector_view);

      if (vector_is_ghosted)
        update_ghost_values();
    }



    template <typename Number>
    inline
    Vector<Number> &
    Vector<Number>::operator *= (const Number factor)
    {
      // ::Vector does not allow empty fields but this might happen on
      // some processors for parallel implementation
      if (local_size())
        vector_view *= factor;

      if (vector_is_ghosted)
        update_ghost_values();

      return *this;
    }



    template <typename Number>
    inline
    Vector<Number> &
    Vector<Number>::operator /= (const Number factor)
    {
      operator *= (1./factor);
      return *this;
    }



    template <typename Number>
    inline
    void
    Vector<Number>::scale (const Vector<Number> &scaling_factors)
    {
      // ::Vector does not allow empty fields but this might happen on
      // some processors for parallel implementation
      if (local_size())
        vector_view.scale (scaling_factors.vector_view);

      if (vector_is_ghosted)
        update_ghost_values();
    }



    template <typename Number>
    template <typename Number2>
    inline
    void
    Vector<Number>::scale (const Vector<Number2> &scaling_factors)
    {
      if (local_size())
        vector_view.template scale<Number2> (scaling_factors.vector_view);

      if (vector_is_ghosted)
        update_ghost_values();
    }



    template <typename Number>
    inline
    void
    Vector<Number>::equ (const Number a,
                         const Vector<Number> &v)
    {
      // ::Vector does not allow empty fields but this might happen on
      // some processors for parallel implementation
      if (local_size())
        vector_view.equ (a, v.vector_view);

      if (vector_is_ghosted)
        update_ghost_values();
    }



    template <typename Number>
    template <typename Number2>
    inline
    void
    Vector<Number>::equ (const Number a,
                         const Vector<Number2> &v)
    {
      // ::Vector does not allow empty fields but this might happen on
      // some processors for parallel implementation
      if (local_size())
        vector_view.equ (a, v.vector_view);

      if (vector_is_ghosted)
        update_ghost_values();
    }



    template <typename Number>
    inline
    void
    Vector<Number>::equ (const Number a,
                         const Vector<Number> &v,
                         const Number b,
                         const Vector<Number> &w)
    {
      // ::Vector does not allow empty fields but this might happen on
      // some processors for parallel implementation
      if (local_size())
        vector_view.equ (a, v.vector_view, b, w.vector_view);

      if (vector_is_ghosted)
        update_ghost_values();
    }



    template <typename Number>
    inline
    void
    Vector<Number>::equ (const Number a,
                         const Vector<Number> &v,
                         const Number b,
                         const Vector<Number> &w,
                         const Number c,
                         const Vector<Number> &x)
    {
      // ::Vector does not allow empty fields but this might happen on
      // some processors for parallel implementation
      if (local_size())
        vector_view.equ (a, v.vector_view, b, w.vector_view,
                         c, x.vector_view);

      if (vector_is_ghosted)
        update_ghost_values();
    }



    template <typename Number>
    inline
    void
    Vector<Number>::ratio (const Vector<Number> &a,
                           const Vector<Number> &b)
    {
      // ::Vector does not allow empty fields but this might happen on
      // some processors for parallel implementation
      if (local_size())
        vector_view.ratio (a.vector_view, b.vector_view);

      if (vector_is_ghosted)
        update_ghost_values();
    }



    template <typename Number>
    inline
    const MPI_Comm &
    Vector<Number>::get_mpi_communicator() const
    {
      return partitioner->get_communicator();
    }



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

#endif  // ifndef DOXYGEN

  } // end of namespace distributed

} // end of namespace parallel



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


DEAL_II_NAMESPACE_CLOSE

#endif