trilinos_tpetra_vector.h

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
// Copyright (C) 2018 - 2026 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Detailed license information governing the source code and contributions
// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
//
// -----------------------------------------------------------------------------
 
#ifndef dealii_trilinos_tpetra_vector_h
#define dealii_trilinos_tpetra_vector_h
 
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/types.h>
 
#include <deal.II/lac/trilinos_tpetra_types.h>
 
#ifdef DEAL_II_TRILINOS_WITH_TPETRA
 
#  include <deal.II/base/index_set.h>
#  include <deal.II/base/mpi_stub.h>
 
#  include <deal.II/lac/read_vector.h>
#  include <deal.II/lac/trilinos_tpetra_communication_pattern.h>
#  include <deal.II/lac/vector.h>
#  include <deal.II/lac/vector_operation.h>
#  include <deal.II/lac/vector_type_traits.h>
 
#  include <Teuchos_Comm.hpp>
#  include <Teuchos_OrdinalTraits.hpp>
#  include <Tpetra_Core.hpp>
#  include <Tpetra_Vector.hpp>
#  include <Tpetra_Version.hpp>
 
#  include <complex>
#  include <memory>
#  include <optional>
#  include <utility>
 
#endif
 
DEAL_II_NAMESPACE_OPEN
 
#ifdef DEAL_II_TRILINOS_WITH_TPETRA
template <typename Number>
struct is_tpetra_type : std::false_type
{};
 
#  ifdef DEAL_II_TRILINOS_WITH_TPETRA_INST_FLOAT
template <>
struct is_tpetra_type<float> : std::true_type
{};
#  endif
 
#  ifdef DEAL_II_TRILINOS_WITH_TPETRA_INST_DOUBLE
template <>
struct is_tpetra_type<double> : std::true_type
{};
#  endif
 
#  ifdef DEAL_II_WITH_COMPLEX_VALUES
#    ifdef DEAL_II_TRILINOS_WITH_TPETRA_INST_COMPLEX_FLOAT
template <>
struct is_tpetra_type<std::complex<float>> : std::true_type
{};
#    endif
 
#    ifdef DEAL_II_TRILINOS_WITH_TPETRA_INST_COMPLEX_DOUBLE
template <>
struct is_tpetra_type<std::complex<double>> : std::true_type
{};
#    endif
#  endif
 
namespace LinearAlgebra
{
  // Forward declaration
#  ifndef DOXYGEN
  template <typename Number>
  class ReadWriteVector;
#  endif


  namespace TpetraWrappers
  {

    class VectorTraits
    {
    public:
      using size_type = types::global_dof_index;
    };


    namespace internal
    {
      using size_type = types::global_dof_index;

      template <typename Number,
                typename MemorySpace = ::MemorySpace::Host>
      class VectorReference
      {
      private:
        VectorReference(Vector<Number, MemorySpace> &vector,
                        const size_type              index);
 
      public:
        VectorReference(const VectorReference &) = default;

        const VectorReference &
        operator=(const VectorReference &r) const;

        VectorReference &
        operator=(const VectorReference &r);

        const VectorReference &
        operator=(const Number &s) const;

        const VectorReference &
        operator+=(const Number &s) const;

        const VectorReference &
        operator-=(const Number &s) const;

        const VectorReference &
        operator*=(const Number &s) const;

        const VectorReference &
        operator/=(const Number &s) const;

        operator Number() const;

        DeclException1(ExcTrilinosError,
                       int,
                       << "An error with error number " << arg1
                       << " occurred while calling a Trilinos function");
 
        /*
         * Access to a an element that is not (locally-)owned.
         *
         * @ingroup Exceptions
         */
        DeclException4(
          ExcAccessToNonLocalElement,
          size_type,
          size_type,
          size_type,
          size_type,
          << "You are trying to access element " << arg1
          << " of a distributed vector, but this element is not stored "
          << "on the current processor. Note: There are " << arg2
          << " elements stored "
          << "on the current processor from within the range [" << arg3 << ','
          << arg4 << "] but Trilinos vectors need not store contiguous "
          << "ranges on each processor, and not every element in "
          << "this range may in fact be stored locally."
          << "\n\n"
          << "A common source for this kind of problem is that you "
          << "are passing a 'fully distributed' vector into a function "
          << "that needs read access to vector elements that correspond "
          << "to degrees of freedom on ghost cells (or at least to "
          << "'locally active' degrees of freedom that are not also "
          << "'locally owned'). You need to pass a vector that has these "
          << "elements as ghost entries.");
 
