chunk_sparsity_pattern.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2008 - 2019 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_chunk_sparsity_pattern_h
#define dealii_chunk_sparsity_pattern_h


#include <deal.II/base/config.h>

#include <deal.II/base/exceptions.h>
#include <deal.II/base/subscriptor.h>

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

#include <iostream>
#include <vector>

DEAL_II_NAMESPACE_OPEN


template <typename>
class ChunkSparseMatrix;


namespace ChunkSparsityPatternIterators
{
  // forward declaration
  class Iterator;

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

    Accessor(const ChunkSparsityPattern *matrix, const size_type row);

    Accessor(const ChunkSparsityPattern *matrix);

    size_type
    row() const;

    std::size_t
    reduced_index() const;

    size_type
    column() const;

    bool
    is_valid_entry() const;


    bool
    operator==(const Accessor &) const;


    bool
    operator<(const Accessor &) const;

  protected:
    const ChunkSparsityPattern *sparsity_pattern;

    SparsityPatternIterators::Accessor reduced_accessor;

    size_type chunk_row;

    size_type chunk_col;

    void
    advance();

    friend class Iterator;
  };



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

    Iterator(const ChunkSparsityPattern *sp, const size_type row);

    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;

  private:
    Accessor accessor;
  };
} // namespace ChunkSparsityPatternIterators



class ChunkSparsityPattern : public Subscriptor
{
public:
  using size_type = types::global_dof_index;
  using const_iterator = ChunkSparsityPatternIterators::Iterator;

  using iterator = ChunkSparsityPatternIterators::Iterator;

  static const size_type invalid_entry = SparsityPattern::invalid_entry;

  ChunkSparsityPattern();

  ChunkSparsityPattern(const ChunkSparsityPattern &);

  ChunkSparsityPattern(const size_type m,
                       const size_type n,
                       const size_type max_chunks_per_row,
                       const size_type chunk_size);

  ChunkSparsityPattern(const size_type               m,
                       const size_type               n,
                       const std::vector<size_type> &row_lengths,
                       const size_type               chunk_size);

  ChunkSparsityPattern(const size_type n,
                       const size_type max_per_row,
                       const size_type chunk_size);

  ChunkSparsityPattern(const size_type               m,
                       const std::vector<size_type> &row_lengths,
                       const size_type               chunk_size);

  ~ChunkSparsityPattern() override = default;

  ChunkSparsityPattern &
  operator=(const ChunkSparsityPattern &);

  void
  reinit(const size_type m,
         const size_type n,
         const size_type max_per_row,
         const size_type chunk_size);

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

  void
  reinit(const size_type                   m,
         const size_type                   n,
         const ArrayView<const size_type> &row_lengths,
         const size_type                   chunk_size);

  void
  compress();

  template <typename ForwardIterator>
  void
  copy_from(const size_type       n_rows,
            const size_type       n_cols,
            const ForwardIterator begin,
            const ForwardIterator end,
            const size_type       chunk_size);

  template <typename SparsityPatternType>
  void
  copy_from(const SparsityPatternType &dsp, const size_type chunk_size);

  template <typename number>
  void
  copy_from(const FullMatrix<number> &matrix, const size_type chunk_size);

  template <typename Sparsity>
  void
  create_from(const size_type m,
              const size_type n,
              const Sparsity &sparsity_pattern_for_chunks,
              const size_type chunk_size,
              const bool      optimize_diagonal = true);

  bool
  empty() const;

  size_type
  get_chunk_size() const;

  size_type
  max_entries_per_row() const;

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

  void
  symmetrize();

  inline size_type
  n_rows() const;

  inline size_type
  n_cols() const;

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

  size_type
  row_length(const size_type row) const;

  size_type
  bandwidth() const;

  size_type
  n_nonzero_elements() const;

  bool
  is_compressed() const;

  bool
  stores_only_added_elements() const;

  iterator
  begin() const;

  iterator
  end() const;

  iterator
  begin(const size_type r) const;

  iterator
  end(const size_type r) const;

  void
  block_write(std::ostream &out) const;

  void
  block_read(std::istream &in);

  void
  print(std::ostream &out) const;

  void
  print_gnuplot(std::ostream &out) const;

  std::size_t
  memory_consumption() const;

