trilinos_sparse_matrix.h

Go to the documentation of this file.

// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 2008 - 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_sparse_matrix_h
#  define dealii_trilinos_sparse_matrix_h
 
 
#  include <deal.II/base/config.h>
 
#  ifdef DEAL_II_WITH_TRILINOS
 
#    include <deal.II/base/enable_observer_pointer.h>
#    include <deal.II/base/index_set.h>
#    include <deal.II/base/mpi_stub.h>
 
#    include <deal.II/lac/exceptions.h>
#    include <deal.II/lac/full_matrix.h>
#    include <deal.II/lac/trilinos_epetra_vector.h>
#    include <deal.II/lac/trilinos_index_access.h>
#    include <deal.II/lac/trilinos_tpetra_vector.h>
#    include <deal.II/lac/trilinos_vector.h>
#    include <deal.II/lac/vector_memory.h>
#    include <deal.II/lac/vector_operation.h>
 
#    include <Epetra_Comm.h>
#    include <Epetra_CrsGraph.h>
#    include <Epetra_Export.h>
#    include <Epetra_FECrsMatrix.h>
#    include <Epetra_Map.h>
#    include <Epetra_MpiComm.h>
#    include <Epetra_MultiVector.h>
#    include <Epetra_Operator.h>
 
#    include <cmath>
#    include <iterator>
#    include <memory>
#    include <type_traits>
#    include <vector>
 
DEAL_II_NAMESPACE_OPEN
 
// forward declarations
#    ifndef DOXYGEN
template <typename MatrixType>
class BlockMatrixBase;
 
template <typename number>
class SparseMatrix;
class SparsityPattern;
class DynamicSparsityPattern;
 
namespace TrilinosWrappers
{
  class SparseMatrix;
  class SparsityPattern;
 
  namespace SparseMatrixIterators
  {
    template <bool Constness>
    class Iterator;
  }
} // namespace TrilinosWrappers
#    endif
 
namespace TrilinosWrappers
{
  namespace SparseMatrixIterators
  {
    DeclException0(ExcBeyondEndOfMatrix);
 
    DeclException3(ExcAccessToNonlocalRow,
                   std::size_t,
                   std::size_t,
                   std::size_t,
                   << "You tried to access row " << arg1
                   << " of a distributed sparsity pattern, "
                   << " but only rows " << arg2 << " through " << arg3
                   << " are stored locally and can be accessed.");
 
    class AccessorBase
    {
    public:
      using size_type = ::types::global_dof_index;
 
      AccessorBase(SparseMatrix   *matrix,
                   const size_type row,
                   const size_type index);
 
      size_type
      row() const;
 
      size_type
      index() const;
 
      size_type
      column() const;
 
    protected:
      mutable SparseMatrix *matrix;
      size_type a_row;
 
      size_type a_index;
 
      void
      visit_present_row();
 
      std::shared_ptr<std::vector<size_type>> colnum_cache;
 
      std::shared_ptr<std::vector<TrilinosScalar>> value_cache;
    };
 
    template <bool Constess>
    class Accessor : public AccessorBase
    {
      TrilinosScalar
      value() const;
 
      TrilinosScalar &
      value();
    };
 
 
 
    template <>
    class Accessor<true> : public AccessorBase
    {
    public:
      using MatrixType = const SparseMatrix;
 
      Accessor(MatrixType *matrix, const size_type row, const size_type index);
 
      template <bool Other>
      Accessor(const Accessor<Other> &a);
 
      TrilinosScalar
      value() const;
 
    private:
      // Make iterator class a friend.
      template <bool>
      friend class Iterator;
    };
 
    template <>
    class Accessor<false> : public AccessorBase
    {
      class Reference
      {
      public:
        Reference(const Accessor<false> &accessor);
 
        operator TrilinosScalar() const;
 
        const Reference &
        operator=(const TrilinosScalar n) const;
 
        const Reference &
        operator+=(const TrilinosScalar n) const;
 
        const Reference &
        operator-=(const TrilinosScalar n) const;
 
        const Reference &
        operator*=(const TrilinosScalar n) const;
 
        const Reference &
        operator/=(const TrilinosScalar n) const;
 
      private:
        Accessor &accessor;
      };
 
    public:
      using MatrixType = SparseMatrix;
 
      Accessor(MatrixType *matrix, const size_type row, const size_type index);
 
      Reference
      value() const;
 
    private:
      // Make iterator class a friend.
      template <bool>
      friend class Iterator;
 
      // Make Reference object a friend.
      friend class Reference;
    };
 
    template <bool Constness>
    class Iterator
    {
    public:
      using size_type = ::types::global_dof_index;
 
      using difference_type = ::types::global_dof_index;
 
      using value_type = TrilinosScalar;
 
      using MatrixType = typename Accessor<Constness>::MatrixType;
 
      Iterator(MatrixType *matrix, const size_type row, const size_type index);
 
      template <bool Other>
      Iterator(const Iterator<Other> &other);
 
      Iterator<Constness> &
      operator++();
 
      Iterator<Constness>
      operator++(int);
 
      const Accessor<Constness> &
      operator*() const;
 
      const Accessor<Constness> *
      operator->() const;
 
      template <bool OtherConstness>
      bool
      operator==(const Iterator<OtherConstness> &) const;
 
      template <bool OtherConstness>
      bool
      operator!=(const Iterator<OtherConstness> &) const;
 
