trilinos_tpetra_vector.h

Go to the documentation of this file.

// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 2018 - 2024 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Part of the source code is dual licensed under Apache-2.0 WITH
// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
// governing the source code and code contributions can be found in
// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
//
// ------------------------------------------------------------------------
 
#ifndef dealii_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 <memory>
#  include <optional>
 
DEAL_II_NAMESPACE_OPEN
 
template <typename Number>
struct is_tpetra_type : std::false_type
{};
 
#  ifdef HAVE_TPETRA_INST_FLOAT
template <>
struct is_tpetra_type<float> : std::true_type
{};
#  endif
 
#  ifdef HAVE_TPETRA_INST_DOUBLE
template <>
struct is_tpetra_type<double> : std::true_type
{};
#  endif
 
#  ifdef DEAL_II_WITH_COMPLEX_VALUES
#    ifdef HAVE_TPETRA_INST_COMPLEX_FLOAT
template <>
struct is_tpetra_type<std::complex<float>> : std::true_type
{};
#    endif
 
#    ifdef HAVE_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);
 
      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);
 
      DEAL_II_DEPRECATED
      void
      import_elements(
        const ReadWriteVector<Number> &V,
        VectorOperation::values        operation,
        const std::shared_ptr<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);
 
      DEAL_II_DEPRECATED
      void
      import(const ReadWriteVector<Number> &V,
             VectorOperation::values        operation,
             std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
               communication_pattern = {})
      {
        import_elements(V, operation, communication_pattern);
      }
 
      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 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
      set(const size_type  n_elements,
          const size_type *indices,
          const 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;
 
      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);
 
      const TpetraTypes::VectorType<Number, MemorySpace> &
      trilinos_vector() const;
 
      TpetraTypes::VectorType<Number, MemorySpace> &
      trilinos_vector();
 
      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;
 
      Teuchos::RCP<TpetraTypes::VectorType<Number, MemorySpace>> vector;
 
      Teuchos::RCP<TpetraTypes::VectorType<Number, MemorySpace>>
        nonlocal_vector;
 
      ::IndexSet source_stored_elements;
 
      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
    {
      vector.swap(v.vector);
    }
 
 
    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());
 
#  if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
      auto vector_2d_local = vector->template getLocalView<Kokkos::HostSpace>(
        Tpetra::Access::ReadWrite);
#  else
      vector->template sync<Kokkos::HostSpace>();
 
      auto vector_2d_local = vector->template getLocalView<Kokkos::HostSpace>();
#  endif
 
      // Having extracted a view into the multivectors above, now also
      // extract a view into the one vector we actually store. We can
      // do this right away for the locally owned part. We defer creating
      // the view into the nonlocal part to when we know that we actually
      // need it; this also makes sure that we correctly deal with the
      // case where we do not actually store a nonlocal part.
      auto vector_1d_local = Kokkos::subview(vector_2d_local, Kokkos::ALL(), 0);
      using ViewType1d     = decltype(vector_1d_local);
      std::optional<ViewType1d> vector_1d_nonlocal;
 
#  if !DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
      // Mark vector as to-be-modified. We may do the same with
      // the nonlocal part too if we end up writing into it.
      vector->template modify<Kokkos::HostSpace>();
#  endif
 
      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->getMap()->getLocalElement(row);
              local_row != Teuchos::OrdinalTraits<int>::invalid())
            {
              vector_1d_local(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 the element was not in the locally owned part,
              // we need to figure out whether it is in the nonlocal
              // part. It better be:
              Assert(nonlocal_vector.get() != nullptr, ExcInternalError());
              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
 
              // Having asserted that it is, write into the nonlocal part.
              // To do so, we first need to make sure that we have a view
              // of the nonlocal part of the vectors, since we have
              // deferred creating this view previously:
              if (!vector_1d_nonlocal)
                {
#  if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
                  auto vector_2d_nonlocal =
                    nonlocal_vector->template getLocalView<Kokkos::HostSpace>(
                      Tpetra::Access::ReadWrite);
#  else
                  auto vector_2d_nonlocal =
                    nonlocal_vector->template getLocalView<Kokkos::HostSpace>();
#  endif
 
                  vector_1d_nonlocal =
                    Kokkos::subview(vector_2d_nonlocal, Kokkos::ALL(), 0);
 
#  if !DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
                  // Mark the nonlocal vector as to-be-modified as well.
                  nonlocal_vector->template modify<Kokkos::HostSpace>();
#  endif
                }
              (*vector_1d_nonlocal)(nonlocal_row) += values[i];
              compressed = false;
            }
        }
 
