trilinos_tpetra_sparsity_pattern.h

Go to the documentation of this file.

// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 2024 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Part of the source code is dual licensed under Apache-2.0 WITH
// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
// governing the source code and code contributions can be found in
// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
//
// ------------------------------------------------------------------------
 
#ifndef dealii_trilinos_tpetra_sparsity_pattern_h
#define dealii_trilinos_tpetra_sparsity_pattern_h
 
#include <deal.II/base/config.h>
 
#include "deal.II/base/types.h"
 
#include "deal.II/lac/trilinos_tpetra_types.h"
 
#ifdef DEAL_II_TRILINOS_WITH_TPETRA
 
#  include <deal.II/base/index_set.h>
#  include <deal.II/base/mpi_stub.h>
#  include <deal.II/base/subscriptor.h>
 
#  include <deal.II/lac/exceptions.h>
#  include <deal.II/lac/sparsity_pattern_base.h>
#  include <deal.II/lac/trilinos_tpetra_sparse_matrix.h>
 
#  include <Tpetra_CrsGraph.hpp>
 
#  include <cmath>
#  include <memory>
#  include <vector>
 
 
DEAL_II_NAMESPACE_OPEN
 
// forward declarations
#  ifndef DOXYGEN
class DynamicSparsityPattern;
 
namespace LinearAlgebra
{
  namespace TpetraWrappers
  {
    template <typename MemorySpace>
    class SparsityPattern;
 
    template <typename Number, typename MemorySpace>
    class SparseMatrix;
 
    namespace SparsityPatternIterators
    {
      template <typename MemorySpace>
      class Iterator;
    }
  } // namespace TpetraWrappers
} // namespace LinearAlgebra
#  endif
 
namespace LinearAlgebra
{
  namespace TpetraWrappers
  {
    namespace SparsityPatternIterators
    {
      template <typename MemorySpace = ::MemorySpace::Host>
      class Accessor
      {
      public:
        using size_type = ::types::global_dof_index;
 
        Accessor(const SparsityPattern<MemorySpace> *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:
        SparsityPattern<MemorySpace> *sparsity_pattern;
 
        size_type a_row;
 
        size_type a_index;
 
        std::shared_ptr<std::vector<::types::signed_global_dof_index>>
          colnum_cache;
 
        void
        visit_present_row();
 
        // Make enclosing class a friend.
        friend class Iterator<MemorySpace>;
      };
 
      template <typename MemorySpace = ::MemorySpace::Host>
      class Iterator
      {
      public:
        using size_type = types::global_dof_index;
 
        Iterator(const SparsityPattern<MemorySpace> *sparsity_pattern,
                 const size_type                     row,
                 const size_type                     index);
 
        Iterator(const Iterator<MemorySpace> &i);
 
        Iterator<MemorySpace> &
        operator++();
 
        Iterator
        operator++(int);
 
        const Accessor<MemorySpace> &
        operator*() const;
 
        const Accessor<MemorySpace> *
        operator->() const;
 
        bool
        operator==(const Iterator<MemorySpace> &) const;
 
        bool
        operator!=(const Iterator<MemorySpace> &) const;
 
        bool
        operator<(const Iterator<MemorySpace> &) const;
 
        DeclException2(ExcInvalidIndexWithinRow,
                       size_type,
                       size_type,
                       << "Attempt to access element " << arg2 << " of row "
                       << arg1 << " which doesn't have that many elements.");
 
      private:
        Accessor<MemorySpace> accessor;
 
        friend class TpetraWrappers::SparsityPattern<MemorySpace>;
      };
 
    } // namespace SparsityPatternIterators
 
 
    template <typename MemorySpace = ::MemorySpace::Host>
    class SparsityPattern : public SparsityPatternBase
    {
    public:
      using size_type = ::types::global_dof_index;
      using const_iterator = SparsityPatternIterators::Iterator<MemorySpace>;
 
      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<MemorySpace> &&other) noexcept;
 
      SparsityPattern(
        const SparsityPattern<MemorySpace> &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<MemorySpace> &input_sparsity_pattern);
 
      template <typename SparsityPatternType>
      void
      copy_from(const SparsityPatternType &nontrilinos_sparsity_pattern);
 
      SparsityPattern<MemorySpace> &
      operator=(const SparsityPattern<MemorySpace> &input_sparsity_pattern);
 
      void
      clear();
 
      void
      compress();
      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;
 
      unsigned int
      local_size() const;
 
      std::pair<size_type, size_type>
      local_range() const;
 
      bool
      in_local_range(const size_type index) const;
 
      std::uint64_t
      n_nonzero_elements() const;
 
      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);
 
      virtual void
      add_row_entries(
        const ::types::global_dof_index                  &row,
        const ArrayView<const ::types::global_dof_index> &columns,
        const bool indices_are_sorted = false) override;
 
      using SparsityPatternBase::add_entries;
 
      Teuchos::RCP<TpetraTypes::GraphType<MemorySpace>>
      trilinos_sparsity_pattern() const;
 