      template <bool OtherConstness>
      bool
      operator<(const Iterator<OtherConstness> &) const;
 
      template <bool OtherConstness>
      bool
      operator>(const Iterator<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<Constness> accessor;
 
      template <bool Other>
      friend class Iterator;
    };
 
  } // namespace SparseMatrixIterators
} // namespace TrilinosWrappers
 
DEAL_II_NAMESPACE_CLOSE
 
namespace std
{
  template <bool Constness>
  struct iterator_traits<
    ::TrilinosWrappers::SparseMatrixIterators::Iterator<Constness>>
  {
    using iterator_category = forward_iterator_tag;
    using value_type =
      typename ::TrilinosWrappers::SparseMatrixIterators::Iterator<
        Constness>::value_type;
    using difference_type =
      typename ::TrilinosWrappers::SparseMatrixIterators::Iterator<
        Constness>::difference_type;
  };
} // namespace std
 
DEAL_II_NAMESPACE_OPEN
 
 
namespace TrilinosWrappers
{
  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<false>;
 
    using const_iterator = SparseMatrixIterators::Iterator<true>;
 
    using value_type = TrilinosScalar;
 
    SparseMatrix();
 
    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(const SparsityPattern &InputSparsityPattern);
 
    SparseMatrix(SparseMatrix &&other) noexcept;
 
    SparseMatrix(const SparseMatrix &) = delete;
 
    SparseMatrix &
    operator=(const SparseMatrix &) = delete;
 
    virtual ~SparseMatrix() override = default;
 
    template <typename SparsityPatternType>
    void
    reinit(const SparsityPatternType &sparsity_pattern);
 
    void
    reinit(const SparsityPattern &sparsity_pattern);
 
    void
    reinit(const SparseMatrix &sparse_matrix);
 
    template <typename number>
    void
    reinit(const ::SparseMatrix<number> &dealii_sparse_matrix,
           const double                          drop_tolerance    = 1e-13,
           const bool                            copy_values       = true,
           const ::SparsityPattern      *use_this_sparsity = nullptr);
 
