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.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------

#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/index_set.h>
#    include <deal.II/base/subscriptor.h>

#    include <deal.II/lac/exceptions.h>

#    include <Epetra_FECrsGraph.h>
#    include <Epetra_Map.h>

#    include <cmath>
#    include <memory>
#    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
class SparsityPattern;
class DynamicSparsityPattern;

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

  namespace SparsityPatternIterators
  {
    // forward declaration
    class Iterator;

    class Accessor
    {
    public:
      using size_type = ::types::global_dof_index;

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

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

  } // namespace SparsityPatternIterators


  class SparsityPattern : public Subscriptor
  {
  public:
    using size_type = ::types::global_dof_index;

    using const_iterator = SparsityPatternIterators::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() override = 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;

    bool
    row_is_stored_locally(const size_type i) 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 that is stored locally.
      if (accessor.a_index >= accessor.colnum_cache->size())
        {
          accessor.a_index = 0;
          ++accessor.a_row;

          while (accessor.a_row < accessor.sparsity_pattern->n_rows())
            {
              const auto row_length =
                accessor.sparsity_pattern->row_length(accessor.a_row);
              if (row_length == 0 ||
                  !accessor.sparsity_pattern->row_is_stored_locally(
                    accessor.a_row))
                ++accessor.a_row;
              else
                break;
            }

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

  } // namespace SparsityPatternIterators



  inline SparsityPattern::const_iterator
  SparsityPattern::begin() const
  {
    const size_type first_valid_row = this->local_range().first;
    return const_iterator(this, first_valid_row, 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 =
      reinterpret_cast<TrilinosWrappers::types::int_type *>(
        const_cast<typename std::decay<decltype(*begin)>::type *>(&*begin));
    // Check at least for the first index that the conversion actually works
    AssertDimension(*col_index_ptr, *begin);
    TrilinosWrappers::types::int_type trilinos_row_index = row;
    const int                         n_cols = static_cast<int>(end - begin);

    int ierr;
    if (row_is_stored_locally(row))
      ierr =
        graph->InsertGlobalIndices(trilinos_row_index, 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(trilinos_row_index,
                                                   n_cols,
                                                   col_index_ptr);
      }
    else
      ierr = graph->InsertGlobalIndices(1,
                                        &trilinos_row_index,
                                        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(graph->DomainMap());
  }



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

#    endif // DOXYGEN
} // namespace TrilinosWrappers


DEAL_II_NAMESPACE_CLOSE


#  endif // DEAL_II_WITH_TRILINOS


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

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