      Teuchos::RCP<const TpetraTypes::MapType<MemorySpace>>
      domain_partitioner() const;
 
      Teuchos::RCP<const TpetraTypes::MapType<MemorySpace>>
      range_partitioner() const;
 
      MPI_Comm
      get_mpi_communicator() const;
 
      Teuchos::RCP<const Teuchos::Comm<int>>
      get_teuchos_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
      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.");
 
      DeclException4(ExcAccessToNonLocalElement,
                     size_type,
                     size_type,
                     size_type,
                     size_type,
                     << "You tried to access element (" << arg1 << '/' << arg2
                     << ')'
                     << " of a distributed matrix, but only rows in range ["
                     << arg3 << ',' << arg4
                     << "] are stored locally and can be accessed.");
 
      DeclException2(ExcAccessToNonPresentElement,
                     size_type,
                     size_type,
                     << "You tried to access element (" << arg1 << '/' << arg2
                     << ')' << " of a sparse matrix, but it appears to not"
                     << " exist in the Trilinos sparsity pattern.");
    private:
      Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> column_space_map;
 
      Teuchos::RCP<TpetraTypes::GraphType<MemorySpace>> graph;
 
      Teuchos::RCP<TpetraTypes::GraphType<MemorySpace>> nonlocal_graph;
 
      // TODO: currently only for double
      friend class SparseMatrix<double, MemorySpace>;
      friend class SparsityPatternIterators::Accessor<MemorySpace>;
      friend class SparsityPatternIterators::Iterator<MemorySpace>;
    };
 
 
 
    // ---------------- inline and template functions -----------------
 
 
#  ifndef DOXYGEN
 
    namespace SparsityPatternIterators
    {
      template <typename MemorySpace>
      inline Accessor<MemorySpace>::Accessor(
        const SparsityPattern<MemorySpace> *sp,
        const size_type                     row,
        const size_type                     index)
        : sparsity_pattern(const_cast<SparsityPattern<MemorySpace> *>(sp))
        , a_row(row)
        , a_index(index)
      {
        visit_present_row();
      }
 
 
 
      template <typename MemorySpace>
      inline typename Accessor<MemorySpace>::size_type
      Accessor<MemorySpace>::row() const
      {
        Assert(a_row < sparsity_pattern->n_rows(),
               ExcBeyondEndOfSparsityPattern());
        return a_row;
      }
 
 
 
      template <typename MemorySpace>
      inline typename Accessor<MemorySpace>::size_type
      Accessor<MemorySpace>::column() const
      {
        Assert(a_row < sparsity_pattern->n_rows(),
               ExcBeyondEndOfSparsityPattern());
        return (*colnum_cache)[a_index];
      }
 
 
 
      template <typename MemorySpace>
      inline typename Accessor<MemorySpace>::size_type
      Accessor<MemorySpace>::index() const
      {
        Assert(a_row < sparsity_pattern->n_rows(),
               ExcBeyondEndOfSparsityPattern());
        return a_index;
      }
 
 
 
      template <typename MemorySpace>
      inline Iterator<MemorySpace>::Iterator(
        const SparsityPattern<MemorySpace> *sp,
        const size_type                     row,
        const size_type                     index)
        : accessor(sp, row, index)
      {}
 
 
 
      template <typename MemorySpace>
      inline Iterator<MemorySpace>::Iterator(const Iterator<MemorySpace> &) =
        default;
 
 
 