#  if !DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
      vector->template sync<
        typename Tpetra::Vector<Number, int, types::signed_global_dof_index>::
          device_type::memory_space>();
 
      // If we have created a view to the nonlocal part, then we have also
      // written into it. Flush these modifications.
      if (vector_1d_nonlocal)
        nonlocal_vector->template sync<
          typename Tpetra::Vector<Number, int, types::signed_global_dof_index>::
            device_type::memory_space>();
#  endif
    }
 
 
 
    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());
 
#  if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
      auto vector_2d_local = vector->template getLocalView<Kokkos::HostSpace>(
        Tpetra::Access::ReadWrite);
#  else
      vector->template sync<Kokkos::HostSpace>();
 
      auto vector_2d_local = vector->template getLocalView<Kokkos::HostSpace>();
#  endif
 
      // Having extracted a view into the multivectors above, now also
      // extract a view into the one vector we actually store. We can
      // do this right away for the locally owned part. We defer creating
      // the view into the nonlocal part to when we know that we actually
      // need it; this also makes sure that we correctly deal with the
      // case where we do not actually store a nonlocal part.
      auto vector_1d_local = Kokkos::subview(vector_2d_local, Kokkos::ALL(), 0);
      using ViewType1d     = decltype(vector_1d_local);
      std::optional<ViewType1d> vector_1d_nonlocal;
 
#  if !DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
      // Mark vector as to-be-modified. We may do the same with
      // the nonlocal part too if we end up writing into it.
      vector->template modify<Kokkos::HostSpace>();
#  endif
 
      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->getMap()->getLocalElement(row);
              local_row != Teuchos::OrdinalTraits<int>::invalid())
            {
              vector_1d_local(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 the element was not in the locally owned part,
              // we need to figure out whether it is in the nonlocal
              // part. It better be:
              Assert(nonlocal_vector.get() != nullptr, ExcInternalError());
              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
 
              // Having asserted that it is, write into the nonlocal part.
              // To do so, we first need to make sure that we have a view
              // of the nonlocal part of the vectors, since we have
              // deferred creating this view previously:
              if (!vector_1d_nonlocal)
                {
#  if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
                  auto vector_2d_nonlocal =
                    nonlocal_vector->template getLocalView<Kokkos::HostSpace>(
                      Tpetra::Access::ReadWrite);
#  else
                  auto vector_2d_nonlocal =
                    nonlocal_vector->template getLocalView<Kokkos::HostSpace>();
#  endif
 
                  vector_1d_nonlocal =
                    Kokkos::subview(vector_2d_nonlocal, Kokkos::ALL(), 0);
 
#  if !DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
                  // Mark the nonlocal vector as to-be-modified as well.
                  nonlocal_vector->template modify<Kokkos::HostSpace>();
#  endif
                }
              (*vector_1d_nonlocal)(nonlocal_row) = values[i];
              compressed                          = false;
            }
        }
 
#  if !DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
      vector->template sync<
        typename Tpetra::Vector<Number, int, types::signed_global_dof_index>::
          device_type::memory_space>();
 
      // If we have created a view to the nonlocal part, then we have also
      // written into it. Flush these modifications.
      if (vector_1d_nonlocal)
        nonlocal_vector->template sync<
          typename Tpetra::Vector<Number, int, types::signed_global_dof_index>::
            device_type::memory_space>();
#  endif
    }
 
 
 
    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
      {
        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
      {
        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
      {
        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
      {
        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
      {
        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
{};
 
DEAL_II_NAMESPACE_CLOSE
 
#endif
 
#endif