      private:
        Vector<Number, MemorySpace> &vector;

        const size_type index;
 
        // Make the vector class a friend, so that it can create objects of the
        // present type.
        friend class Vector<Number, MemorySpace>;
      }; // class VectorReference
 
    } // namespace internal
 

    template <typename Number, typename MemorySpace = ::MemorySpace::Host>
    class Vector : public ReadVector<Number>
    {
    public:
      using value_type = Number;
      using real_type  = typename numbers::NumberTraits<Number>::real_type;
      using size_type  = types::global_dof_index;
      using reference  = internal::VectorReference<Number, MemorySpace>;
      using const_reference =
        const internal::VectorReference<Number, MemorySpace>;

      Vector();

      Vector(const Vector &V);

      Vector(
        const Teuchos::RCP<TpetraTypes::VectorType<Number, MemorySpace>> V);

      explicit Vector(const IndexSet &parallel_partitioner,
                      const MPI_Comm  communicator = MPI_COMM_WORLD);

      explicit Vector(const IndexSet &locally_owned_entries,
                      const IndexSet &ghost_entries,
                      const MPI_Comm  communicator,
                      const bool      vector_writable = false);

      void
      clear();

      void
      reinit(const IndexSet &parallel_partitioner,
             const MPI_Comm  communicator         = MPI_COMM_WORLD,
             const bool      omit_zeroing_entries = false);

      void
      reinit(const IndexSet &locally_owned_entries,
             const IndexSet &locally_relevant_or_ghost_entries,
             const MPI_Comm  communicator    = MPI_COMM_WORLD,
             const bool      vector_writable = false);

      void
      reinit(const Vector<Number, MemorySpace> &V,
             const bool                         omit_zeroing_entries = false);

      virtual void
      swap(Vector &v) noexcept;

      virtual void
      extract_subvector_to(
        const ArrayView<const types::global_dof_index> &indices,
        const ArrayView<Number> &elements) const override;

      Vector &
      operator=(const Vector &V);

      template <typename OtherNumber>
      Vector &
      operator=(const ::Vector<OtherNumber> &V);

      Vector &
      operator=(const Number s);

      void
      import_elements(
        const ReadWriteVector<Number> &V,
        VectorOperation::values        operation,
        const Teuchos::RCP<const Utilities::MPI::CommunicationPatternBase>
          &communication_pattern);
 
      /*
       * Imports all the elements present in the vector's IndexSet from the
       * input vector @p V. VectorOperation::values @p operation is used to decide if
       * the elements in @p V should be added to the current vector or replace the
       * current elements.
       */
      void
      import_elements(const ReadWriteVector<Number> &V,
                      VectorOperation::values        operation);

 


      reference
      operator()(const size_type index);

      Number
      operator()(const size_type index) const;

      reference
      operator[](const size_type index);

      Number
      operator[](const size_type index) const;

 


      Vector &
      operator*=(const Number factor);

      Vector &
      operator/=(const Number factor);

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

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

      Number
      operator*(const Vector<Number, MemorySpace> &V) const;

      void
      add(const Number a);

      void
      add(const Vector<Number, MemorySpace> &V,
          const bool                         allow_different_maps = false);

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

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

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

      void
      add(const size_type  n_elements,
          const size_type *indices,
          const Number    *values);

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

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

      void
      set(const size_type  n_elements,
          const size_type *indices,
          const Number    *values);

      void
      set(const std::vector<size_type> &indices,
          const std::vector<Number>    &values);

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

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

      bool
      all_zero() const;

      bool
      is_non_negative() const;

 


      Number
      mean_value() const;

      Number
      min() const;

      Number
      max() const;

      real_type
      l1_norm() const;