      template <typename MemorySpace>
      inline Iterator<MemorySpace> &
      Iterator<MemorySpace>::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 >=
            static_cast<::types::signed_global_dof_index>(
              accessor.colnum_cache->size()))
          {
            accessor.a_index = 0;
            ++accessor.a_row;
 
            while (accessor.a_row <
                   static_cast<::types::signed_global_dof_index>(
                     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;
      }
 
 
 
      template <typename MemorySpace>
      inline Iterator<MemorySpace>
      Iterator<MemorySpace>::operator++(int)
      {
        const Iterator<MemorySpace> old_state = *this;
        ++(*this);
        return old_state;
      }
 
 
 
      template <typename MemorySpace>
      inline const Accessor<MemorySpace> &
      Iterator<MemorySpace>::operator*() const
      {
        return accessor;
      }
 
 
 
      template <typename MemorySpace>
      inline const Accessor<MemorySpace> *
      Iterator<MemorySpace>::operator->() const
      {
        return &accessor;
      }
 
 
 
      template <typename MemorySpace>
      inline bool
      Iterator<MemorySpace>::operator==(
        const Iterator<MemorySpace> &other) const
      {
        return (accessor.a_row == other.accessor.a_row &&
                accessor.a_index == other.accessor.a_index);
      }
 
 
 
      template <typename MemorySpace>
      inline bool
      Iterator<MemorySpace>::operator!=(
        const Iterator<MemorySpace> &other) const
      {
        return !(*this == other);
      }
 
 
 
      template <typename MemorySpace>
      inline bool
      Iterator<MemorySpace>::operator<(const Iterator<MemorySpace> &other) const
      {
        return (accessor.row() < other.accessor.row() ||
                (accessor.row() == other.accessor.row() &&
                 accessor.index() < other.accessor.index()));
      }
 
    } // namespace SparsityPatternIterators
 
 
 
    template <typename MemorySpace>
    inline typename SparsityPattern<MemorySpace>::const_iterator
    SparsityPattern<MemorySpace>::begin() const
    {
      const size_type first_valid_row = this->local_range().first;
      return const_iterator(this, first_valid_row, 0);
    }
 
 
 
    template <typename MemorySpace>
    inline typename SparsityPattern<MemorySpace>::const_iterator
    SparsityPattern<MemorySpace>::end() const
    {
      return const_iterator(this, n_rows(), 0);
    }
 
 
 
    template <typename MemorySpace>
    inline typename SparsityPattern<MemorySpace>::const_iterator
    SparsityPattern<MemorySpace>::begin(const size_type r) const
    {
      AssertIndexRange(r, n_rows());
      if (row_length(r) > 0)
        return const_iterator(this, r, 0);
      else
        return end(r);
    }
 
 
 
    template <typename MemorySpace>
    inline typename SparsityPattern<MemorySpace>::const_iterator
    SparsityPattern<MemorySpace>::end(const size_type r) const
    {
      AssertIndexRange(r, 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();
    }
 
 
 
    template <typename MemorySpace>
    inline bool
    SparsityPattern<MemorySpace>::in_local_range(const size_type index) const
    {
      const TrilinosWrappers::types::int_type begin =
        graph->getRowMap()->getMinGlobalIndex();
      const TrilinosWrappers::types::int_type end =
        graph->getRowMap()->getMaxGlobalIndex() + 1;
 
      return ((index >= static_cast<size_type>(begin)) &&
              (index < static_cast<size_type>(end)));
    }
 
 
 
    template <typename MemorySpace>
    inline bool
    SparsityPattern<MemorySpace>::is_compressed() const
    {
      return graph->isFillComplete();
    }
 
 
 
    template <typename MemorySpace>
    inline bool
    SparsityPattern<MemorySpace>::empty() const
    {
      return ((n_rows() == 0) && (n_cols() == 0));
    }
 
 
 
    template <typename MemorySpace>
    inline void
    SparsityPattern<MemorySpace>::add(const size_type i, const size_type j)
    {
      add_entries(i, &j, &j + 1);
    }
 
 
 
    template <typename MemorySpace>
    template <typename ForwardIterator>
    inline void
    SparsityPattern<MemorySpace>::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());
 
      const TrilinosWrappers::types::int_type *col_index_ptr_begin =
        reinterpret_cast<TrilinosWrappers::types::int_type *>(
          const_cast<std::decay_t<decltype(*begin)> *>(&*begin));
 
      const TrilinosWrappers::types::int_type *col_index_ptr_end =
        reinterpret_cast<TrilinosWrappers::types::int_type *>(
          const_cast<std::decay_t<decltype(*end)> *>(&*end));
 
      // Check at least for the first index that the conversion actually works
      AssertDimension(*col_index_ptr_begin, *begin);
      AssertDimension(*col_index_ptr_end, *end);
      TrilinosWrappers::types::int_type trilinos_row_index = row;
 
      // TODO: The following line creates an array by copying the entries.
      //       Perhaps there is a way to only create a 'view' of these arrays
      //       and pass that to Tpetra?
      Teuchos::Array<TrilinosWrappers::types::int_type> array(
        col_index_ptr_begin, col_index_ptr_end);
 
      if (row_is_stored_locally(row))
        graph->insertGlobalIndices(trilinos_row_index, array());
      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->getRowMap()->getLocalElement(row) !=
                   Teuchos::OrdinalTraits<
                     ::types::signed_global_dof_index>::invalid(),
                 ExcMessage("Attempted to write into off-processor matrix row "
                            "that has not be specified as being writable upon "
                            "initialization"));
          nonlocal_graph->insertGlobalIndices(trilinos_row_index, array);
        }
      else
        graph->insertGlobalIndices(trilinos_row_index, array);
    }
 
 
 
    template <typename MemorySpace>
    inline Teuchos::RCP<Tpetra::CrsGraph<int,
                                         ::types::signed_global_dof_index,
                                         TpetraTypes::NodeType<MemorySpace>>>
    SparsityPattern<MemorySpace>::trilinos_sparsity_pattern() const
    {
      return graph;
    }
 
 
 
    template <typename MemorySpace>
    inline IndexSet
    SparsityPattern<MemorySpace>::locally_owned_domain_indices() const
    {
      return IndexSet(graph->getDomainMap().getConst());
    }
 
 
 
    template <typename MemorySpace>
    inline IndexSet
    SparsityPattern<MemorySpace>::locally_owned_range_indices() const
    {
      return IndexSet(graph->getRangeMap().getConst());
    }
 
#  endif // DOXYGEN
  }      // namespace TpetraWrappers
 
} // namespace LinearAlgebra
 
 
DEAL_II_NAMESPACE_CLOSE
 
 
#endif // DEAL_II_TRILINOS_WITH_TPETRA
 
#endif