trilinos_tpetra_sparse_matrix.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) 2024 - 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_sparse_matrix_h
#define dealii_trilinos_tpetra_sparse_matrix_h
 
#include <deal.II/base/config.h>
 
#ifdef DEAL_II_TRILINOS_WITH_TPETRA
 
#  include <deal.II/base/enable_observer_pointer.h>
#  include <deal.II/base/index_set.h>
#  include <deal.II/base/trilinos_utilities.h>
 
#  include <deal.II/lac/sparse_matrix.h>
#  include <deal.II/lac/sparsity_pattern.h>
#  include <deal.II/lac/trilinos_tpetra_sparsity_pattern.h>
#  include <deal.II/lac/trilinos_tpetra_vector.h>
#  include <deal.II/lac/vector.h>
 
// Tpetra includes
#  include <Tpetra_Core.hpp>
#  include <Tpetra_CrsMatrix.hpp>
 
#  include <type_traits>
 
#endif // DEAL_II_TRILINOS_WITH_TPETRA
 
DEAL_II_NAMESPACE_OPEN
 
#ifdef DEAL_II_TRILINOS_WITH_TPETRA
// forward declarations
 
template <typename Number>
class SparseMatrix;
 
#  ifndef DOXYGEN
template <typename MatrixType>
class BlockMatrixBase;
 
namespace LinearAlgebra
{
  namespace TpetraWrappers
  {
    template <typename MemorySpace>
    class SparsityPattern;
 
    namespace SparseMatrixIterators
    {
      template <typename Number, typename MemorySpace, bool Constness>
      class Iterator;
    }
  } // namespace TpetraWrappers
} // namespace LinearAlgebra
#  endif
 
namespace LinearAlgebra
{
 
  namespace TpetraWrappers
  {
    template <typename Number, typename MemorySpace = ::MemorySpace::Host>
    class SparseMatrix : public EnableObserverPointer
    {
    public:
      using size_type = ::types::global_dof_index;
 
      DeclException1(ExcAccessToNonlocalRow,
                     std::size_t,
                     << "You tried to access row " << arg1
                     << " of a non-contiguous locally owned row set."
                     << " The row " << arg1
                     << " is not stored locally and can't be accessed.");
 
      struct Traits
      {
        static const bool zero_addition_can_be_elided = true;
      };
 
      using iterator =
        SparseMatrixIterators::Iterator<Number, MemorySpace, false>;
 
      using const_iterator =
        SparseMatrixIterators::Iterator<Number, MemorySpace, true>;
 
      using value_type = Number;
 
      SparseMatrix();
 
      SparseMatrix(const SparsityPattern<MemorySpace> &sparsity_pattern);
 
      SparseMatrix(const size_type    m,
                   const size_type    n,
                   const unsigned int n_max_entries_per_row);
 
      SparseMatrix(const size_type                  m,
                   const size_type                  n,
                   const std::vector<unsigned int> &n_entries_per_row);
 
      SparseMatrix(SparseMatrix<Number, MemorySpace> &&other) noexcept;
 
      SparseMatrix(const SparseMatrix<Number, MemorySpace> &) = delete;
 
      SparseMatrix<Number, MemorySpace> &
      operator=(const SparseMatrix<Number, MemorySpace> &) = delete;
 
      SparseMatrix<Number, MemorySpace> &
      operator=(SparseMatrix<Number, MemorySpace> &&other) noexcept;
 
      virtual ~SparseMatrix() override = default;
 
      template <typename SparsityPatternType>
      void
      reinit(const SparsityPatternType &sparsity_pattern);
 
      void
      reinit(const SparsityPattern<MemorySpace> &sparsity_pattern);
 
      void
      reinit(const SparseMatrix<Number, MemorySpace> &matrix);
      SparseMatrix(const IndexSet    &parallel_partitioning,
                   const MPI_Comm     communicator,
                   const unsigned int n_max_entries_per_row);
 
      SparseMatrix(const IndexSet                  &parallel_partitioning,
                   const MPI_Comm                   communicator,
                   const std::vector<unsigned int> &n_entries_per_row);
 
      SparseMatrix(const IndexSet &row_parallel_partitioning,
                   const IndexSet &col_parallel_partitioning,
                   const MPI_Comm  communicator,
                   const size_type n_max_entries_per_row);
 
      SparseMatrix(const IndexSet                  &row_parallel_partitioning,
                   const IndexSet                  &col_parallel_partitioning,
                   const MPI_Comm                   communicator,
                   const std::vector<unsigned int> &n_entries_per_row);
 
