trilinos_sparsity_pattern.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_sparsity_pattern_h
#define dealii_trilinos_sparsity_pattern_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/exceptions.h>

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


#  include <Epetra_FECrsGraph.h>
#  include <Epetra_Map.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
class SparsityPattern;
class DynamicSparsityPattern;

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

  namespace SparsityPatternIterators
  {
    // forward declaration
    class Iterator;

    class Accessor
    {
    public:
      typedef ::types::global_dof_index size_type;

      Accessor (const SparsityPattern *sparsity_pattern,
                const size_type        row,
                const size_type        index);

      size_type row() const;

      size_type index() const;

      size_type column() const;

      DeclException0 (ExcBeyondEndOfSparsityPattern);

      DeclException3 (ExcAccessToNonlocalRow,
                      size_type, size_type, size_type,
                      << "You tried to access row " << arg1
                      << " of a distributed sparsity pattern, "
                      << " but only rows " << arg2 << " through " << arg3
                      << " are stored locally and can be accessed.");

    private:
      mutable SparsityPattern *sparsity_pattern;

      size_type a_row;

      size_type a_index;

      std::shared_ptr<const std::vector<size_type> > colnum_cache;

      void visit_present_row ();

      friend class Iterator;
    };

    class Iterator
    {
    public:
      typedef ::types::global_dof_index size_type;

      Iterator (const SparsityPattern *sparsity_pattern,
                const size_type        row,
                const size_type        index);

      Iterator (const Iterator &i);

      Iterator &operator++ ();

      Iterator operator++ (int);

      const Accessor &operator* () const;

      const Accessor *operator-> () const;

      bool operator == (const Iterator &) const;

      bool operator != (const Iterator &) const;

      bool operator < (const Iterator &) const;

      DeclException2 (ExcInvalidIndexWithinRow,
                      size_type, size_type,
                      << "Attempt to access element " << arg2
                      << " of row " << arg1
                      << " which doesn't have that many elements.");

    private:
      Accessor accessor;

      friend class TrilinosWrappers::SparsityPattern;
    };

  }


  class SparsityPattern : public Subscriptor
  {
  public:

    typedef ::types::global_dof_index size_type;

    typedef SparsityPatternIterators::Iterator const_iterator;

    SparsityPattern ();

    SparsityPattern (const size_type  m,
                     const size_type  n,
                     const size_type  n_entries_per_row = 0);

    SparsityPattern (const size_type               m,
                     const size_type               n,
                     const std::vector<size_type> &n_entries_per_row);

    SparsityPattern (SparsityPattern &&other) noexcept;

    SparsityPattern (const SparsityPattern &input_sparsity_pattern);

    virtual ~SparsityPattern () = default;

    void
    reinit (const size_type  m,
            const size_type  n,
            const size_type  n_entries_per_row = 0);

    void
    reinit (const size_type               m,
            const size_type               n,
            const std::vector<size_type> &n_entries_per_row);

    void
    copy_from (const SparsityPattern &input_sparsity_pattern);

    template <typename SparsityPatternType>
    void
    copy_from (const SparsityPatternType &nontrilinos_sparsity_pattern);

    SparsityPattern &operator = (const SparsityPattern &input_sparsity_pattern);

    void clear ();

    void compress ();


    DEAL_II_DEPRECATED
    SparsityPattern (const Epetra_Map &parallel_partitioning,
                     const size_type   n_entries_per_row = 0);

    DEAL_II_DEPRECATED
    SparsityPattern (const Epetra_Map             &parallel_partitioning,
                     const std::vector<size_type> &n_entries_per_row);

    DEAL_II_DEPRECATED
    SparsityPattern (const Epetra_Map   &row_parallel_partitioning,
                     const Epetra_Map   &col_parallel_partitioning,
                     const size_type     n_entries_per_row = 0);

    DEAL_II_DEPRECATED
    SparsityPattern (const Epetra_Map             &row_parallel_partitioning,
                     const Epetra_Map             &col_parallel_partitioning,
                     const std::vector<size_type> &n_entries_per_row);

    DEAL_II_DEPRECATED
    void
    reinit (const Epetra_Map &parallel_partitioning,
            const size_type   n_entries_per_row = 0);

    DEAL_II_DEPRECATED
    void
    reinit (const Epetra_Map             &parallel_partitioning,
            const std::vector<size_type> &n_entries_per_row);

    DEAL_II_DEPRECATED
    void
    reinit (const Epetra_Map   &row_parallel_partitioning,
            const Epetra_Map   &col_parallel_partitioning,
            const size_type     n_entries_per_row = 0);

    DEAL_II_DEPRECATED
    void
    reinit (const Epetra_Map             &row_parallel_partitioning,
            const Epetra_Map             &col_parallel_partitioning,
            const std::vector<size_type> &n_entries_per_row);

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

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


    SparsityPattern (const IndexSet  &parallel_partitioning,
                     const MPI_Comm  &communicator = MPI_COMM_WORLD,
                     const size_type  n_entries_per_row = 0);

    SparsityPattern (const IndexSet                  &parallel_partitioning,
                     const MPI_Comm                  &communicator,
                     const std::vector<size_type> &n_entries_per_row);

    SparsityPattern (const IndexSet  &row_parallel_partitioning,
                     const IndexSet  &col_parallel_partitioning,
                     const MPI_Comm  &communicator = MPI_COMM_WORLD,
                     const size_type  n_entries_per_row = 0);

    SparsityPattern (const IndexSet               &row_parallel_partitioning,
                     const IndexSet               &col_parallel_partitioning,
                     const MPI_Comm               &communicator,
                     const std::vector<size_type> &n_entries_per_row);

    SparsityPattern (const IndexSet  &row_parallel_partitioning,
                     const IndexSet  &col_parallel_partitioning,
                     const IndexSet  &writable_rows,
                     const MPI_Comm  &communicator = MPI_COMM_WORLD,
                     const size_type  n_entries_per_row = 0);

    void
    reinit (const IndexSet  &parallel_partitioning,
            const MPI_Comm  &communicator = MPI_COMM_WORLD,
            const size_type  n_entries_per_row = 0);

    void
    reinit (const IndexSet               &parallel_partitioning,
            const MPI_Comm               &communicator,
            const std::vector<size_type> &n_entries_per_row);

    void
    reinit (const IndexSet  &row_parallel_partitioning,
            const IndexSet  &col_parallel_partitioning,
            const MPI_Comm  &communicator = MPI_COMM_WORLD,
            const size_type  n_entries_per_row = 0);

    void
    reinit (const IndexSet  &row_parallel_partitioning,
            const IndexSet  &col_parallel_partitioning,
            const IndexSet  &writeable_rows,
            const MPI_Comm  &communicator = MPI_COMM_WORLD,
            const size_type  n_entries_per_row = 0);

    void
    reinit (const IndexSet               &row_parallel_partitioning,
            const IndexSet               &col_parallel_partitioning,
            const MPI_Comm               &communicator,
            const std::vector<size_type> &n_entries_per_row);

    template <typename SparsityPatternType>
    void
    reinit (const IndexSet            &row_parallel_partitioning,
            const IndexSet            &col_parallel_partitioning,
            const SparsityPatternType &nontrilinos_sparsity_pattern,
            const MPI_Comm            &communicator  = MPI_COMM_WORLD,
            const bool                 exchange_data = false);

    template <typename SparsityPatternType>
    void
    reinit (const IndexSet            &parallel_partitioning,
            const SparsityPatternType &nontrilinos_sparsity_pattern,
            const MPI_Comm            &communicator  = MPI_COMM_WORLD,
            const bool                 exchange_data = false);


    bool is_compressed () const;

    unsigned int max_entries_per_row () const;

    size_type n_rows () const;

    size_type n_cols () 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;

    size_type row_length (const size_type row) const;

    size_type bandwidth () const;

    bool empty () const;

    bool exists (const size_type i,
                 const size_type j) const;

    std::size_t memory_consumption () const;


    void add (const size_type i,
              const size_type j);


    template <typename ForwardIterator>
    void add_entries (const size_type  row,
                      ForwardIterator  begin,
                      ForwardIterator  end,
                      const bool       indices_are_sorted = false);


    const Epetra_FECrsGraph &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;

    DEAL_II_DEPRECATED
    const Epetra_Comm &trilinos_communicator () const;

    MPI_Comm get_mpi_communicator () const;


    IndexSet locally_owned_domain_indices() const;

    IndexSet locally_owned_range_indices() const;



    const_iterator begin () const;

    const_iterator end () const;

    const_iterator begin (const size_type r) const;

    const_iterator end (const size_type r) const;



    void write_ascii ();

    void print (std::ostream &out,
                const bool    write_extended_trilinos_info = false) const;

    void print_gnuplot (std::ostream &out) 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 sparsity patterns that "
                      "are the same object, but the operation requires that the "
                      "two objects are in fact different.");

    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.");
  private:

    std::unique_ptr<Epetra_Map> column_space_map;

    std::unique_ptr<Epetra_FECrsGraph> graph;

    std::unique_ptr<Epetra_CrsGraph> nonlocal_graph;

    friend class TrilinosWrappers::SparseMatrix;
    friend class SparsityPatternIterators::Accessor;
    friend class SparsityPatternIterators::Iterator;
  };