    void
    reinit(const Epetra_CrsMatrix &input_matrix, const bool copy_values = true);
    SparseMatrix(const IndexSet    &parallel_partitioning,
                 const MPI_Comm     communicator          = MPI_COMM_WORLD,
                 const unsigned int n_max_entries_per_row = 0);
 
    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          = MPI_COMM_WORLD,
                 const size_type n_max_entries_per_row = 0);
 
    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>
    void
    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);
 
    template <typename number>
    void
    reinit(const IndexSet                       &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);
 
    template <typename number>
    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;
 
    std::uint64_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);
 
    void
    clear();
 
    void
    compress(VectorOperation::values operation);
 
    void
    set(const size_type i, const size_type j, const TrilinosScalar value);
 
    void
    set(const std::vector<size_type>     &indices,
        const FullMatrix<TrilinosScalar> &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<TrilinosScalar> &full_matrix,
        const bool                        elide_zero_values = false);
 
    void
    set(const size_type                    row,
        const std::vector<size_type>      &col_indices,
        const std::vector<TrilinosScalar> &values,
        const bool                         elide_zero_values = false);
 
    template <typename Number>
    void
    set(const size_type  row,
        const size_type  n_cols,
        const size_type *col_indices,
        const Number    *values,
        const bool       elide_zero_values = false);
 
    void
    add(const size_type i, const size_type j, const TrilinosScalar value);
 
    void
    add(const std::vector<size_type>     &indices,
        const FullMatrix<TrilinosScalar> &full_matrix,
        const bool                        elide_zero_values = true);
 
    void
    add(const std::vector<size_type>     &row_indices,
        const std::vector<size_type>     &col_indices,
        const FullMatrix<TrilinosScalar> &full_matrix,
        const bool                        elide_zero_values = true);
 
    void
    add(const size_type                    row,
        const std::vector<size_type>      &col_indices,
        const std::vector<TrilinosScalar> &values,
        const bool                         elide_zero_values = true);
 
    void
    add(const size_type       row,
        const size_type       n_cols,
        const size_type      *col_indices,
        const TrilinosScalar *values,
        const bool            elide_zero_values      = true,
        const bool            col_indices_are_sorted = false);
 
    SparseMatrix &
    operator*=(const TrilinosScalar factor);
 
    SparseMatrix &
    operator/=(const TrilinosScalar factor);
 
    void
    copy_from(const SparseMatrix &source);
 
    void
    add(const TrilinosScalar factor, const SparseMatrix &matrix);
 
    void
    clear_row(const size_type row, const TrilinosScalar new_diag_value = 0);
 
    void
    clear_rows(const ArrayView<const size_type> &rows,
               const TrilinosScalar              new_diag_value = 0);
 
    void
    transpose();
 
    TrilinosScalar
    operator()(const size_type i, const size_type j) const;
 
    TrilinosScalar
    el(const size_type i, const size_type j) const;
 
    TrilinosScalar
    diag_element(const size_type i) const;
 
    template <typename VectorType>
    std::enable_if_t<
      std::is_same_v<typename VectorType::value_type, TrilinosScalar>>
    vmult(VectorType &dst, const VectorType &src) const;
 
    template <typename VectorType>
    std::enable_if_t<
      !std::is_same_v<typename VectorType::value_type, TrilinosScalar>>
    vmult(VectorType &dst, const VectorType &src) const;
 
    template <typename VectorType>
    std::enable_if_t<
      std::is_same_v<typename VectorType::value_type, TrilinosScalar>>
    Tvmult(VectorType &dst, const VectorType &src) const;
 
    template <typename VectorType>
    std::enable_if_t<
      !std::is_same_v<typename VectorType::value_type, TrilinosScalar>>
    Tvmult(VectorType &dst, const VectorType &src) const;
 
    template <typename VectorType>
    void
    vmult_add(VectorType &dst, const VectorType &src) const;
 
    template <typename VectorType>
    void
    Tvmult_add(VectorType &dst, const VectorType &src) const;
 
    TrilinosScalar
    matrix_norm_square(const MPI::Vector &v) const;
 
    TrilinosScalar
    matrix_scalar_product(const MPI::Vector &u, const MPI::Vector &v) const;
 
    template <typename VectorType>
    TrilinosScalar
    residual(VectorType &dst, const VectorType &x, const VectorType &b) const;
 
    void
    mmult(SparseMatrix       &C,
          const SparseMatrix &B,
          const MPI::Vector  &V = MPI::Vector()) const;
 
 
    void
    Tmmult(SparseMatrix       &C,
           const SparseMatrix &B,
           const MPI::Vector  &V = MPI::Vector()) const;
 
    TrilinosScalar
    l1_norm() const;
 
    TrilinosScalar
    linfty_norm() const;
 
    TrilinosScalar
    frobenius_norm() const;
 
    const Epetra_CrsMatrix &
    trilinos_matrix() const;
 
    const Epetra_CrsGraph &
    trilinos_sparsity_pattern() const;
 
    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);
 
    void
    write_ascii();
 
    void
    print(std::ostream &out,
          const bool    write_extended_trilinos_info = false) const;
 
    DeclException1(
      ExcTrilinosError,
      int,
      << "An error with error number " << arg1
      << " occurred while calling a Trilinos function. "
         "\n\n"
         "For historical reasons, many Trilinos functions express "
         "errors by returning specific integer values to indicate "
         "certain errors. Unfortunately, different Trilinos functions "
         "often use the same integer values for different kinds of "
         "errors, and in most cases it is also not documented what "
         "each error code actually means. As a consequence, it is often "
         "difficult to say what a particular error (in this case, "
         "the error with integer code '"
      << arg1
      << "') represents and how one should fix a code to avoid it. "
         "The best one can often do is to look up the call stack to "
         "see which deal.II function generated the error, and which "
         "Trilinos function the error code had originated from; "
         "then look up the Trilinos source code of that function (for "
         "example on github) to see what code path set that error "
         "code. Short of going through all of that, the only other "
         "option is to guess the cause of the error from "
         "the context in which the error appeared.");
 
 
    DeclException2(ExcInvalidIndex,
                   size_type,
                   size_type,
                   << "The entry with index <" << arg1 << ',' << arg2
                   << "> does not exist.");
 
    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.");
 
    DeclException0(ExcMatrixNotCompressed);
 
    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.");
 
    DeclException2(ExcAccessToNonPresentElement,
                   size_type,
                   size_type,
                   << "You tried to access element (" << arg1 << '/' << arg2
                   << ')' << " of a sparse matrix, but it appears to not"
                   << " exist in the Trilinos sparsity pattern.");
  protected:
  private:
    std::unique_ptr<Epetra_Map> column_space_map;
 
    std::unique_ptr<Epetra_FECrsMatrix> matrix;
 
    std::unique_ptr<Epetra_CrsMatrix> nonlocal_matrix;
 
    std::unique_ptr<Epetra_Export> nonlocal_matrix_exporter;
 
    Epetra_CombineMode last_action;
 
    bool compressed;
 
    void
    prepare_add();
 
    void
    prepare_set();
 
    // To allow calling protected prepare_add() and prepare_set().
    friend class BlockMatrixBase<SparseMatrix>;
  };
 
 
 
  // forwards declarations
  class SolverBase;
  class PreconditionBase;
 
  namespace internal
  {
    inline void
    check_vector_map_equality(const Epetra_CrsMatrix   &mtrx,
                              const Epetra_MultiVector &src,
                              const Epetra_MultiVector &dst,
                              const bool                transpose)
    {
      if (transpose == false)
        {
          Assert(src.Map().SameAs(mtrx.DomainMap()) == true,
                 ExcMessage(
                   "Column map of matrix does not fit with vector map!"));
          Assert(dst.Map().SameAs(mtrx.RangeMap()) == true,
                 ExcMessage("Row map of matrix does not fit with vector map!"));
        }
      else
        {
          Assert(src.Map().SameAs(mtrx.RangeMap()) == true,
                 ExcMessage(
                   "Column map of matrix does not fit with vector map!"));
          Assert(dst.Map().SameAs(mtrx.DomainMap()) == true,
                 ExcMessage("Row map of matrix does not fit with vector map!"));
        }
      (void)mtrx; // removes -Wunused-variable in optimized mode
      (void)src;
      (void)dst;
    }
 
    inline void
    check_vector_map_equality(const Epetra_Operator    &op,
                              const Epetra_MultiVector &src,
                              const Epetra_MultiVector &dst,
                              const bool                transpose)
    {
      if (transpose == false)
        {
          Assert(src.Map().SameAs(op.OperatorDomainMap()) == true,
                 ExcMessage(
                   "Column map of operator does not fit with vector map!"));
          Assert(dst.Map().SameAs(op.OperatorRangeMap()) == true,
                 ExcMessage(
                   "Row map of operator does not fit with vector map!"));
        }
      else
        {
          Assert(src.Map().SameAs(op.OperatorRangeMap()) == true,
                 ExcMessage(
                   "Column map of operator does not fit with vector map!"));
          Assert(dst.Map().SameAs(op.OperatorDomainMap()) == true,
                 ExcMessage(
                   "Row map of operator does not fit with vector map!"));
        }
      (void)op; // removes -Wunused-variable in optimized mode
      (void)src;
      (void)dst;
    }
 
 
    namespace LinearOperatorImplementation
    {
      class TrilinosPayload : public Epetra_Operator
      {
      public:
        using VectorType = Epetra_MultiVector;
 
        using Range = VectorType;
 
        using Domain = VectorType;
 
        using size_type = ::types::global_dof_index;
 
        TrilinosPayload();
 
        TrilinosPayload(const TrilinosWrappers::SparseMatrix &matrix_exemplar,
                        const TrilinosWrappers::SparseMatrix &matrix);
 
        TrilinosPayload(const TrilinosPayload                &payload_exemplar,
                        const TrilinosWrappers::SparseMatrix &matrix);
 
        TrilinosPayload(
          const TrilinosWrappers::SparseMatrix     &matrix_exemplar,
          const TrilinosWrappers::PreconditionBase &preconditioner);
 
        TrilinosPayload(
          const TrilinosWrappers::PreconditionBase &preconditioner_exemplar,
          const TrilinosWrappers::PreconditionBase &preconditioner);
 
        TrilinosPayload(
          const TrilinosPayload                    &payload_exemplar,
          const TrilinosWrappers::PreconditionBase &preconditioner);
 
        TrilinosPayload(const TrilinosPayload &payload);
 
        TrilinosPayload(const TrilinosPayload &first_op,
                        const TrilinosPayload &second_op);
 
        virtual ~TrilinosPayload() override = default;
 
        TrilinosPayload
        identity_payload() const;
 
        TrilinosPayload
        null_payload() const;
 
        TrilinosPayload
        transpose_payload() const;
 
        template <typename Solver, typename Preconditioner>
        std::enable_if_t<
          std::is_base_of_v<TrilinosWrappers::SolverBase, Solver> &&
            std::is_base_of_v<TrilinosWrappers::PreconditionBase,
                              Preconditioner>,
          TrilinosPayload>
        inverse_payload(Solver &, const Preconditioner &) const;
 
        template <typename Solver, typename Preconditioner>
        std::enable_if_t<
          !(std::is_base_of_v<TrilinosWrappers::SolverBase, Solver> &&
            std::is_base_of_v<TrilinosWrappers::PreconditionBase,
                              Preconditioner>),
          TrilinosPayload>
        inverse_payload(Solver &, const Preconditioner &) const;
 
        IndexSet
        locally_owned_domain_indices() const;
 
        IndexSet
        locally_owned_range_indices() const;
 
        MPI_Comm
        get_mpi_communicator() const;
 
        void
        transpose();
 
        std::function<void(VectorType &, const VectorType &)> vmult;
 
        std::function<void(VectorType &, const VectorType &)> Tvmult;
 
        std::function<void(VectorType &, const VectorType &)> inv_vmult;
 
        std::function<void(VectorType &, const VectorType &)> inv_Tvmult;
 
        virtual bool
        UseTranspose() const override;
 
        virtual int
        SetUseTranspose(bool UseTranspose) override;
 
        virtual int
        Apply(const VectorType &X, VectorType &Y) const override;
 
        virtual int
        ApplyInverse(const VectorType &Y, VectorType &X) const override;
        virtual const char *
        Label() const override;
 
        virtual const Epetra_Comm &
        Comm() const override;
 
        virtual const Epetra_Map &
        OperatorDomainMap() const override;
 
        virtual const Epetra_Map &
        OperatorRangeMap() const override;
      private:
        template <typename EpetraOpType>
        TrilinosPayload(EpetraOpType   &op,
                        const bool      supports_inverse_operations,
                        const bool      use_transpose,
                        const MPI_Comm  mpi_communicator,
                        const IndexSet &locally_owned_domain_indices,
                        const IndexSet &locally_owned_range_indices);
 
        bool use_transpose;
 
        Epetra_MpiComm communicator;
 
        Epetra_Map domain_map;
 
        Epetra_Map range_map;
 
        virtual bool
        HasNormInf() const override;
 
        virtual double
        NormInf() const override;
      };
 
      TrilinosPayload
      operator+(const TrilinosPayload &first_op,
                const TrilinosPayload &second_op);
 
      TrilinosPayload
      operator*(const TrilinosPayload &first_op,
                const TrilinosPayload &second_op);
 
    } // namespace LinearOperatorImplementation
  }   /* namespace internal */
 
 
 
  // ----------------------- inline and template functions --------------------
 