      template <typename SparsityPatternType>
      std::enable_if_t<
        !std::is_same_v<SparsityPatternType, ::SparseMatrix<double>>>
      reinit(const IndexSet            &parallel_partitioning,
             const SparsityPatternType &sparsity_pattern,
             const MPI_Comm             communicator  = MPI_COMM_WORLD,
             const bool                 exchange_data = false);
 
      template <typename SparsityPatternType>
      std::enable_if_t<
        !std::is_same_v<SparsityPatternType, ::SparseMatrix<double>>>
      reinit(const IndexSet            &row_parallel_partitioning,
             const IndexSet            &col_parallel_partitioning,
             const SparsityPatternType &sparsity_pattern,
             const MPI_Comm             communicator  = MPI_COMM_WORLD,
             const bool                 exchange_data = false);
 
      void
      reinit(const IndexSet                     &row_parallel_partitioning,
             const IndexSet                     &col_parallel_partitioning,
             const ::SparseMatrix<Number> &dealii_sparse_matrix,
             const MPI_Comm                      communicator = MPI_COMM_WORLD,
             const double                        drop_tolerance    = 1e-13,
             const bool                          copy_values       = true,
             const ::SparsityPattern      *use_this_sparsity = nullptr);
 
      size_type
      m() const;
 
      size_type
      n() const;
 
 
      unsigned int
      local_size() const;
 
      std::pair<size_type, size_type>
      local_range() const;
 
      bool
      in_local_range(const size_type index) const;
 
      size_t
      n_nonzero_elements() const;
 
      unsigned int
      row_length(const size_type row) const;
 
      bool
      is_compressed() const;
 
      size_type
      memory_consumption() const;
 
      MPI_Comm
      get_mpi_communicator() const;
      SparseMatrix &
      operator=(const double d);
 
      SparseMatrix &
      operator*=(const Number factor);
 
      SparseMatrix &
      operator/=(const Number factor);
 
      void
      copy_from(const SparseMatrix<Number, MemorySpace> &source);
 
      void
      add(const size_type i, const size_type j, const Number value);
 
      void
      add(const size_type  row,
          const size_type  n_cols,
          const size_type *col_indices,
          const Number    *values,
          const bool       elide_zero_values      = true,
          const bool       col_indices_are_sorted = false);
 
      void
      add(const Number factor, const SparseMatrix<Number, MemorySpace> &matrix);
 
      void
      set(const size_type i, const size_type j, const Number value);
 
      void
      set(const std::vector<size_type> &indices,
          const FullMatrix<Number>     &full_matrix,
          const bool                    elide_zero_values = false);
 
      void
      set(const std::vector<size_type> &row_indices,
          const std::vector<size_type> &col_indices,
          const FullMatrix<Number>     &full_matrix,
          const bool                    elide_zero_values = false);
 
      void
      set(const size_type               row,
          const std::vector<size_type> &col_indices,
          const std::vector<Number>    &values,
          const bool                    elide_zero_values = false);
 
      template <typename OtherNumber>
      void
      set(const size_type    row,
          const size_type    n_cols,
          const size_type   *col_indices,
          const OtherNumber *values,
          const bool         elide_zero_values = false);
 
      void
      clear_row(const size_type row, const Number new_diag_value = 0);
 
      void
      clear_rows(const ArrayView<const size_type> &rows,
                 const Number                      new_diag_value = 0);
 
      void
      clear();
 
      Number
      operator()(const size_type i, const size_type j) const;
 
      Number
      el(const size_type i, const size_type j) const;
 
      Number
      diag_element(const size_type i) const;
 
      template <typename InputVectorType>
      void
      vmult(InputVectorType &dst, const InputVectorType &src) const;
 
      template <typename InputVectorType>
      void
      Tvmult(InputVectorType &dst, const InputVectorType &src) const;
 
      template <typename InputVectorType>
      void
      vmult_add(InputVectorType &dst, const InputVectorType &src) const;
 
      template <typename InputVectorType>
      void
      Tvmult_add(InputVectorType &dst, const InputVectorType &src) const;
 
      Number
      matrix_norm_square(const Vector<Number, MemorySpace> &v) const;
 
      Number
      matrix_scalar_product(const Vector<Number, MemorySpace> &u,
                            const Vector<Number, MemorySpace> &v) const;
 
      Number
      residual(Vector<Number, MemorySpace>       &dst,
               const Vector<Number, MemorySpace> &x,
               const Vector<Number, MemorySpace> &b) const;
 