      real_type
      l2_norm() const;

      real_type
      linfty_norm() const;

      real_type
      norm_sqr() const;

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

 


      bool
      has_ghost_elements() const;

      bool
      operator==(const Vector<Number, MemorySpace> &v) const;

      bool
      operator!=(const Vector<Number, MemorySpace> &v) const;

      virtual size_type
      size() const override;

      size_type
      locally_owned_size() const;

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

      bool
      in_local_range(const size_type index) const;

      bool
      is_compressed() const;

      MPI_Comm
      get_mpi_communicator() const;

      ::IndexSet
      locally_owned_elements() const;

 

      void
      compress(const VectorOperation::values operation);

      void
      update_ghost_values() const;

      const TpetraTypes::VectorType<Number, MemorySpace> &
      trilinos_vector() const;

      TpetraTypes::VectorType<Number, MemorySpace> &
      trilinos_vector();

      const TpetraTypes::MapType<MemorySpace> &
      trilinos_partitioner() const;

      Teuchos::RCP<const TpetraTypes::VectorType<Number, MemorySpace>>
      trilinos_rcp() const;

      Teuchos::RCP<TpetraTypes::VectorType<Number, MemorySpace>>
      trilinos_rcp();

      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;

      MPI_Comm
      mpi_comm() const;

 

      DeclException0(ExcDifferentParallelPartitioning);

      DeclException0(ExcVectorTypeNotCompatible);
 
      /*
       * Access to a an element that is not (locally-)owned.
       *
       * @ingroup Exceptions
       */
      DeclException4(
        ExcAccessToNonLocalElement,
        size_type,
        size_type,
        size_type,
        size_type,
        << "You are trying to access element " << arg1
        << " of a distributed vector, but this element is not stored "
        << "on the current processor. Note: There are " << arg2
        << " elements stored "
        << "on the current processor from within the range [" << arg3 << ','
        << arg4 << "] but Trilinos vectors need not store contiguous "
        << "ranges on each processor, and not every element in "
        << "this range may in fact be stored locally."
        << "\n\n"
        << "A common source for this kind of problem is that you "
        << "are passing a 'fully distributed' vector into a function "
        << "that needs read access to vector elements that correspond "
        << "to degrees of freedom on ghost cells (or at least to "
        << "'locally active' degrees of freedom that are not also "
        << "'locally owned'). You need to pass a vector that has these "
        << "elements as ghost entries.");

      DeclExceptionMsg(ExcMissingIndexSet,
                       "To compress a vector, a locally_relevant_dofs "
                       "index set, and a locally_owned_dofs index set "
                       "must be provided. These index sets must be "
                       "provided either when the vector is initialized "
                       "or when compress is called. See the documentation "
                       "of compress() for more information.");

      DeclException1(ExcTrilinosError,
                     int,
                     << "An error with error number " << arg1
                     << " occurred while calling a Trilinos function");
 
    private:
      void
      create_tpetra_comm_pattern(const IndexSet &source_index_set,
                                 const MPI_Comm  mpi_comm);

      bool compressed;

      bool has_ghost;

      VectorOperation::values last_action;

      Teuchos::RCP<TpetraTypes::VectorType<Number, MemorySpace>> vector;

      Teuchos::RCP<TpetraTypes::VectorType<Number, MemorySpace>>
        nonlocal_vector;

      ::IndexSet source_stored_elements;

      ::IndexSet local_entries;
 
      using dual_view_type =
        typename TpetraTypes::VectorType<Number, MemorySpace>::dual_view_type;
      using host_view_type = typename dual_view_type::t_host;
 
      using array_view_type =
        Kokkos::Subview<host_view_type,
                        std::remove_const_t<decltype(Kokkos::ALL)>,
                        unsigned>;
      std::vector<types::global_dof_index> nonlocal_cached_indices;

      std::vector<Number> nonlocal_cached_values;

      mutable std::optional<array_view_type> vector_1d_view;

      mutable std::optional<array_view_type> nonlocal_vector_1d_view;

      mutable std::mutex mutex;

      Teuchos::RCP<const TpetraWrappers::CommunicationPattern<MemorySpace>>
        tpetra_comm_pattern;
 