#    ifndef DOXYGEN
 
  namespace SparseMatrixIterators
  {
    inline AccessorBase::AccessorBase(SparseMatrix *matrix,
                                      size_type     row,
                                      size_type     index)
      : matrix(matrix)
      , a_row(row)
      , a_index(index)
    {
      visit_present_row();
    }
 
 
    inline AccessorBase::size_type
    AccessorBase::row() const
    {
      Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
      return a_row;
    }
 
 
    inline AccessorBase::size_type
    AccessorBase::column() const
    {
      Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
      return (*colnum_cache)[a_index];
    }
 
 
    inline AccessorBase::size_type
    AccessorBase::index() const
    {
      Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
      return a_index;
    }
 
 
    inline Accessor<true>::Accessor(MatrixType     *matrix,
                                    const size_type row,
                                    const size_type index)
      : AccessorBase(const_cast<SparseMatrix *>(matrix), row, index)
    {}
 
 
    template <bool Other>
    inline Accessor<true>::Accessor(const Accessor<Other> &other)
      : AccessorBase(other)
    {}
 
 
    inline TrilinosScalar
    Accessor<true>::value() const
    {
      Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
      return (*value_cache)[a_index];
    }
 
 
    inline Accessor<false>::Reference::Reference(const Accessor<false> &acc)
      : accessor(const_cast<Accessor<false> &>(acc))
    {}
 
 
    inline Accessor<false>::Reference::operator TrilinosScalar() const
    {
      return (*accessor.value_cache)[accessor.a_index];
    }
 
    inline const Accessor<false>::Reference &
    Accessor<false>::Reference::operator=(const TrilinosScalar n) const
    {
      (*accessor.value_cache)[accessor.a_index] = n;
      accessor.matrix->set(accessor.row(),
                           accessor.column(),
                           static_cast<TrilinosScalar>(*this));
      return *this;
    }
 
 
    inline const Accessor<false>::Reference &
    Accessor<false>::Reference::operator+=(const TrilinosScalar n) const
    {
      (*accessor.value_cache)[accessor.a_index] += n;
      accessor.matrix->set(accessor.row(),
                           accessor.column(),
                           static_cast<TrilinosScalar>(*this));
      return *this;
    }
 
 
    inline const Accessor<false>::Reference &
    Accessor<false>::Reference::operator-=(const TrilinosScalar n) const
    {
      (*accessor.value_cache)[accessor.a_index] -= n;
      accessor.matrix->set(accessor.row(),
                           accessor.column(),
                           static_cast<TrilinosScalar>(*this));
      return *this;
    }
 
 
    inline const Accessor<false>::Reference &
    Accessor<false>::Reference::operator*=(const TrilinosScalar n) const
    {
      (*accessor.value_cache)[accessor.a_index] *= n;
      accessor.matrix->set(accessor.row(),
                           accessor.column(),
                           static_cast<TrilinosScalar>(*this));
      return *this;
    }
 
 
    inline const Accessor<false>::Reference &
    Accessor<false>::Reference::operator/=(const TrilinosScalar n) const
    {
      (*accessor.value_cache)[accessor.a_index] /= n;
      accessor.matrix->set(accessor.row(),
                           accessor.column(),
                           static_cast<TrilinosScalar>(*this));
      return *this;
    }
 
 
    inline Accessor<false>::Accessor(MatrixType     *matrix,
                                     const size_type row,
                                     const size_type index)
      : AccessorBase(matrix, row, index)
    {}
 
 
    inline Accessor<false>::Reference
    Accessor<false>::value() const
    {
      Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
      return {*this};
    }
 
 
 
    template <bool Constness>
    inline Iterator<Constness>::Iterator(MatrixType     *matrix,
                                         const size_type row,
                                         const size_type index)
      : accessor(matrix, row, index)
    {}
 
 
    template <bool Constness>
    template <bool Other>
    inline Iterator<Constness>::Iterator(const Iterator<Other> &other)
      : accessor(other.accessor)
    {}
 
 
    template <bool Constness>
    inline Iterator<Constness> &
    Iterator<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 <bool Constness>
    inline Iterator<Constness>
    Iterator<Constness>::operator++(int)
    {
      const Iterator<Constness> old_state = *this;
      ++(*this);
      return old_state;
    }
 
 
 
    template <bool Constness>
    inline const Accessor<Constness> &
    Iterator<Constness>::operator*() const
    {
      return accessor;
    }
 
 
 
    template <bool Constness>
    inline const Accessor<Constness> *
    Iterator<Constness>::operator->() const
    {
      return &accessor;
    }
 
 
 
    template <bool Constness>
    template <bool OtherConstness>
    inline bool
    Iterator<Constness>::operator==(const Iterator<OtherConstness> &other) const
    {
      return (accessor.a_row == other.accessor.a_row &&
              accessor.a_index == other.accessor.a_index);
    }
 
 
 
    template <bool Constness>
    template <bool OtherConstness>
    inline bool
    Iterator<Constness>::operator!=(const Iterator<OtherConstness> &other) const
    {
      return !(*this == other);
    }
 
 
 
    template <bool Constness>
    template <bool OtherConstness>
    inline bool
    Iterator<Constness>::operator<(const Iterator<OtherConstness> &other) const
    {
      return (accessor.row() < other.accessor.row() ||
              (accessor.row() == other.accessor.row() &&
               accessor.index() < other.accessor.index()));
    }
 
 
    template <bool Constness>
    template <bool OtherConstness>
    inline bool
    Iterator<Constness>::operator>(const Iterator<OtherConstness> &other) const
    {
      return (other < *this);
    }
 
  } // namespace SparseMatrixIterators
 
 
 
  inline SparseMatrix::const_iterator
  SparseMatrix::begin() const
  {
    return begin(0);
  }
 
 
 
  inline SparseMatrix::const_iterator
  SparseMatrix::end() const
  {
    return const_iterator(this, m(), 0);
  }
 
 
 
  inline SparseMatrix::const_iterator
  SparseMatrix::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);
  }
 
 
 
  inline SparseMatrix::const_iterator
  SparseMatrix::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();
  }
 
 
 
  inline SparseMatrix::iterator
  SparseMatrix::begin()
  {
    return begin(0);
  }
 
 
 
  inline SparseMatrix::iterator
  SparseMatrix::end()
  {
    return iterator(this, m(), 0);
  }
 
 
 
  inline SparseMatrix::iterator
  SparseMatrix::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);
  }
 
 
 
  inline SparseMatrix::iterator
  SparseMatrix::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();
  }
 
 
 
  inline bool
  SparseMatrix::in_local_range(const size_type index) const
  {
    TrilinosWrappers::types::int_type begin, end;
#      ifndef DEAL_II_WITH_64BIT_INDICES
    begin = matrix->RowMap().MinMyGID();
    end   = matrix->RowMap().MaxMyGID() + 1;
#      else
    begin = matrix->RowMap().MinMyGID64();
    end   = matrix->RowMap().MaxMyGID64() + 1;
#      endif
 
    return ((index >= static_cast<size_type>(begin)) &&
            (index < static_cast<size_type>(end)));
  }
 
 
 
  inline bool
  SparseMatrix::is_compressed() const
  {
    return compressed;
  }
 
 
 
  // Inline the set() and add() functions, since they will be called
  // frequently, and the compiler can optimize away some unnecessary loops
  // when the sizes are given at compile time.
  template <>
  void
  SparseMatrix::set<TrilinosScalar>(const size_type       row,
                                    const size_type       n_cols,
                                    const size_type      *col_indices,
                                    const TrilinosScalar *values,
                                    const bool            elide_zero_values);
 
 
 
  template <typename Number>
  void
  SparseMatrix::set(const size_type  row,
                    const size_type  n_cols,
                    const size_type *col_indices,
                    const Number    *values,
                    const bool       elide_zero_values)
  {
    std::vector<TrilinosScalar> trilinos_values(n_cols);
    std::copy(values, values + n_cols, trilinos_values.begin());
    this->set(
      row, n_cols, col_indices, trilinos_values.data(), elide_zero_values);
  }
 
 
 
  inline void
  SparseMatrix::set(const size_type      i,
                    const size_type      j,
                    const TrilinosScalar value)
  {
    AssertIsFinite(value);
 
    set(i, 1, &j, &value, false);
  }
 
 
 
  inline void
  SparseMatrix::set(const std::vector<size_type>     &indices,
                    const FullMatrix<TrilinosScalar> &values,
                    const bool                        elide_zero_values)
  {
    Assert(indices.size() == values.m(),
           ExcDimensionMismatch(indices.size(), values.m()));
    Assert(values.m() == values.n(), ExcNotQuadratic());
 
    for (size_type i = 0; i < indices.size(); ++i)
      set(indices[i],
          indices.size(),
          indices.data(),
          &values(i, 0),
          elide_zero_values);
  }
 
 
 
  inline void
  SparseMatrix::add(const size_type      i,
                    const size_type      j,
                    const TrilinosScalar value)
  {
    AssertIsFinite(value);
 
    if (value == 0)
      {
        // we have to check after Insert/Add in any case to be consistent
        // with the MPI communication model, but we can save some
        // work if the addend is zero. However, these actions are done in case
        // we pass on to the other function.
 
        // TODO: fix this (do not run compress here, but fail)
        if (last_action == Insert)
          {
            const int ierr = matrix->GlobalAssemble(*column_space_map,
                                                    matrix->RowMap(),
                                                    false);
 
            AssertThrow(ierr == 0, ExcTrilinosError(ierr));
          }
 
        last_action = Add;
 
        return;
      }
    else
      add(i, 1, &j, &value, false);
  }
 
 
 
  // inline "simple" functions that are called frequently and do only involve
  // a call to some Trilinos function.
  inline SparseMatrix::size_type
  SparseMatrix::m() const
  {
#      ifndef DEAL_II_WITH_64BIT_INDICES
    return matrix->NumGlobalRows();
#      else
    return matrix->NumGlobalRows64();
#      endif
  }
 
 
 
  inline SparseMatrix::size_type
  SparseMatrix::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 n_global_elements(*column_space_map);
  }
 
 
 
  inline unsigned int
  SparseMatrix::local_size() const
  {
    return matrix->NumMyRows();
  }
 
 
 
  inline std::pair<SparseMatrix::size_type, SparseMatrix::size_type>
  SparseMatrix::local_range() const
  {
    size_type begin, end;
#      ifndef DEAL_II_WITH_64BIT_INDICES
    begin = matrix->RowMap().MinMyGID();
    end   = matrix->RowMap().MaxMyGID() + 1;
#      else
    begin = matrix->RowMap().MinMyGID64();
    end   = matrix->RowMap().MaxMyGID64() + 1;
#      endif
 
    return std::make_pair(begin, end);
  }
 
 
 
  inline std::uint64_t
  SparseMatrix::n_nonzero_elements() const
  {
    // Trilinos uses 64bit functions internally for attribute access, which
    // return `long long`. They also offer 32bit variants that return `int`,
    // however those call the 64bit version and convert the values to 32bit.
    // There is no necessity in using the 32bit versions at all.
    return static_cast<std::uint64_t>(matrix->NumGlobalNonzeros64());
  }
 
 
 
  template <typename SparsityPatternType>
  inline void
  SparseMatrix::reinit(const IndexSet            &parallel_partitioning,
                       const SparsityPatternType &sparsity_pattern,
                       const MPI_Comm             communicator,
                       const bool                 exchange_data)
  {
    reinit(parallel_partitioning,
           parallel_partitioning,
           sparsity_pattern,
           communicator,
           exchange_data);
  }
 
 
 
  template <typename number>
  inline void
  SparseMatrix::reinit(const IndexSet &parallel_partitioning,
                       const ::SparseMatrix<number> &sparse_matrix,
                       const MPI_Comm                        communicator,
                       const double                          drop_tolerance,
                       const bool                            copy_values,
                       const ::SparsityPattern      *use_this_sparsity)
  {
    Epetra_Map map =
      parallel_partitioning.make_trilinos_map(communicator, false);
    reinit(parallel_partitioning,
           parallel_partitioning,
           sparse_matrix,
           drop_tolerance,
           copy_values,
           use_this_sparsity);
  }
 
 
 
  inline const Epetra_CrsMatrix &
  SparseMatrix::trilinos_matrix() const
  {
    return static_cast<const Epetra_CrsMatrix &>(*matrix);
  }
 
 
 
  inline const Epetra_CrsGraph &
  SparseMatrix::trilinos_sparsity_pattern() const
  {
    return matrix->Graph();
  }
 
 
 
  inline IndexSet
  SparseMatrix::locally_owned_domain_indices() const
  {
    return IndexSet(matrix->DomainMap());
  }
 
 
 
  inline IndexSet
  SparseMatrix::locally_owned_range_indices() const
  {
    return IndexSet(matrix->RangeMap());
  }
 
 
 
  inline void
  SparseMatrix::prepare_add()
  {
    // nothing to do here
  }
 
 
 
  inline void
  SparseMatrix::prepare_set()
  {
    // nothing to do here
  }
 
 
 
  template <typename VectorType>
  inline TrilinosScalar
  SparseMatrix::residual(VectorType       &dst,
                         const VectorType &x,
                         const VectorType &b) const
  {
    vmult(dst, x);
    dst -= b;
    dst *= -1.;
 
    return dst.l2_norm();
  }
 
 
  namespace internal
  {
    namespace LinearOperatorImplementation
    {
      template <typename EpetraOpType>
      TrilinosPayload::TrilinosPayload(
        EpetraOpType   &op,
        const bool      supports_inverse_operations,
        const bool      use_transpose,
        const MPI_Comm  mpi_communicator,
        const IndexSet &locally_owned_domain_indices,
        const IndexSet &locally_owned_range_indices)
        : use_transpose(use_transpose)
        , communicator(mpi_communicator)
        , domain_map(
            locally_owned_domain_indices.make_trilinos_map(communicator.Comm()))
        , range_map(
            locally_owned_range_indices.make_trilinos_map(communicator.Comm()))
      {
        vmult = [&op](Range &tril_dst, const Domain &tril_src) {
          // Duplicated from TrilinosWrappers::PreconditionBase::vmult
          // as well as from TrilinosWrappers::SparseMatrix::Tvmult
          Assert(&tril_src != &tril_dst,
                 TrilinosWrappers::SparseMatrix::ExcSourceEqualsDestination());
          internal::check_vector_map_equality(op,
                                              tril_src,
                                              tril_dst,
                                              op.UseTranspose());
 
          const int ierr = op.Apply(tril_src, tril_dst);
          AssertThrow(ierr == 0, ExcTrilinosError(ierr));
        };
 
        Tvmult = [&op](Domain &tril_dst, const Range &tril_src) {
          // Duplicated from TrilinosWrappers::PreconditionBase::vmult
          // as well as from TrilinosWrappers::SparseMatrix::Tvmult
          Assert(&tril_src != &tril_dst,
                 TrilinosWrappers::SparseMatrix::ExcSourceEqualsDestination());
          internal::check_vector_map_equality(op,
                                              tril_src,
                                              tril_dst,
                                              !op.UseTranspose());
 
          op.SetUseTranspose(!op.UseTranspose());
          const int ierr = op.Apply(tril_src, tril_dst);
          AssertThrow(ierr == 0, ExcTrilinosError(ierr));
          op.SetUseTranspose(!op.UseTranspose());
        };
 
        if (supports_inverse_operations)
          {
            inv_vmult = [&op](Domain &tril_dst, const Range &tril_src) {
              // Duplicated from TrilinosWrappers::PreconditionBase::vmult
              // as well as from TrilinosWrappers::SparseMatrix::Tvmult
              Assert(
                &tril_src != &tril_dst,
                TrilinosWrappers::SparseMatrix::ExcSourceEqualsDestination());
              internal::check_vector_map_equality(op,
                                                  tril_src,
                                                  tril_dst,
                                                  !op.UseTranspose());
 
              const int ierr = op.ApplyInverse(tril_src, tril_dst);
              AssertThrow(ierr == 0, ExcTrilinosError(ierr));
            };
 
            inv_Tvmult = [&op](Range &tril_dst, const Domain &tril_src) {
              // Duplicated from TrilinosWrappers::PreconditionBase::vmult
              // as well as from TrilinosWrappers::SparseMatrix::Tvmult
              Assert(
                &tril_src != &tril_dst,
                TrilinosWrappers::SparseMatrix::ExcSourceEqualsDestination());
              internal::check_vector_map_equality(op,
                                                  tril_src,
                                                  tril_dst,
                                                  op.UseTranspose());
 
              op.SetUseTranspose(!op.UseTranspose());
              const int ierr = op.ApplyInverse(tril_src, tril_dst);
              AssertThrow(ierr == 0, ExcTrilinosError(ierr));
              op.SetUseTranspose(!op.UseTranspose());
            };
          }
        else
          {
            inv_vmult = [](Domain &, const Range &) {
              Assert(false,
                     ExcMessage(
                       "Uninitialized TrilinosPayload::inv_vmult called. "
                       "The operator does not support inverse operations."));
            };
 
            inv_Tvmult = [](Range &, const Domain &) {
              Assert(false,
                     ExcMessage(
                       "Uninitialized TrilinosPayload::inv_Tvmult called. "
                       "The operator does not support inverse operations."));
            };
          }
      }
 
 
      template <typename Solver, typename Preconditioner>
      std::enable_if_t<
        std::is_base_of_v<TrilinosWrappers::SolverBase, Solver> &&
          std::is_base_of_v<TrilinosWrappers::PreconditionBase, Preconditioner>,
        TrilinosPayload>
      TrilinosPayload::inverse_payload(
        Solver               &solver,
        const Preconditioner &preconditioner) const
      {
        const auto &payload = *this;
 
        TrilinosPayload return_op(payload);
 
        // Capture by copy so the payloads are always valid
 
        return_op.inv_vmult = [payload, &solver, &preconditioner](
                                TrilinosPayload::Domain      &tril_dst,
                                const TrilinosPayload::Range &tril_src) {
          // Duplicated from TrilinosWrappers::PreconditionBase::vmult
          // as well as from TrilinosWrappers::SparseMatrix::Tvmult
          Assert(&tril_src != &tril_dst,
                 TrilinosWrappers::SparseMatrix::ExcSourceEqualsDestination());
          internal::check_vector_map_equality(payload,
                                              tril_src,
                                              tril_dst,
                                              !payload.UseTranspose());
          solver.solve(payload, tril_dst, tril_src, preconditioner);
        };
 
        return_op.inv_Tvmult = [payload, &solver, &preconditioner](
                                 TrilinosPayload::Range        &tril_dst,
                                 const TrilinosPayload::Domain &tril_src) {
          // Duplicated from TrilinosWrappers::PreconditionBase::vmult
          // as well as from TrilinosWrappers::SparseMatrix::Tvmult
          Assert(&tril_src != &tril_dst,
                 TrilinosWrappers::SparseMatrix::ExcSourceEqualsDestination());
          internal::check_vector_map_equality(payload,
                                              tril_src,
                                              tril_dst,
                                              payload.UseTranspose());
 
          const_cast<TrilinosPayload &>(payload).transpose();
          solver.solve(payload, tril_dst, tril_src, preconditioner);
          const_cast<TrilinosPayload &>(payload).transpose();
        };
 
        // If the input operator is already setup for transpose operations, then
        // we must do similar with its inverse.
        if (return_op.UseTranspose() == true)
          std::swap(return_op.inv_vmult, return_op.inv_Tvmult);
 
        return return_op;
      }
 
      template <typename Solver, typename Preconditioner>
      std::enable_if_t<
        !(std::is_base_of_v<TrilinosWrappers::SolverBase, Solver> &&
          std::is_base_of_v<TrilinosWrappers::PreconditionBase,
                            Preconditioner>),
        TrilinosPayload>
      TrilinosPayload::inverse_payload(Solver &, const Preconditioner &) const
      {
        TrilinosPayload return_op(*this);
 
        return_op.inv_vmult = [](TrilinosPayload::Domain &,
                                 const TrilinosPayload::Range &) {
          AssertThrow(false,
                      ExcMessage("Payload inv_vmult disabled because of "
                                 "incompatible solver/preconditioner choice."));
        };
 
        return_op.inv_Tvmult = [](TrilinosPayload::Range &,
                                  const TrilinosPayload::Domain &) {
          AssertThrow(false,
                      ExcMessage("Payload inv_vmult disabled because of "
                                 "incompatible solver/preconditioner choice."));
        };
 
        return return_op;
      }
    } // namespace LinearOperatorImplementation
  }   // namespace internal
 
  template <>
  void
  SparseMatrix::set<TrilinosScalar>(const size_type       row,
                                    const size_type       n_cols,
                                    const size_type      *col_indices,
                                    const TrilinosScalar *values,
                                    const bool            elide_zero_values);
#    endif // DOXYGEN
 
} /* namespace TrilinosWrappers */
 
 
DEAL_II_NAMESPACE_CLOSE
 
 
#  endif // DEAL_II_WITH_TRILINOS
 
 
/*-----------------------   trilinos_sparse_matrix.h     --------------------*/
 
#endif
/*-----------------------   trilinos_sparse_matrix.h     --------------------*/