      Number
      frobenius_norm() const;
 
      Number
      l1_norm() const;
 
      Number
      linfty_norm() const;
 
      void
      print(std::ostream &out,
            const bool    print_detailed_trilinos_information = false) const;
 
      void
      compress(VectorOperation::values operation);
 
      void
      resume_fill();
 
      const TpetraTypes::MatrixType<Number, MemorySpace> &
      trilinos_matrix() const;
 
      TpetraTypes::MatrixType<Number, MemorySpace> &
      trilinos_matrix();
 
      Teuchos::RCP<const TpetraTypes::MatrixType<Number, MemorySpace>>
      trilinos_rcp() const;
 
      Teuchos::RCP<TpetraTypes::MatrixType<Number, MemorySpace>>
      trilinos_rcp();
      IndexSet
      locally_owned_domain_indices() const;
 
      IndexSet
      locally_owned_range_indices() const;
 
      const_iterator
      begin() const;
 
      iterator
      begin();
 
      const_iterator
      end() const;
 
      iterator
      end();
 
      const_iterator
      begin(const size_type r) const;
 
      iterator
      begin(const size_type r);
 
      const_iterator
      end(const size_type r) const;
 
      iterator
      end(const size_type r);
 
      DeclException0(ExcMatrixNotCompressed);
 
      DeclExceptionMsg(
        ExcSourceEqualsDestination,
        "You are attempting an operation on two vectors that "
        "are the same object, but the operation requires that the "
        "two objects are in fact different.");
 
      /*
       * Exception
       */
      DeclExceptionMsg(ExcColMapMismatch,
                       "The column partitioning of a matrix does not match "
                       "the partitioning of a vector you are trying to "
                       "multiply it with. Are you multiplying the "
                       "matrix with a vector that has ghost elements?");
 
      /*
       * Exception
       */
      DeclExceptionMsg(ExcDomainMapMismatch,
                       "The row partitioning of a matrix does not match "
                       "the partitioning of a vector you are trying to "
                       "put the result of a matrix-vector product in. "
                       "Are you trying to put the product of the "
                       "matrix with a vector into a vector that has "
                       "ghost elements?");
 
      DeclException2(ExcInvalidIndex,
                     size_type,
                     size_type,
                     << "The entry with index <" << arg1 << ',' << arg2
                     << "> does not exist.");
 
      DeclException4(ExcAccessToNonLocalElement,
                     size_type,
                     size_type,
                     size_type,
                     size_type,
                     << "You tried to access element (" << arg1 << '/' << arg2
                     << ')'
                     << " of a distributed matrix, but only rows in range ["
                     << arg3 << ',' << arg4
                     << "] are stored locally and can be accessed.");
 
    private:
      Number
      element(const size_type i, const size_type j, const bool no_error) const;
 
      Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> column_space_map;
 
      Teuchos::RCP<TpetraTypes::MatrixType<Number, MemorySpace>> matrix;
 
      bool compressed;
 
      void
      prepare_add();
 
      void
      prepare_set();
 
      // To allow calling protected prepare_add() and prepare_set().
      friend class BlockMatrixBase<SparseMatrix<Number, MemorySpace>>;
    }; // class SparseMatrix
 
    namespace SparseMatrixIterators
    {
      DeclException0(ExcBeyondEndOfMatrix);
 
      template <typename Number, typename MemorySpace>
      class AccessorBase
      {
      public:
        using size_type = ::types::global_dof_index;
 
        AccessorBase(SparseMatrix<Number, MemorySpace> *matrix,
                     const size_type                    row,
                     const size_type                    index);
 
        size_type
        row() const;
 
        size_type
        index() const;
 
        size_type
        column() const;
 
      protected:
        mutable SparseMatrix<Number, MemorySpace> *matrix;
        size_type a_row;
 
        size_type a_index;
 
        void
        visit_present_row();
 
        std::shared_ptr<std::vector<::types::signed_global_dof_index>>
          colnum_cache;
 
        std::shared_ptr<std::vector<Number>> value_cache;
 
      private:
        friend class Iterator<Number, MemorySpace, false>;
        friend class Iterator<Number, MemorySpace, true>;
      };
 
      template <typename Number, typename MemorySpace, bool Constness>
      class Accessor : public AccessorBase<Number, MemorySpace>
      {
        Number
        value() const;
 
        Number &
        value();
      };
 
