trilinos_sparse_matrix.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2008 - 2018 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 at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#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/subscriptor.h>
#  include <deal.II/base/index_set.h>
#  include <deal.II/lac/full_matrix.h>
#  include <deal.II/lac/exceptions.h>
#  include <deal.II/lac/trilinos_vector.h>
#  include <deal.II/lac/trilinos_epetra_vector.h>
#  include <deal.II/lac/vector_memory.h>
#include <deal.II/lac/vector_operation.h>

#  include <type_traits>
#  include <vector>
#  include <cmath>
#  include <memory>

#  include <Epetra_FECrsMatrix.h>
#  include <Epetra_Export.h>
#  include <Epetra_Map.h>
#  include <Epetra_CrsGraph.h>
#  include <Epetra_MultiVector.h>
#  include <Epetra_Operator.h>
#  include <Epetra_Comm.h>
#  ifdef DEAL_II_WITH_MPI
#    include <Epetra_MpiComm.h>
#    include <mpi.h>
#  else
#    include <Epetra_SerialComm.h>
#  endif

DEAL_II_NAMESPACE_OPEN

// forward declarations
template <typename MatrixType> class BlockMatrixBase;

template <typename number> class SparseMatrix;
class SparsityPattern;
class DynamicSparsityPattern;



namespace TrilinosWrappers
{
  // forward declarations
  class SparseMatrix;
  class SparsityPattern;

  namespace SparseMatrixIterators
  {
    // forward declaration
    template <bool Constness> class Iterator;

    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:
      typedef ::types::global_dof_index size_type;

      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:
      typedef const SparseMatrix MatrixType;

      Accessor (MatrixType *matrix,
                const size_type  row,
                const size_type  index);

      template <bool Other>
      Accessor (const Accessor<Other> &a);

      TrilinosScalar value() const;

    private:
      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:
      typedef SparseMatrix MatrixType;

      Accessor (MatrixType *matrix,
                const size_type  row,
                const size_type  index);

      Reference value() const;

    private:
      template <bool> friend class Iterator;
      friend class Reference;
    };

    template <bool Constness>
    class Iterator
    {
    public:
      typedef ::types::global_dof_index size_type;

      typedef typename Accessor<Constness>::MatrixType 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;
    };

  }


  class SparseMatrix : public Subscriptor
  {
  public:
    typedef ::types::global_dof_index size_type;

    struct Traits
    {
      static const bool zero_addition_can_be_elided = true;
    };

    typedef SparseMatrixIterators::Iterator<false> iterator;

    typedef SparseMatrixIterators::Iterator<true> const_iterator;

    typedef TrilinosScalar value_type;

    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 () = 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);

    DEAL_II_DEPRECATED
    SparseMatrix (const Epetra_Map  &parallel_partitioning,
                  const size_type    n_max_entries_per_row = 0);

    DEAL_II_DEPRECATED
    SparseMatrix (const Epetra_Map                &parallel_partitioning,
                  const std::vector<unsigned int> &n_entries_per_row);

    DEAL_II_DEPRECATED
    SparseMatrix (const Epetra_Map &row_parallel_partitioning,
                  const Epetra_Map &col_parallel_partitioning,
                  const size_type   n_max_entries_per_row = 0);

    DEAL_II_DEPRECATED
    SparseMatrix (const Epetra_Map                &row_parallel_partitioning,
                  const Epetra_Map                &col_parallel_partitioning,
                  const std::vector<unsigned int> &n_entries_per_row);

    template <typename SparsityPatternType>
    DEAL_II_DEPRECATED
    void reinit (const Epetra_Map          &parallel_partitioning,
                 const SparsityPatternType &sparsity_pattern,
                 const bool                 exchange_data = false);

    template <typename SparsityPatternType>
    DEAL_II_DEPRECATED
    void reinit (const Epetra_Map          &row_parallel_partitioning,
                 const Epetra_Map          &col_parallel_partitioning,
                 const SparsityPatternType &sparsity_pattern,
                 const bool                 exchange_data = false);

    template <typename number>
    DEAL_II_DEPRECATED
    void reinit (const Epetra_Map                     &parallel_partitioning,
                 const ::SparseMatrix<number> &dealii_sparse_matrix,
                 const double                          drop_tolerance=1e-13,
                 const bool                            copy_values=true,
                 const ::SparsityPattern      *use_this_sparsity=nullptr);

    template <typename number>
    DEAL_II_DEPRECATED
    void reinit (const Epetra_Map                      &row_parallel_partitioning,
                 const Epetra_Map                      &col_parallel_partitioning,
                 const ::SparseMatrix<number>  &dealii_sparse_matrix,
                 const double                           drop_tolerance=1e-13,
                 const bool                             copy_values=true,
                 const ::SparsityPattern      *use_this_sparsity=nullptr);

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

    void set (const size_type       row,
              const size_type       n_cols,
              const size_type      *col_indices,
              const TrilinosScalar *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>
    void vmult (VectorType       &dst,
                const VectorType &src) const;

    template <typename VectorType>
    void 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;

    DEAL_II_DEPRECATED
    const Epetra_Map &domain_partitioner ()  const;

    DEAL_II_DEPRECATED
    const Epetra_Map &range_partitioner () const;

    DEAL_II_DEPRECATED
    const Epetra_Map &row_partitioner () const;

    DEAL_II_DEPRECATED
    const Epetra_Map &col_partitioner () 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;

    friend class BlockMatrixBase<SparseMatrix>;
  };



  // forwards declarations
  class SolverBase;
  class PreconditionBase;

  namespace internal
  {
    namespace
    {
      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:

        typedef Epetra_MultiVector VectorType;

        typedef VectorType Range;

        typedef VectorType Domain;


        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() = 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;

        virtual int
        SetUseTranspose (bool UseTranspose);

        virtual int
        Apply(const VectorType &X,
              VectorType       &Y) const;

        virtual int
        ApplyInverse(const VectorType &Y,
                     VectorType       &X) const;


        virtual const char *
        Label () const;

        virtual const Epetra_Comm &
        Comm () const;

        virtual const Epetra_Map &
        OperatorDomainMap () const;

        virtual const Epetra_Map &
        OperatorRangeMap () const;

      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;

        virtual double
        NormInf () const;
      };

      TrilinosPayload operator+(const TrilinosPayload &first_op,
                                const TrilinosPayload &second_op);

      TrilinosPayload operator*(const TrilinosPayload &first_op,
                                const TrilinosPayload &second_op);

    } /* namespace LinearOperator */
  } /* 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 Reference(*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);
    }

  }



  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
  {
    Assert (r < m(), ExcIndexRange(r, 0, 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
  {
    Assert (r < m(), ExcIndexRange(r, 0, 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)
  {
    Assert (r < m(), ExcIndexRange(r, 0, 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)
  {
    Assert (r < m(), ExcIndexRange(r, 0, 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.
  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());
#ifndef DEAL_II_WITH_64BIT_INDICES
    return column_space_map->NumGlobalElements();
#else
    return column_space_map->NumGlobalElements64();
#endif
  }



  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 LinearOperator
  } // namespace internal

#endif // DOXYGEN

} /* namespace TrilinosWrappers */


DEAL_II_NAMESPACE_CLOSE


#endif // DEAL_II_WITH_TRILINOS


/*-----------------------   trilinos_sparse_matrix.h     --------------------*/

#endif
/*-----------------------   trilinos_sparse_matrix.h     --------------------*/