      // Make the reference class a friend.
      friend class internal::VectorReference<Number, MemorySpace>;
    };
 
 
    /* ------------------------- Inline functions ---------------------- */
 
    template <typename Number, typename MemorySpace>
    inline void
    swap(Vector<Number, MemorySpace> &u,
         Vector<Number, MemorySpace> &v) noexcept
    {
      u.swap(v);
    }
 
 
    template <typename Number, typename MemorySpace>
    inline bool
    Vector<Number, MemorySpace>::has_ghost_elements() const
    {
      return has_ghost;
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline bool
    Vector<Number, MemorySpace>::is_compressed() const
    {
      return compressed;
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline void
    Vector<Number, MemorySpace>::swap(Vector<Number, MemorySpace> &v) noexcept
    {
      std::swap(compressed, v.compressed);
      std::swap(has_ghost, v.has_ghost);
      std::swap(last_action, v.last_action);
      vector.swap(v.vector);
      nonlocal_vector.swap(v.nonlocal_vector);
      std::swap(vector_1d_view, v.vector_1d_view);
      std::swap(nonlocal_vector_1d_view, v.nonlocal_vector_1d_view);
      std::swap(source_stored_elements, v.source_stored_elements);
      std::swap(local_entries, v.local_entries);
      std::swap(nonlocal_cached_indices, v.nonlocal_cached_indices);
      std::swap(nonlocal_cached_values, v.nonlocal_cached_values);
      tpetra_comm_pattern.swap(v.tpetra_comm_pattern);
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline void
    Vector<Number, MemorySpace>::add(const Vector<Number, MemorySpace> &V,
                                     const bool allow_different_maps)
    {
      Assert(!allow_different_maps, ExcNotImplemented());
      this->add(1, V);
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline void
    Vector<Number, MemorySpace>::add(const std::vector<size_type> &indices,
                                     const std::vector<Number>    &values)
    {
      // if we have ghost values, do not allow
      // writing to this vector at all.
      AssertDimension(indices.size(), values.size());
 
      add(indices.size(), indices.data(), values.data());
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline void
    Vector<Number, MemorySpace>::add(const size_type  n_elements,
                                     const size_type *indices,
                                     const Number    *values)
    {
      // if we have ghost values, do not allow
      // writing to this vector at all.
      Assert(!has_ghost_elements(), ExcGhostsPresent());
 
      Assert(
        nonlocal_vector.is_null() ||
          (last_action == VectorOperation::unknown) ||
          (last_action == VectorOperation::add),
        ExcMessage(
          "Cannot mix add and insert operations on a Tpetra vector "
          "with non-locally owned entries without calling compress() in between."));
 
      last_action = VectorOperation::add;
 
      // Create a non-owning pointer to the vector map. This way we avoid thread
      // safety problems with owning pointers.
      Teuchos::Ptr<
        const typename TpetraTypes::VectorType<Number, MemorySpace>::map_type>
        vector_map;
      {
#  ifndef HAVE_TEUCHOS_THREAD_SAFE
        // Make this part of the function thread safe
        std::scoped_lock lock(mutex);
#  endif
 
        vector_map = vector->getMap().ptr();
 
        // Make sure we have a view available to access the vector entries.
        if (!vector_1d_view)
          {
            auto vector_2d_view =
              vector->template getLocalView<Kokkos::HostSpace>(
                Tpetra::Access::ReadWriteStruct{});
 
            vector_1d_view = Kokkos::subview(vector_2d_view, Kokkos::ALL(), 0);
          }
 
        if (!nonlocal_vector_1d_view && !nonlocal_vector.is_null())
          {
            auto nonlocal_vector_2d_view =
              nonlocal_vector->template getLocalView<Kokkos::HostSpace>(
                Tpetra::Access::ReadWriteStruct{});
 
            nonlocal_vector_1d_view =
              Kokkos::subview(nonlocal_vector_2d_view, Kokkos::ALL(), 0);
          }
      }
 
      for (size_type i = 0; i < n_elements; ++i)
        {
          const size_type row = indices[i];
 
          // Check if the index is in the locally owned index set.
          // If so, we can write right into the locally owned
          // part of the vector.
          if (TrilinosWrappers::types::int_type local_row =
                vector_map->getLocalElement(row);
              local_row != Teuchos::OrdinalTraits<int>::invalid())
            {
              (*vector_1d_view)(local_row) += values[i];
 
              // Set the compressed state to false only if there is nonlocal
              // part in this distributed vector, otherwise it's always
              // compressed.
              if (nonlocal_vector.get() != nullptr)
                compressed = false;
            }
          else if (nonlocal_vector.get() == nullptr)
            {
#  ifndef HAVE_TEUCHOS_THREAD_SAFE
              // Make this part of the function thread safe
              std::scoped_lock lock(mutex);
#  endif
 
              // The element is not in the locally owned part, and
              // there is no specified nonlocal_vector, i.e. this is
              // a fully distributed vector with a nonlocal cache.
              // Add the element to the cache.
              nonlocal_cached_indices.push_back(row);
              nonlocal_cached_values.push_back(values[i]);
              compressed = false;
            }
          else
            {
#  ifndef HAVE_TEUCHOS_THREAD_SAFE
              // Make this part of the function thread safe
              std::scoped_lock lock(mutex);
#  endif
 
              // If the element was not in the locally owned part,
              // and we have a predetermined nonlocal buffer, we need
              // to figure out whether it is in the nonlocal
              // part. It better be:
              TrilinosWrappers::types::int_type nonlocal_row =
                nonlocal_vector->getMap()->getLocalElement(row);
 
#  if DEAL_II_TRILINOS_VERSION_GTE(14, 0, 0)
              Assert(nonlocal_row != Teuchos::OrdinalTraits<int>::invalid(),
                     ExcAccessToNonLocalElement(
                       row,
                       vector->getMap()->getLocalNumElements(),
                       vector->getMap()->getMinLocalIndex(),
                       vector->getMap()->getMaxLocalIndex()));
#  else
              Assert(nonlocal_row != Teuchos::OrdinalTraits<int>::invalid(),
                     ExcAccessToNonLocalElement(
                       row,
                       vector->getMap()->getNodeNumElements(),
                       vector->getMap()->getMinLocalIndex(),
                       vector->getMap()->getMaxLocalIndex()));
 
#  endif
 
              (*nonlocal_vector_1d_view)(nonlocal_row) += values[i];
              compressed = false;
            }
        }
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline void
    Vector<Number, MemorySpace>::set(const std::vector<size_type> &indices,
                                     const std::vector<Number>    &values)
    {
      AssertDimension(indices.size(), values.size());
      set(indices.size(), indices.data(), values.data());
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline void
    Vector<Number, MemorySpace>::set(const size_type  n_elements,
                                     const size_type *indices,
                                     const Number    *values)
    {
      // if we have ghost values, do not allow
      // writing to this vector at all.
      Assert(!has_ghost_elements(), ExcGhostsPresent());
 
      Assert(
        nonlocal_vector.is_null() ||
          (last_action == VectorOperation::unknown) ||
          (last_action == VectorOperation::insert),
        ExcMessage(
          "Cannot mix add and insert operations on a Tpetra vector "
          "with non-locally owned entries without calling compress() in between."));
 
      last_action = VectorOperation::insert;
 
      // Create a non-owning pointer to the vector map. This way we avoid thread
      // safety problems with owning pointers.
      Teuchos::Ptr<
        const typename TpetraTypes::VectorType<Number, MemorySpace>::map_type>
        vector_map;
      {
#  ifndef HAVE_TEUCHOS_THREAD_SAFE
        // Make this part of the function thread safe
        std::scoped_lock lock(mutex);
#  endif
 
        vector_map = vector->getMap().ptr();
 
        // Make sure we have a view available to access the vector entries.
        if (!vector_1d_view)
          {
            auto vector_2d_view =
              vector->template getLocalView<Kokkos::HostSpace>(
                Tpetra::Access::ReadWriteStruct{});
 
            vector_1d_view = Kokkos::subview(vector_2d_view, Kokkos::ALL(), 0);
          }
 
        if (!nonlocal_vector_1d_view && !nonlocal_vector.is_null())
          {
            auto nonlocal_vector_2d_view =
              nonlocal_vector->template getLocalView<Kokkos::HostSpace>(
                Tpetra::Access::ReadWriteStruct{});
 
            nonlocal_vector_1d_view =
              Kokkos::subview(nonlocal_vector_2d_view, Kokkos::ALL(), 0);
          }
      }
 
      for (size_type i = 0; i < n_elements; ++i)
        {
          const size_type row = indices[i];
 
          // Check if the index is in the locally owned index set.
          // If so, we can write right into the locally owned
          // part of the vector.
          if (TrilinosWrappers::types::int_type local_row =
                vector_map->getLocalElement(row);
              local_row != Teuchos::OrdinalTraits<int>::invalid())
            {
              (*vector_1d_view)(local_row) = values[i];
 
              // Set the compressed state to false only if there is nonlocal
              // part in this distributed vector, otherwise it's always
              // compressed.
              if (nonlocal_vector.get() != nullptr)
                compressed = false;
            }
          else if (nonlocal_vector.get() == nullptr)
            {
#  ifndef HAVE_TEUCHOS_THREAD_SAFE
              // Make this part of the function thread safe
              std::scoped_lock lock(mutex);
#  endif
 
              // The element is not in the locally owned part, and
              // there is no specified nonlocal_vector, i.e. this is
              // a fully distributed vector with a nonlocal cache.
              // Add the element to the cache.
              nonlocal_cached_indices.push_back(row);
              nonlocal_cached_values.push_back(values[i]);
              compressed = false;
            }
          else
            {
#  ifndef HAVE_TEUCHOS_THREAD_SAFE
              // Make this part of the function thread safe
              std::scoped_lock lock(mutex);
#  endif
 
              // If the element was not in the locally owned part,
              // and we have a predetermined nonlocal buffer, we need
              // to figure out whether it is in the nonlocal
              // part. It better be:
              TrilinosWrappers::types::int_type nonlocal_row =
                nonlocal_vector->getMap()->getLocalElement(row);
 
#  if DEAL_II_TRILINOS_VERSION_GTE(14, 0, 0)
              Assert(nonlocal_row != Teuchos::OrdinalTraits<int>::invalid(),
                     ExcAccessToNonLocalElement(
                       row,
                       vector->getMap()->getLocalNumElements(),
                       vector->getMap()->getMinLocalIndex(),
                       vector->getMap()->getMaxLocalIndex()));
#  else
              Assert(nonlocal_row != Teuchos::OrdinalTraits<int>::invalid(),
                     ExcAccessToNonLocalElement(
                       row,
                       vector->getMap()->getNodeNumElements(),
                       vector->getMap()->getMinLocalIndex(),
                       vector->getMap()->getMaxLocalIndex()));
 
#  endif
 
              (*nonlocal_vector_1d_view)(nonlocal_row) = values[i];
              compressed                               = false;
            }
        }
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline internal::VectorReference<Number, MemorySpace>
    Vector<Number, MemorySpace>::operator()(const size_type index)
    {
      return internal::VectorReference(*this, index);
    }
 
    template <typename Number, typename MemorySpace>
    inline internal::VectorReference<Number, MemorySpace>
    Vector<Number, MemorySpace>::operator[](const size_type index)
    {
      return operator()(index);
    }
 
    template <typename Number, typename MemorySpace>
    inline Number
    Vector<Number, MemorySpace>::operator[](const size_type index) const
    {
      return operator()(index);
    }
 
#  ifndef DOXYGEN
 
    // VectorReference
    namespace internal
    {
      template <typename Number, typename MemorySpace>
      inline VectorReference<Number, MemorySpace>::VectorReference(
        Vector<Number, MemorySpace> &vector,
        const size_type              index)
        : vector(vector)
        , index(index)
      {}
 
 
 
      template <typename Number, typename MemorySpace>
      inline const VectorReference<Number, MemorySpace> &
      VectorReference<Number, MemorySpace>::operator=(
        const VectorReference<Number, MemorySpace> &r) const
      {
        // as explained in the class
        // documentation, this is not the copy
        // operator. so simply pass on to the
        // "correct" assignment operator
        *this = static_cast<Number>(r);
 
        return *this;
      }
 
 
 
      template <typename Number, typename MemorySpace>
      inline VectorReference<Number, MemorySpace> &
      VectorReference<Number, MemorySpace>::operator=(
        const VectorReference<Number, MemorySpace> &r)
      {
        // as above
        *this = static_cast<Number>(r);
 
        return *this;
      }
 
 
 
      template <typename Number, typename MemorySpace>
      inline const VectorReference<Number, MemorySpace> &
      VectorReference<Number, MemorySpace>::operator=(const Number &value) const
      {
        Assert(!vector.has_ghost_elements(), ExcGhostsPresent());
 
        vector.set(1, &index, &value);
 
        return *this;
      }
 
 
 
      template <typename Number, typename MemorySpace>
      inline const VectorReference<Number, MemorySpace> &
      VectorReference<Number, MemorySpace>::operator+=(
        const Number &value) const
      {
        Assert(!vector.has_ghost_elements(), ExcGhostsPresent());
 
        vector.add(1, &index, &value);
 
        return *this;
      }
 
 
 
      template <typename Number, typename MemorySpace>
      inline const VectorReference<Number, MemorySpace> &
      VectorReference<Number, MemorySpace>::operator-=(
        const Number &value) const
      {
        Assert(!vector.has_ghost_elements(), ExcGhostsPresent());
 
        const Number new_value = -value;
        vector.add(1, &index, &new_value);
 
        return *this;
      }
 
 
 
      template <typename Number, typename MemorySpace>
      inline const VectorReference<Number, MemorySpace> &
      VectorReference<Number, MemorySpace>::operator*=(
        const Number &value) const
      {
        Assert(!vector.has_ghost_elements(), ExcGhostsPresent());
 
        const Number new_value = static_cast<Number>(*this) * value;
        vector.set(1, &index, &new_value);
 
        return *this;
      }
 
 
 
      template <typename Number, typename MemorySpace>
      inline const VectorReference<Number, MemorySpace> &
      VectorReference<Number, MemorySpace>::operator/=(
        const Number &value) const
      {
        Assert(!vector.has_ghost_elements(), ExcGhostsPresent());
 
        const Number new_value = static_cast<Number>(*this) / value;
        vector.set(1, &index, &new_value);
 
        return *this;
      }
    } // namespace internal
 
#  endif /* DOXYGEN */
 
  } // namespace TpetraWrappers

 
} // namespace LinearAlgebra
 
 
 
namespace internal
{
  namespace LinearOperatorImplementation
  {
    template <typename>
    class ReinitHelper;

    template <typename Number, typename MemorySpace>
    class ReinitHelper<
      LinearAlgebra::TpetraWrappers::Vector<Number, MemorySpace>>
    {
    public:
      template <typename Matrix>
      static void
      reinit_range_vector(
        const Matrix                                               &matrix,
        LinearAlgebra::TpetraWrappers::Vector<Number, MemorySpace> &v,
        bool omit_zeroing_entries)
      {
        v.reinit(matrix.locally_owned_range_indices(),
                 matrix.get_mpi_communicator(),
                 omit_zeroing_entries);
      }
 
      template <typename Matrix>
      static void
      reinit_domain_vector(
        const Matrix                                               &matrix,
        LinearAlgebra::TpetraWrappers::Vector<Number, MemorySpace> &v,
        bool omit_zeroing_entries)
      {
        v.reinit(matrix.locally_owned_domain_indices(),
                 matrix.get_mpi_communicator(),
                 omit_zeroing_entries);
      }
    };
 
  } // namespace LinearOperatorImplementation
} /* namespace internal */

template <typename Number, typename MemorySpace>
struct is_serial_vector<
  LinearAlgebra::TpetraWrappers::Vector<Number, MemorySpace>> : std::false_type
{};
 
#endif
 
DEAL_II_NAMESPACE_CLOSE
 
#endif