      template <typename Number, typename MemorySpace>
      class Accessor<Number, MemorySpace, true>
        : public AccessorBase<Number, MemorySpace>
      {
      public:
        using MatrixType = const SparseMatrix<Number, MemorySpace>;
 
        using size_type = typename AccessorBase<Number, MemorySpace>::size_type;
 
        Accessor(MatrixType     *matrix,
                 const size_type row,
                 const size_type index);
 
        template <bool Other>
        Accessor(const Accessor<Number, MemorySpace, Other> &a);
 
        Number
        value() const;
      };
 
      template <typename Number, typename MemorySpace>
      class Accessor<Number, MemorySpace, false>
        : public AccessorBase<Number, MemorySpace>
      {
        class Reference
        {
        public:
          Reference(const Accessor<Number, MemorySpace, false> &accessor);
 
          operator Number() const;
 
          const Reference &
          operator=(const Number n) const;
 
          const Reference &
          operator+=(const Number n) const;
 
          const Reference &
          operator-=(const Number n) const;
 
          const Reference &
          operator*=(const Number n) const;
 
          const Reference &
          operator/=(const Number n) const;
 
        private:
          Accessor &accessor;
        };
 
      public:
        using MatrixType = SparseMatrix<Number, MemorySpace>;
 
        using size_type = typename AccessorBase<Number, MemorySpace>::size_type;
 
        Accessor(MatrixType     *matrix,
                 const size_type row,
                 const size_type index);
 
        Reference
        value() const;
 
      private:
        // Make Reference object a friend.
        friend class Reference;
      };
 
      template <typename Number, typename MemorySpace, bool Constness>
      class Iterator
      {
      public:
        using size_type = ::types::global_dof_index;
 
        using difference_type = ::types::global_dof_index;
 
        using value_type = Number;
 
        using MatrixType =
          typename Accessor<Number, MemorySpace, Constness>::MatrixType;
 
        Iterator(MatrixType     *matrix,
                 const size_type row,
                 const size_type index);
 
        template <bool Other>
        Iterator(const Iterator<Number, MemorySpace, Other> &other);
 
        Iterator<Number, MemorySpace, Constness> &
        operator++();
 
        Iterator<Number, MemorySpace, Constness>
        operator++(int);
 
        const Accessor<Number, MemorySpace, Constness> &
        operator*() const;
 
        const Accessor<Number, MemorySpace, Constness> *
        operator->() const;
 
        template <bool OtherConstness>
        bool
        operator==(const Iterator<Number, MemorySpace, OtherConstness> &) const;
 
        template <bool OtherConstness>
        bool
        operator!=(const Iterator<Number, MemorySpace, OtherConstness> &) const;
 
        template <bool OtherConstness>
        bool
        operator<(const Iterator<Number, MemorySpace, OtherConstness> &) const;
 
        template <bool OtherConstness>
        bool
        operator>(const Iterator<Number, MemorySpace, OtherConstness> &) const;
 
        DeclException2(ExcInvalidIndexWithinRow,
                       size_type,
                       size_type,
                       << "Attempt to access element " << arg2 << " of row "
                       << arg1 << " which doesn't have that many elements.");
 
      private:
        Accessor<Number, MemorySpace, Constness> accessor;
 
        friend class Iterator<Number, MemorySpace, true>;
        friend class Iterator<Number, MemorySpace, false>;
      };
    } // namespace SparseMatrixIterators
 
  } // namespace TpetraWrappers
 
} // namespace LinearAlgebra
 
 
DEAL_II_NAMESPACE_CLOSE // Do not convert for module purposes
 
  namespace std
{
  template <typename Number, typename MemorySpace, bool Constness>
  struct iterator_traits<
    ::LinearAlgebra::TpetraWrappers::SparseMatrixIterators::
      Iterator<Number, MemorySpace, Constness>>
  {
    using iterator_category = forward_iterator_tag;
    using value_type =
      typename ::LinearAlgebra::TpetraWrappers::SparseMatrixIterators::
        Iterator<Number, MemorySpace, Constness>::value_type;
    using difference_type =
      typename ::LinearAlgebra::TpetraWrappers::SparseMatrixIterators::
        Iterator<Number, MemorySpace, Constness>::difference_type;
  };
} // namespace std
 