  DeclException1(ExcInvalidNumber,
                 size_type,
                 << "The provided number is invalid here: " << arg1);
  DeclException2(ExcInvalidIndex,
                 size_type,
                 size_type,
                 << "The given index " << arg1 << " should be less than "
                 << arg2 << ".");
  DeclException2(ExcNotEnoughSpace,
                 size_type,
                 size_type,
                 << "Upon entering a new entry to row " << arg1
                 << ": there was no free entry any more. " << std::endl
                 << "(Maximum number of entries for this row: " << arg2
                 << "; maybe the matrix is already compressed?)");
  DeclExceptionMsg(
    ExcNotCompressed,
    "The operation you attempted is only allowed after the SparsityPattern "
    "has been set up and compress() was called.");
  DeclExceptionMsg(
    ExcMatrixIsCompressed,
    "The operation you attempted changes the structure of the SparsityPattern "
    "and is not possible after compress() has been called.");
  DeclException0(ExcEmptyObject);
  DeclException2(ExcIteratorRange,
                 size_type,
                 size_type,
                 << "The iterators denote a range of " << arg1
                 << " elements, but the given number of rows was " << arg2);
  DeclException0(ExcMETISNotInstalled);
  DeclException1(ExcInvalidNumberOfPartitions,
                 size_type,
                 << "The number of partitions you gave is " << arg1
                 << ", but must be greater than zero.");
  DeclException2(ExcInvalidArraySize,
                 size_type,
                 size_type,
                 << "The array has size " << arg1 << " but should have size "
                 << arg2);
private:
  size_type rows;

  size_type cols;

  size_type chunk_size;

  SparsityPattern sparsity_pattern;

  template <typename>
  friend class ChunkSparseMatrix;

  friend class ChunkSparsityPatternIterators::Accessor;
};


/*---------------------- Inline functions -----------------------------------*/

#ifndef DOXYGEN

namespace ChunkSparsityPatternIterators
{
  inline Accessor::Accessor(const ChunkSparsityPattern *sparsity_pattern,
                            const size_type             row)
    : sparsity_pattern(sparsity_pattern)
    , reduced_accessor(row == sparsity_pattern->n_rows() ?
                         *sparsity_pattern->sparsity_pattern.end() :
                         *sparsity_pattern->sparsity_pattern.begin(
                           row / sparsity_pattern->get_chunk_size()))
    , chunk_row(row == sparsity_pattern->n_rows() ?
                  0 :
                  row % sparsity_pattern->get_chunk_size())
    , chunk_col(0)
  {}



  inline Accessor::Accessor(const ChunkSparsityPattern *sparsity_pattern)
    : sparsity_pattern(sparsity_pattern)
    , reduced_accessor(*sparsity_pattern->sparsity_pattern.end())
    , chunk_row(0)
    , chunk_col(0)
  {}



  inline bool
  Accessor::is_valid_entry() const
  {
    return reduced_accessor.is_valid_entry() &&
           sparsity_pattern->get_chunk_size() * reduced_accessor.row() +
               chunk_row <
             sparsity_pattern->n_rows() &&
           sparsity_pattern->get_chunk_size() * reduced_accessor.column() +
               chunk_col <
             sparsity_pattern->n_cols();
  }



  inline Accessor::size_type
  Accessor::row() const
  {
    Assert(is_valid_entry() == true, ExcInvalidIterator());

    return sparsity_pattern->get_chunk_size() * reduced_accessor.row() +
           chunk_row;
  }



  inline Accessor::size_type
  Accessor::column() const
  {
    Assert(is_valid_entry() == true, ExcInvalidIterator());

    return sparsity_pattern->get_chunk_size() * reduced_accessor.column() +
           chunk_col;
  }



  inline std::size_t
  Accessor::reduced_index() const
  {
    Assert(is_valid_entry() == true, ExcInvalidIterator());

    return reduced_accessor.linear_index;
  }



  inline bool
  Accessor::operator==(const Accessor &other) const
  {
    // no need to check for equality of sparsity patterns as this is done in
    // the reduced case already and every ChunkSparsityPattern has its own
    // reduced sparsity pattern
    return (reduced_accessor == other.reduced_accessor &&
            chunk_row == other.chunk_row && chunk_col == other.chunk_col);
  }



  inline bool
  Accessor::operator<(const Accessor &other) const
  {
    Assert(sparsity_pattern == other.sparsity_pattern, ExcInternalError());

    if (chunk_row != other.chunk_row)
      {
        if (reduced_accessor.linear_index ==
            reduced_accessor.container->n_nonzero_elements())
          return false;
        if (other.reduced_accessor.linear_index ==
            reduced_accessor.container->n_nonzero_elements())
          return true;

        const auto global_row = sparsity_pattern->get_chunk_size() *
                                  reduced_accessor.row() +
                                chunk_row,
                   other_global_row = sparsity_pattern->get_chunk_size() *
                                        other.reduced_accessor.row() +
                                      other.chunk_row;
        if (global_row < other_global_row)
          return true;
        else if (global_row > other_global_row)
          return false;
      }

    return (
      reduced_accessor.linear_index < other.reduced_accessor.linear_index ||
      (reduced_accessor.linear_index == other.reduced_accessor.linear_index &&
       chunk_col < other.chunk_col));
  }


