trilinos_sparse_matrix.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 2008 - 2020 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.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/index_set.h>
#    include <deal.II/base/subscriptor.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_MultiVector.h>
#    include <Epetra_Operator.h>

#    include <cmath>
#    include <memory>
#    include <type_traits>
#    include <vector>
#    ifdef DEAL_II_WITH_MPI
#      include <Epetra_MpiComm.h>
#      include <mpi.h>
#    else
#      include <Epetra_SerialComm.h>
#    endif

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

      bool
      operator==(const Iterator<Constness> &) const;

      bool
      operator!=(const Iterator<Constness> &) const;

      bool
      operator<(const Iterator<Constness> &) const;

      bool
      operator>(const Iterator<Constness> &) 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


  class SparseMatrix : public Subscriptor
  {
  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>
    typename std::enable_if<
      !std::is_same<SparsityPatternType,
                    ::SparseMatrix<double>>::value>::type
    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;

    size_type
    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 std::vector<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>
    typename std::enable_if<std::is_same<typename VectorType::value_type,
                                         TrilinosScalar>::value>::type
    vmult(VectorType &dst, const VectorType &src) const;

    template <typename VectorType>
    typename std::enable_if<!std::is_same<typename VectorType::value_type,
                                          TrilinosScalar>::value>::type
    vmult(VectorType &dst, const VectorType &src) const;

    template <typename VectorType>
    typename std::enable_if<std::is_same<typename VectorType::value_type,
                                         TrilinosScalar>::value>::type
    Tvmult(VectorType &dst, const VectorType &src) const;

    template <typename VectorType>
    typename std::enable_if<!std::is_same<typename VectorType::value_type,
                                          TrilinosScalar>::value>::type
    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;

    TrilinosScalar
    residual(MPI::Vector &      dst,
             const MPI::Vector &x,
             const MPI::Vector &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");

    DeclException2(ExcInvalidIndex,
                   size_type,
                   size_type,
                   << "The entry with index <" << arg1 << ',' << arg2
                   << "> does not exist.");

    DeclExceptionMsg(ExcSourceEqualsDestination,
                     "You are attempting an operation on two matrices 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 " << arg3
                   << " through " << 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:
    void
    prepare_add();

    void
    prepare_set();



  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;

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


        TrilinosPayload();

        TrilinosPayload(const TrilinosWrappers::SparseMatrix &matrix_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);

        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>
        typename std::enable_if<
          std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
            std::is_base_of<TrilinosWrappers::PreconditionBase,
                            Preconditioner>::value,
          TrilinosPayload>::type
        inverse_payload(Solver &, const Preconditioner &) const;

        template <typename Solver, typename Preconditioner>
        typename std::enable_if<
          !(std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
            std::is_base_of<TrilinosWrappers::PreconditionBase,
                            Preconditioner>::value),
          TrilinosPayload>::type
        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:
        bool use_transpose;

#    ifdef DEAL_II_WITH_MPI
        Epetra_MpiComm communicator;
#    else
        Epetra_SerialComm communicator;
#    endif

        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>
    inline bool
    Iterator<Constness>::operator==(const Iterator<Constness> &other) const
    {
      return (accessor.a_row == other.accessor.a_row &&
              accessor.a_index == other.accessor.a_index);
    }



    template <bool Constness>
    inline bool
    Iterator<Constness>::operator!=(const Iterator<Constness> &other) const
    {
      return !(*this == other);
    }



    template <bool Constness>
    inline bool
    Iterator<Constness>::operator<(const Iterator<Constness> &other) const
    {
      return (accessor.row() < other.accessor.row() ||
              (accessor.row() == other.accessor.row() &&
               accessor.index() < other.accessor.index()));
    }


    template <bool Constness>
    inline bool
    Iterator<Constness>::operator>(const Iterator<Constness> &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)
          {
            int ierr;
            ierr = matrix->GlobalAssemble(*column_space_map,
                                          matrix->RowMap(),
                                          false);

            Assert(ierr == 0, ExcTrilinosError(ierr));
            (void)ierr; // removes -Wunused-but-set-variable in optimized mode
          }

        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 SparseMatrix::size_type
  SparseMatrix::n_nonzero_elements() const
  {
#      ifndef DEAL_II_WITH_64BIT_INDICES
    return matrix->NumGlobalNonzeros();
#      else
    return matrix->NumGlobalNonzeros64();
#      endif
  }



  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
  }


  namespace internal
  {
    namespace LinearOperatorImplementation
    {
      template <typename Solver, typename Preconditioner>
      typename std::enable_if<
        std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
          std::is_base_of<TrilinosWrappers::PreconditionBase,
                          Preconditioner>::value,
        TrilinosPayload>::type
      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>
      typename std::enable_if<
        !(std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
          std::is_base_of<TrilinosWrappers::PreconditionBase,
                          Preconditioner>::value),
        TrilinosPayload>::type
      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     --------------------*/