/* ------------------------- Inline functions ---------------------- */
 
 
DEAL_II_NAMESPACE_OPEN // Do not convert for module purposes
 
  namespace LinearAlgebra
{
  namespace TpetraWrappers
  {
    template <typename Number, typename MemorySpace>
    inline void
    SparseMatrix<Number, MemorySpace>::set(const size_type i,
                                           const size_type j,
                                           const Number    value)
    {
      set(i, 1, &j, &value, false);
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline void
    SparseMatrix<Number, MemorySpace>::add(const size_type i,
                                           const size_type j,
                                           const Number    value)
    {
      add(i, 1, &j, &value, false);
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline Number
    SparseMatrix<Number, MemorySpace>::residual(
      Vector<Number, MemorySpace>       &dst,
      const Vector<Number, MemorySpace> &x,
      const Vector<Number, MemorySpace> &b) const
    {
      vmult(dst, x);
      dst -= b;
      dst *= -1.;
 
      return dst.l2_norm();
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline typename SparseMatrix<Number, MemorySpace>::size_type
    SparseMatrix<Number, MemorySpace>::m() const
    {
      return matrix->getRowMap()->getGlobalNumElements();
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline typename SparseMatrix<Number, MemorySpace>::size_type
    SparseMatrix<Number, MemorySpace>::n() const
    {
      // If the matrix structure has not been fixed (i.e., we did not have a
      // sparsity pattern), it does not know about the number of columns, so we
      // must always take this from the additional column space map
      Assert(column_space_map.get() != nullptr, ExcInternalError());
      return column_space_map->getGlobalNumElements();
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline typename SparseMatrix<Number, MemorySpace>::size_type
    SparseMatrix<Number, MemorySpace>::memory_consumption() const
    {
      size_type static_memory =
        sizeof(*this) + sizeof(*matrix) + sizeof(matrix->getGraph().get());
      return ((sizeof(Number) + sizeof(TrilinosWrappers::types::int_type)) *
#  if DEAL_II_TRILINOS_VERSION_GTE(13, 4, 0)
                matrix->getLocalNumEntries() +
#  else
                matrix->getNodeNumEntries() +
#  endif
              sizeof(int) * local_size() + static_memory);
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline bool
    SparseMatrix<Number, MemorySpace>::is_compressed() const
    {
      return compressed;
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline typename SparseMatrix<Number, MemorySpace>::const_iterator
    SparseMatrix<Number, MemorySpace>::begin() const
    {
      return begin(0);
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline typename SparseMatrix<Number, MemorySpace>::const_iterator
    SparseMatrix<Number, MemorySpace>::end() const
    {
      return const_iterator(this, m(), 0);
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline typename SparseMatrix<Number, MemorySpace>::const_iterator
    SparseMatrix<Number, MemorySpace>::begin(const size_type r) const
    {
      AssertIndexRange(r, m());
      if (in_local_range(r) && (row_length(r) > 0))
        return const_iterator(this, r, 0);
      else
        return end(r);
    }
 
 
    template <typename Number, typename MemorySpace>
    inline typename SparseMatrix<Number, MemorySpace>::const_iterator
    SparseMatrix<Number, MemorySpace>::end(const size_type r) const
    {
      AssertIndexRange(r, m());
 
      // place the iterator on the first entry
      // past this line, or at the end of the
      // matrix
      for (size_type i = r + 1; i < m(); ++i)
        if (in_local_range(i) && (row_length(i) > 0))
          return const_iterator(this, i, 0);
 
      // if there is no such line, then take the
      // end iterator of the matrix
      return end();
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline typename SparseMatrix<Number, MemorySpace>::iterator
    SparseMatrix<Number, MemorySpace>::begin()
    {
      return begin(0);
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline typename SparseMatrix<Number, MemorySpace>::iterator
    SparseMatrix<Number, MemorySpace>::end()
    {
      return iterator(this, m(), 0);
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline typename SparseMatrix<Number, MemorySpace>::iterator
    SparseMatrix<Number, MemorySpace>::begin(const size_type r)
    {
      AssertIndexRange(r, m());
      if (in_local_range(r) && (row_length(r) > 0))
        return iterator(this, r, 0);
      else
        return end(r);
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline typename SparseMatrix<Number, MemorySpace>::iterator
    SparseMatrix<Number, MemorySpace>::end(const size_type r)
    {
      AssertIndexRange(r, m());
 
      // place the iterator on the first entry
      // past this line, or at the end of the
      // matrix
      for (size_type i = r + 1; i < m(); ++i)
        if (in_local_range(i) && (row_length(i) > 0))
          return iterator(this, i, 0);
 
      // if there is no such line, then take the
      // end iterator of the matrix
      return end();
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline bool
    SparseMatrix<Number, MemorySpace>::in_local_range(
      const size_type index) const
    {
      const size_type begin = matrix->getRowMap()->getMinGlobalIndex();
      const size_type end   = matrix->getRowMap()->getMaxGlobalIndex() + 1;
 
      return ((index >= begin) && (index < end));
    }
 
 
 
    template <typename Number, typename MemorySpace>
    unsigned int
    SparseMatrix<Number, MemorySpace>::row_length(const size_type row) const
    {
      auto n_entries = matrix->getNumEntriesInGlobalRow(row);
      Assert(n_entries !=
               Teuchos::OrdinalTraits<decltype(n_entries)>::invalid(),
             ExcAccessToNonlocalRow(row));
 
      return n_entries;
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline void
    SparseMatrix<Number, MemorySpace>::prepare_add()
    {
      // nothing to do here
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline void
    SparseMatrix<Number, MemorySpace>::prepare_set()
    {
      // nothing to do here
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline const TpetraTypes::MatrixType<Number, MemorySpace> &
    SparseMatrix<Number, MemorySpace>::trilinos_matrix() const
    {
      return *matrix;
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline TpetraTypes::MatrixType<Number, MemorySpace> &
    SparseMatrix<Number, MemorySpace>::trilinos_matrix()
    {
      return *matrix;
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline Teuchos::RCP<const TpetraTypes::MatrixType<Number, MemorySpace>>
    SparseMatrix<Number, MemorySpace>::trilinos_rcp() const
    {
      return matrix.getConst();
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline Teuchos::RCP<TpetraTypes::MatrixType<Number, MemorySpace>>
    SparseMatrix<Number, MemorySpace>::trilinos_rcp()
    {
      return matrix;
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline IndexSet
    SparseMatrix<Number, MemorySpace>::locally_owned_domain_indices() const
    {
      return IndexSet(matrix->getDomainMap());
    }
 
 
 
    template <typename Number, typename MemorySpace>
    inline IndexSet
    SparseMatrix<Number, MemorySpace>::locally_owned_range_indices() const
    {
      return IndexSet(matrix->getRangeMap());
    }
 
 
    namespace SparseMatrixIterators
    {
      template <typename Number, typename MemorySpace>
      inline AccessorBase<Number, MemorySpace>::AccessorBase(
        SparseMatrix<Number, MemorySpace> *matrix,
        size_type                          row,
        size_type                          index)
        : matrix(matrix)
        , a_row(row)
        , a_index(index)
      {
        visit_present_row();
      }
 
 
      template <typename Number, typename MemorySpace>
      inline typename AccessorBase<Number, MemorySpace>::size_type
      AccessorBase<Number, MemorySpace>::row() const
      {
        Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
        return a_row;
      }
 
 
      template <typename Number, typename MemorySpace>
      inline typename AccessorBase<Number, MemorySpace>::size_type
      AccessorBase<Number, MemorySpace>::column() const
      {
        Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
        return (*colnum_cache)[a_index];
      }
 
 
      template <typename Number, typename MemorySpace>
      inline typename AccessorBase<Number, MemorySpace>::size_type
      AccessorBase<Number, MemorySpace>::index() const
      {
        Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
        return a_index;
      }
 
 
 
      template <typename Number, typename MemorySpace>
      void
      AccessorBase<Number, MemorySpace>::visit_present_row()
      {
        // if we are asked to visit the past-the-end line, then simply
        // release all our caches and go on with life.
        //
        // do the same if the row we're supposed to visit is not locally
        // owned. this is simply going to make non-locally owned rows
        // look like they're empty
        if ((this->a_row == matrix->m()) ||
            (matrix->in_local_range(this->a_row) == false))
          {
            colnum_cache.reset();
            value_cache.reset();
 
            return;
          }
 
        // get a representation of the present row
        size_t                            ncols;
        TrilinosWrappers::types::int_type colnums =
          matrix->row_length(this->a_row);
        if (value_cache.get() == nullptr)
          {
            value_cache  = std::make_shared<std::vector<Number>>(colnums);
            colnum_cache = std::make_shared<
              std::vector<::types::signed_global_dof_index>>(colnums);
          }
        else
          {
            value_cache->resize(colnums);
            colnum_cache->resize(colnums);
          }
 
        typename TpetraTypes::MatrixType<Number, MemorySpace>::
          nonconst_global_inds_host_view_type col_indices(colnum_cache->data(),
                                                          colnums);
        typename TpetraTypes::MatrixType<Number, MemorySpace>::
          nonconst_values_host_view_type values(value_cache->data(), colnums);
 
        matrix->trilinos_matrix().getGlobalRowCopy(this->a_row,
                                                   col_indices,
                                                   values,
                                                   ncols);
 
        AssertDimension(ncols, colnums);
 
        // copy it into our caches if the
        // line isn't empty. if it is, then
        // we've done something wrong, since
        // we shouldn't have initialized an
        // iterator for an empty line (what
        // would it point to?)
      }
 
 
 
      template <typename Number, typename MemorySpace>
      inline Accessor<Number, MemorySpace, true>::Accessor(
        MatrixType     *matrix,
        const size_type row,
        const size_type index)
        : AccessorBase<Number, MemorySpace>(
            const_cast<SparseMatrix<Number, MemorySpace> *>(matrix),
            row,
            index)
      {}
 
 
 
      template <typename Number, typename MemorySpace>
      template <bool Other>
      inline Accessor<Number, MemorySpace, true>::Accessor(
        const Accessor<Number, MemorySpace, Other> &other)
        : AccessorBase<Number, MemorySpace>(other)
      {}
 
 
 
      template <typename Number, typename MemorySpace>
      inline Number
      Accessor<Number, MemorySpace, true>::value() const
      {
        Assert((AccessorBase<Number, MemorySpace>::a_row <
                AccessorBase<Number, MemorySpace>::matrix->m()),
               ExcBeyondEndOfMatrix());
        return (*AccessorBase<Number, MemorySpace>::value_cache)
          [AccessorBase<Number, MemorySpace>::a_index];
      }
 
 
 
      template <typename Number, typename MemorySpace>
      inline Accessor<Number, MemorySpace, false>::Reference::Reference(
        const Accessor<Number, MemorySpace, false> &acc)
        : accessor(const_cast<Accessor<Number, MemorySpace, false> &>(acc))
      {}
 
 
 
      template <typename Number, typename MemorySpace>
      inline Accessor<Number, MemorySpace, false>::Reference::operator Number()
        const
      {
        return (*accessor.value_cache)[accessor.a_index];
      }
 
 
 
      template <typename Number, typename MemorySpace>
      inline const typename Accessor<Number, MemorySpace, false>::Reference &
      Accessor<Number, MemorySpace, false>::Reference::operator=(
        const Number n) const
      {
        (*accessor.value_cache)[accessor.a_index] = n;
        accessor.matrix->set(accessor.row(),
                             accessor.column(),
                             static_cast<Number>(*this));
        return *this;
      }
 
 
 
      template <typename Number, typename MemorySpace>
      inline const typename Accessor<Number, MemorySpace, false>::Reference &
      Accessor<Number, MemorySpace, false>::Reference::operator+=(
        const Number n) const
      {
        (*accessor.value_cache)[accessor.a_index] += n;
        accessor.matrix->set(accessor.row(),
                             accessor.column(),
                             static_cast<Number>(*this));
        return *this;
      }
 
 
 
      template <typename Number, typename MemorySpace>
      inline const typename Accessor<Number, MemorySpace, false>::Reference &
      Accessor<Number, MemorySpace, false>::Reference::operator-=(
        const Number n) const
      {
        (*accessor.value_cache)[accessor.a_index] -= n;
        accessor.matrix->set(accessor.row(),
                             accessor.column(),
                             static_cast<Number>(*this));
        return *this;
      }
 
 
      template <typename Number, typename MemorySpace>
      inline const typename Accessor<Number, MemorySpace, false>::Reference &
      Accessor<Number, MemorySpace, false>::Reference::operator*=(
        const Number n) const
      {
        (*accessor.value_cache)[accessor.a_index] *= n;
        accessor.matrix->set(accessor.row(),
                             accessor.column(),
                             static_cast<Number>(*this));
        return *this;
      }
 
 
      template <typename Number, typename MemorySpace>
      inline const typename Accessor<Number, MemorySpace, false>::Reference &
      Accessor<Number, MemorySpace, false>::Reference::operator/=(
        const Number n) const
      {
        (*accessor.value_cache)[accessor.a_index] /= n;
        accessor.matrix->set(accessor.row(),
                             accessor.column(),
                             static_cast<Number>(*this));
        return *this;
      }
 
 
      template <typename Number, typename MemorySpace>
      inline Accessor<Number, MemorySpace, false>::Accessor(
        MatrixType     *matrix,
        const size_type row,
        const size_type index)
        : AccessorBase<Number, MemorySpace>(matrix, row, index)
      {}
 
 
      template <typename Number, typename MemorySpace>
      inline typename Accessor<Number, MemorySpace, false>::Reference
      Accessor<Number, MemorySpace, false>::value() const
      {
        Assert((AccessorBase<Number, MemorySpace>::a_row <
                AccessorBase<Number, MemorySpace>::matrix->m()),
               ExcBeyondEndOfMatrix());
        return {*this};
      }
 
 
 
      template <typename Number, typename MemorySpace, bool Constness>
      inline Iterator<Number, MemorySpace, Constness>::Iterator(
        MatrixType     *matrix,
        const size_type row,
        const size_type index)
        : accessor(matrix, row, index)
      {}
 
 
      template <typename Number, typename MemorySpace, bool Constness>
      template <bool Other>
      inline Iterator<Number, MemorySpace, Constness>::Iterator(
        const Iterator<Number, MemorySpace, Other> &other)
        : accessor(other.accessor)
      {}
 
 
 
      template <typename Number, typename MemorySpace, bool Constness>
      inline Iterator<Number, MemorySpace, Constness> &
      Iterator<Number, MemorySpace, Constness>::operator++()
      {
        Assert(accessor.a_row < accessor.matrix->m(), ExcIteratorPastEnd());
 
        ++accessor.a_index;
 
        // If at end of line: do one
        // step, then cycle until we
        // find a row with a nonzero
        // number of entries.
        if (accessor.a_index >= accessor.colnum_cache->size())
          {
            accessor.a_index = 0;
            ++accessor.a_row;
 
            while (
              (accessor.a_row < accessor.matrix->m()) &&
              ((accessor.matrix->in_local_range(accessor.a_row) == false) ||
               (accessor.matrix->row_length(accessor.a_row) == 0)))
              ++accessor.a_row;
 
            accessor.visit_present_row();
          }
        return *this;
      }
 
 
      template <typename Number, typename MemorySpace, bool Constness>
      inline Iterator<Number, MemorySpace, Constness>
      Iterator<Number, MemorySpace, Constness>::operator++(int)
      {
        const Iterator<Number, MemorySpace, Constness> old_state = *this;
        ++(*this);
        return old_state;
      }
 
 
 
      template <typename Number, typename MemorySpace, bool Constness>
      inline const Accessor<Number, MemorySpace, Constness> &
      Iterator<Number, MemorySpace, Constness>::operator*() const
      {
        return accessor;
      }
 
 
 
      template <typename Number, typename MemorySpace, bool Constness>
      inline const Accessor<Number, MemorySpace, Constness> *
      Iterator<Number, MemorySpace, Constness>::operator->() const
      {
        return &accessor;
      }
 
 
 
      template <typename Number, typename MemorySpace, bool Constness>
      template <bool OtherConstness>
      inline bool
      Iterator<Number, MemorySpace, Constness>::operator==(
        const Iterator<Number, MemorySpace, OtherConstness> &other) const
      {
        return (accessor.a_row == other.accessor.a_row &&
                accessor.a_index == other.accessor.a_index);
      }
 
 
 
      template <typename Number, typename MemorySpace, bool Constness>
      template <bool OtherConstness>
      inline bool
      Iterator<Number, MemorySpace, Constness>::operator!=(
        const Iterator<Number, MemorySpace, OtherConstness> &other) const
      {
        return !(*this == other);
      }
 
 
 
      template <typename Number, typename MemorySpace, bool Constness>
      template <bool OtherConstness>
      inline bool
      Iterator<Number, MemorySpace, Constness>::operator<(
        const Iterator<Number, MemorySpace, OtherConstness> &other) const
      {
        return (accessor.row() < other.accessor.row() ||
                (accessor.row() == other.accessor.row() &&
                 accessor.index() < other.accessor.index()));
      }
 
 
      template <typename Number, typename MemorySpace, bool Constness>
      template <bool OtherConstness>
      inline bool
      Iterator<Number, MemorySpace, Constness>::operator>(
        const Iterator<Number, MemorySpace, OtherConstness> &other) const
      {
        return (other < *this);
      }
 
    } // namespace SparseMatrixIterators
  }   // namespace TpetraWrappers
 
} // namespace LinearAlgebra
 
#endif // DEAL_II_TRILINOS_WITH_TPETRA
 
DEAL_II_NAMESPACE_CLOSE
 
#endif // dealii_trilinos_tpetra_sparse_matrix_h