// -------------------------- inline and template functions ----------------------


#ifndef DOXYGEN

  namespace SparsityPatternIterators
  {

    inline
    Accessor::Accessor (const SparsityPattern *sp,
                        const size_type        row,
                        const size_type        index)
      :
      sparsity_pattern(const_cast<SparsityPattern *>(sp)),
      a_row(row),
      a_index(index)
    {
      visit_present_row ();
    }



    inline
    Accessor::size_type
    Accessor::row() const
    {
      Assert (a_row < sparsity_pattern->n_rows(), ExcBeyondEndOfSparsityPattern());
      return a_row;
    }



    inline
    Accessor::size_type
    Accessor::column() const
    {
      Assert (a_row < sparsity_pattern->n_rows(), ExcBeyondEndOfSparsityPattern());
      return (*colnum_cache)[a_index];
    }



    inline
    Accessor::size_type
    Accessor::index() const
    {
      Assert (a_row < sparsity_pattern->n_rows(), ExcBeyondEndOfSparsityPattern());
      return a_index;
    }



    inline
    Iterator::Iterator(const SparsityPattern *sp,
                       const size_type        row,
                       const size_type        index)
      :
      accessor(sp, row, index)
    {}


    inline
    Iterator::Iterator(const Iterator &) = default;



    inline
    Iterator &
    Iterator::operator++ ()
    {
      Assert (accessor.a_row < accessor.sparsity_pattern->n_rows(),
              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.sparsity_pattern->n_rows())
                 &&
                 (accessor.sparsity_pattern->row_length(accessor.a_row) == 0))
            ++accessor.a_row;

          accessor.visit_present_row();
        }
      return *this;
    }



    inline
    Iterator
    Iterator::operator++ (int)
    {
      const Iterator old_state = *this;
      ++(*this);
      return old_state;
    }



    inline
    const Accessor &
    Iterator::operator* () const
    {
      return accessor;
    }



    inline
    const Accessor *
    Iterator::operator-> () const
    {
      return &accessor;
    }



    inline
    bool
    Iterator::operator == (const Iterator &other) const
    {
      return (accessor.a_row == other.accessor.a_row &&
              accessor.a_index == other.accessor.a_index);
    }



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



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

  }



  inline
  SparsityPattern::const_iterator
  SparsityPattern::begin() const
  {
    return const_iterator(this, 0, 0);
  }



  inline
  SparsityPattern::const_iterator
  SparsityPattern::end() const
  {
    return const_iterator(this, n_rows(), 0);
  }



  inline
  SparsityPattern::const_iterator
  SparsityPattern::begin(const size_type r) const
  {
    Assert (r < n_rows(), ExcIndexRangeType<size_type>(r, 0, n_rows()));
    if (row_length(r) > 0)
      return const_iterator(this, r, 0);
    else
      return end (r);
  }



  inline
  SparsityPattern::const_iterator
  SparsityPattern::end(const size_type r) const
  {
    Assert (r < n_rows(), ExcIndexRangeType<size_type>(r, 0, n_rows()));

    // place the iterator on the first entry
    // past this line, or at the end of the
    // matrix
    for (size_type i=r+1; i<n_rows(); ++i)
      if (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
  bool
  SparsityPattern::in_local_range (const size_type index) const
  {
    TrilinosWrappers::types::int_type begin, end;
#ifndef DEAL_II_WITH_64BIT_INDICES
    begin = graph->RowMap().MinMyGID();
    end = graph->RowMap().MaxMyGID()+1;
#else
    begin = graph->RowMap().MinMyGID64();
    end = graph->RowMap().MaxMyGID64()+1;
#endif

    return ((index >= static_cast<size_type>(begin)) &&
            (index < static_cast<size_type>(end)));
  }



  inline
  bool
  SparsityPattern::is_compressed () const
  {
    return graph->Filled();
  }



  inline
  bool
  SparsityPattern::empty () const
  {
    return ((n_rows() == 0) && (n_cols() == 0));
  }



  inline
  void
  SparsityPattern::add (const size_type i,
                        const size_type j)
  {
    add_entries (i, &j, &j+1);
  }



  template <typename ForwardIterator>
  inline
  void
  SparsityPattern::add_entries (const size_type row,
                                ForwardIterator begin,
                                ForwardIterator end,
                                const bool      /*indices_are_sorted*/)
  {
    if (begin == end)
      return;

    // verify that the size of the data type Trilinos expects matches that the
    // iterator points to. we allow for some slippage between signed and
    // unsigned and only compare that they are both either 32 or 64 bit. to
    // write this test properly, not that we cannot compare the size of
    // '*begin' because 'begin' may be an iterator and '*begin' may be an
    // accessor class. consequently, we need to somehow get an actual value
    // from it which we can by evaluating an expression such as when
    // multiplying the value produced by 2
    Assert (sizeof(TrilinosWrappers::types::int_type) ==
            sizeof((*begin)*2),
            ExcNotImplemented());

    TrilinosWrappers::types::int_type *col_index_ptr =
      (TrilinosWrappers::types::int_type *)(&*begin);
    const int n_cols = static_cast<int>(end - begin);

    int ierr;
    if ( graph->RowMap().LID(static_cast<TrilinosWrappers::types::int_type>(row)) != -1)
      ierr = graph->InsertGlobalIndices (row, n_cols, col_index_ptr);
    else if (nonlocal_graph.get() != nullptr)
      {
        // this is the case when we have explicitly set the off-processor rows
        // and want to create a separate matrix object for them (to retain
        // thread-safety)
        Assert (nonlocal_graph->RowMap().LID(static_cast<TrilinosWrappers::types::int_type>(row)) != -1,
                ExcMessage("Attempted to write into off-processor matrix row "
                           "that has not be specified as being writable upon "
                           "initialization"));
        ierr = nonlocal_graph->InsertGlobalIndices (row, n_cols, col_index_ptr);
      }
    else
      ierr = graph->InsertGlobalIndices
             (1, (TrilinosWrappers::types::int_type *)&row, n_cols, col_index_ptr);

    AssertThrow (ierr >= 0, ExcTrilinosError(ierr));
  }



  inline
  const Epetra_FECrsGraph &
  SparsityPattern::trilinos_sparsity_pattern () const
  {
    return *graph;
  }



  inline
  IndexSet
  SparsityPattern::locally_owned_domain_indices () const
  {
    return IndexSet(static_cast<const Epetra_Map &>(graph->DomainMap()));
  }



  inline
  IndexSet
  SparsityPattern::locally_owned_range_indices () const
  {
    return IndexSet(static_cast<const Epetra_Map &>(graph->RangeMap()));
  }

#endif // DOXYGEN
}


DEAL_II_NAMESPACE_CLOSE


#endif // DEAL_II_WITH_TRILINOS


/*--------------------   trilinos_sparsity_pattern.h     --------------------*/

#endif
/*--------------------   trilinos_sparsity_pattern.h     --------------------*/