  inline void
  Accessor::advance()
  {
    const auto chunk_size = sparsity_pattern->get_chunk_size();
    Assert(chunk_row < chunk_size && chunk_col < chunk_size,
           ExcIteratorPastEnd());
    Assert(reduced_accessor.row() * chunk_size + chunk_row <
               sparsity_pattern->n_rows() &&
             reduced_accessor.column() * chunk_size + chunk_col <
               sparsity_pattern->n_cols(),
           ExcIteratorPastEnd());
    if (chunk_size == 1)
      {
        reduced_accessor.advance();
        return;
      }

    ++chunk_col;

    // end of chunk
    if (chunk_col == chunk_size ||
        reduced_accessor.column() * chunk_size + chunk_col ==
          sparsity_pattern->n_cols())
      {
        const auto reduced_row = reduced_accessor.row();
        // end of row
        if (reduced_accessor.linear_index + 1 ==
            reduced_accessor.container->rowstart[reduced_row + 1])
          {
            ++chunk_row;

            chunk_col = 0;

            // end of chunk rows or end of matrix
            if (chunk_row == chunk_size ||
                (reduced_row * chunk_size + chunk_row ==
                 sparsity_pattern->n_rows()))
              {
                chunk_row = 0;
                reduced_accessor.advance();
              }
            // go back to the beginning of the same reduced row but with
            // chunk_row increased by one
            else
              reduced_accessor.linear_index =
                reduced_accessor.container->rowstart[reduced_row];
          }
        // advance within chunk
        else
          {
            reduced_accessor.advance();
            chunk_col = 0;
          }
      }
  }



  inline Iterator::Iterator(const ChunkSparsityPattern *sparsity_pattern,
                            const size_type             row)
    : accessor(sparsity_pattern, row)
  {}



  inline Iterator &
  Iterator::operator++()
  {
    accessor.advance();
    return *this;
  }



  inline Iterator
  Iterator::operator++(int)
  {
    const Iterator iter = *this;
    accessor.advance();
    return iter;
  }



  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 == other.accessor);
  }



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


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

} // namespace ChunkSparsityPatternIterators



inline ChunkSparsityPattern::iterator
ChunkSparsityPattern::begin() const
{
  return {this, 0};
}


inline ChunkSparsityPattern::iterator
ChunkSparsityPattern::end() const
{
  return {this, n_rows()};
}



inline ChunkSparsityPattern::iterator
ChunkSparsityPattern::begin(const size_type r) const
{
  Assert(r < n_rows(), ExcIndexRange(r, 0, n_rows()));
  return {this, r};
}



inline ChunkSparsityPattern::iterator
ChunkSparsityPattern::end(const size_type r) const
{
  Assert(r < n_rows(), ExcIndexRange(r, 0, n_rows()));
  return {this, r + 1};
}



inline ChunkSparsityPattern::size_type
ChunkSparsityPattern::n_rows() const
{
  return rows;
}


inline ChunkSparsityPattern::size_type
ChunkSparsityPattern::n_cols() const
{
  return cols;
}



inline ChunkSparsityPattern::size_type
ChunkSparsityPattern::get_chunk_size() const
{
  return chunk_size;
}



inline bool
ChunkSparsityPattern::is_compressed() const
{
  return sparsity_pattern.compressed;
}



template <typename ForwardIterator>
void
ChunkSparsityPattern::copy_from(const size_type       n_rows,
                                const size_type       n_cols,
                                const ForwardIterator begin,
                                const ForwardIterator end,
                                const size_type       chunk_size)
{
  Assert(static_cast<size_type>(std::distance(begin, end)) == n_rows,
         ExcIteratorRange(std::distance(begin, end), n_rows));

  // first determine row lengths for each row. if the matrix is quadratic,
  // then we might have to add an additional entry for the diagonal, if that
  // is not yet present. as we have to call compress anyway later on, don't
  // bother to check whether that diagonal entry is in a certain row or not
  const bool             is_square = (n_rows == n_cols);
  std::vector<size_type> row_lengths;
  row_lengths.reserve(n_rows);
  for (ForwardIterator i = begin; i != end; ++i)
    row_lengths.push_back(std::distance(i->begin(), i->end()) +
                          (is_square ? 1 : 0));
  reinit(n_rows, n_cols, row_lengths, chunk_size);

  // now enter all the elements into the matrix
  size_type row = 0;
  using inner_iterator =
    typename std::iterator_traits<ForwardIterator>::value_type::const_iterator;
  for (ForwardIterator i = begin; i != end; ++i, ++row)
    {
      const inner_iterator end_of_row = i->end();
      for (inner_iterator j = i->begin(); j != end_of_row; ++j)
        {
          const size_type col =
            internal::SparsityPatternTools::get_column_index_from_iterator(*j);
          Assert(col < n_cols, ExcInvalidIndex(col, n_cols));

          add(row, col);
        }
    }

  // finally compress everything. this also sorts the entries within each row
  compress();
}


#endif // DOXYGEN

DEAL_II_NAMESPACE_CLOSE

#endif