sparsity_pattern.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2000 - 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_sparsity_pattern_h
#define dealii_sparsity_pattern_h


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

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

// boost::serialization::make_array used to be in array.hpp, but was
// moved to a different file in BOOST 1.64
#include <boost/version.hpp>
#if BOOST_VERSION >= 106400
#  include <boost/serialization/array_wrapper.hpp>
#else
#  include <boost/serialization/array.hpp>
#endif
#include <boost/serialization/split_member.hpp>

#include <algorithm>
#include <iostream>
#include <vector>

DEAL_II_NAMESPACE_OPEN

class SparsityPattern;
class DynamicSparsityPattern;
class ChunkSparsityPattern;
template <typename number>
class FullMatrix;
template <typename number>
class SparseMatrix;
template <typename number>
class SparseLUDecomposition;
template <typename number>
class SparseILU;
template <typename VectorType>
class VectorSlice;

namespace ChunkSparsityPatternIterators
{
  class Accessor;
}


namespace internals
{
  namespace SparsityPatternTools
  {
    using size_type = types::global_dof_index;

    size_type
    get_column_index_from_iterator(const size_type i);

    template <typename value>
    size_type
    get_column_index_from_iterator(const std::pair<size_type, value> &i);

    template <typename value>
    size_type
    get_column_index_from_iterator(const std::pair<const size_type, value> &i);

  } // namespace SparsityPatternTools
} // namespace internals


namespace SparsityPatternIterators
{
  // forward declaration
  class Iterator;

  using size_type = types::global_dof_index;

  class Accessor
  {
  public:
    Accessor(const SparsityPattern *matrix,
             const std::size_t      index_within_sparsity);

    Accessor(const SparsityPattern *matrix);

    size_type
    row() const;

    size_type
    index() const;

    size_type
    global_index() const;

    size_type
    column() const;

    bool
    is_valid_entry() const;

    bool
    operator==(const Accessor &) const;

    bool
    operator<(const Accessor &) const;

  protected:
    const SparsityPattern *sparsity_pattern;

    std::size_t index_within_sparsity;

    void
    advance();

    friend class Iterator;

    friend class ChunkSparsityPatternIterators::Accessor;
  };



  class Iterator
  {
  public:
    Iterator(const SparsityPattern *sp,
             const std::size_t      index_within_sparsity);

    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;

    int
    operator-(const Iterator &p) const;

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



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

  using const_iterator = SparsityPatternIterators::Iterator;

  using iterator = SparsityPatternIterators::Iterator;


  static const size_type invalid_entry = numbers::invalid_size_type;

  // @{
  SparsityPattern();

  SparsityPattern(const SparsityPattern &);

  SparsityPattern(const size_type    m,
                  const size_type    n,
                  const unsigned int max_per_row);


  SparsityPattern(const size_type                  m,
                  const size_type                  n,
                  const std::vector<unsigned int> &row_lengths);

  SparsityPattern(const size_type m, const unsigned int max_per_row);

  SparsityPattern(const size_type                  m,
                  const std::vector<unsigned int> &row_lengths);

  SparsityPattern(const SparsityPattern &original,
                  const unsigned int     max_per_row,
                  const size_type        extra_off_diagonals);

  ~SparsityPattern() override = default;

  SparsityPattern &
  operator=(const SparsityPattern &);

  void
  reinit(const size_type m, const size_type n, const unsigned int max_per_row);


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


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

  void
  compress();


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

  void
  copy_from(const DynamicSparsityPattern &dsp);

  void
  copy_from(const SparsityPattern &sp);

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

  void
  symmetrize();

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

  // @}



  // @{

  iterator
  begin() const;

  iterator
  end() const;

  iterator
  begin(const size_type r) const;

  iterator
  end(const size_type r) const;


  // @}
  // @{
  bool
  operator==(const SparsityPattern &) const;

  bool
  empty() const;

  size_type
  max_entries_per_row() const;

  size_type
  bandwidth() const;

  std::size_t
  n_nonzero_elements() const;

  bool
  is_compressed() const;

  size_type
  n_rows() const;

  size_type
  n_cols() const;

  unsigned int
  row_length(const size_type row) const;

  bool
  stores_only_added_elements() const;

  std::size_t
  memory_consumption() const;

  // @}
  // @{
  size_type
  operator()(const size_type i, const size_type j) const;

  std::pair<size_type, size_type>
  matrix_position(const std::size_t global_index) const;

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

  size_type
  row_position(const size_type i, const size_type j) const;

  size_type
  column_number(const size_type row, const unsigned int index) 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;

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


  template <class Archive>
  void
  save(Archive &ar, const unsigned int version) const;

  template <class Archive>
  void
  load(Archive &ar, const unsigned int version);

  BOOST_SERIALIZATION_SPLIT_MEMBER()

  // @}

  
  DeclException2(ExcNotEnoughSpace,
                 int,
                 int,
                 << "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.");
  DeclException2(ExcIteratorRange,
                 int,
                 int,
                 << "The iterators denote a range of " << arg1
                 << " elements, but the given number of rows was " << arg2);
  DeclException1(ExcInvalidNumberOfPartitions,
                 int,
                 << "The number of partitions you gave is " << arg1
                 << ", but must be greater than zero.");
private:
  size_type max_dim;

  size_type rows;

  size_type cols;

  std::size_t max_vec_len;

  unsigned int max_row_length;

  std::unique_ptr<std::size_t[]> rowstart;

  std::unique_ptr<size_type[]> colnums;

  bool compressed;

  bool store_diagonal_first_in_row;

  template <typename number>
  friend class SparseMatrix;
  template <typename number>
  friend class SparseLUDecomposition;
  template <typename number>
  friend class SparseILU;
  template <typename number>
  friend class ChunkSparseMatrix;

  friend class ChunkSparsityPattern;

  friend class SparsityPatternIterators::Iterator;
  friend class SparsityPatternIterators::Accessor;
  friend class ChunkSparsityPatternIterators::Accessor;
};


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

#ifndef DOXYGEN


namespace SparsityPatternIterators
{
  inline Accessor::Accessor(const SparsityPattern *sparsity_pattern,
                            const std::size_t      i)
    : sparsity_pattern(sparsity_pattern)
    , index_within_sparsity(i)
  {}



  inline Accessor::Accessor(const SparsityPattern *sparsity_pattern)
    : sparsity_pattern(sparsity_pattern)
    , index_within_sparsity(sparsity_pattern->rowstart[sparsity_pattern->rows])
  {}



  inline bool
  Accessor::is_valid_entry() const
  {
    return (index_within_sparsity <
              sparsity_pattern->rowstart[sparsity_pattern->rows] &&
            sparsity_pattern->colnums[index_within_sparsity] !=
              SparsityPattern::invalid_entry);
  }



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

    const std::size_t *insert_point =
      std::upper_bound(sparsity_pattern->rowstart.get(),
                       sparsity_pattern->rowstart.get() +
                         sparsity_pattern->rows + 1,
                       index_within_sparsity);
    return insert_point - sparsity_pattern->rowstart.get() - 1;
  }



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

    return (sparsity_pattern->colnums[index_within_sparsity]);
  }



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

    return index_within_sparsity - sparsity_pattern->rowstart[row()];
  }



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

    return index_within_sparsity;
  }



  inline bool
  Accessor::operator==(const Accessor &other) const
  {
    return (sparsity_pattern == other.sparsity_pattern &&
            index_within_sparsity == other.index_within_sparsity);
  }



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

    return index_within_sparsity < other.index_within_sparsity;
  }


  inline void
  Accessor::advance()
  {
    Assert(index_within_sparsity <
             sparsity_pattern->rowstart[sparsity_pattern->rows],
           ExcIteratorPastEnd());
    ++index_within_sparsity;
  }



  inline Iterator::Iterator(const SparsityPattern *sparsity_pattern,
                            const std::size_t      i)
    : accessor(sparsity_pattern, i)
  {}



  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 !(*this == other);
  }



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



  inline int
  Iterator::operator-(const Iterator &other) const
  {
    Assert(accessor.sparsity_pattern == other.accessor.sparsity_pattern,
           ExcInternalError());

    return (*this)->index_within_sparsity - other->index_within_sparsity;
  }
} // namespace SparsityPatternIterators



inline SparsityPattern::iterator
SparsityPattern::begin() const
{
  return iterator(this, rowstart[0]);
}



inline SparsityPattern::iterator
SparsityPattern::end() const
{
  return iterator(this, rowstart[rows]);
}



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

  return iterator(this, rowstart[r]);
}



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

  return iterator(this, rowstart[r + 1]);
}



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



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



inline bool
SparsityPattern::is_compressed() const
{
  return compressed;
}



inline bool
SparsityPattern::stores_only_added_elements() const
{
  return (store_diagonal_first_in_row == false);
}



inline unsigned int
SparsityPattern::row_length(const size_type row) const
{
  Assert(row < rows, ExcIndexRangeType<size_type>(row, 0, rows));
  return rowstart[row + 1] - rowstart[row];
}



inline SparsityPattern::size_type
SparsityPattern::column_number(const size_type    row,
                               const unsigned int index) const
{
  Assert(row < rows, ExcIndexRangeType<size_type>(row, 0, rows));
  Assert(index < row_length(row), ExcIndexRange(index, 0, row_length(row)));

  return colnums[rowstart[row] + index];
}



inline std::size_t
SparsityPattern::n_nonzero_elements() const
{
  Assert(compressed, ExcNotCompressed());

  if ((rowstart != nullptr) && (colnums != nullptr))
    return rowstart[rows] - rowstart[0];
  else
    // the object is empty or has zero size
    return 0;
}



template <class Archive>
inline void
SparsityPattern::save(Archive &ar, const unsigned int) const
{
  // forward to serialization function in the base class.
  ar &static_cast<const Subscriptor &>(*this);

  ar &max_dim &rows &cols &max_vec_len &max_row_length &compressed
    &store_diagonal_first_in_row;

  ar &boost::serialization::make_array(rowstart.get(), max_dim + 1);
  ar &boost::serialization::make_array(colnums.get(), max_vec_len);
}



template <class Archive>
inline void
SparsityPattern::load(Archive &ar, const unsigned int)
{
  // forward to serialization function in the base class.
  ar &static_cast<Subscriptor &>(*this);

  ar &max_dim &rows &cols &max_vec_len &max_row_length &compressed
    &store_diagonal_first_in_row;

  rowstart = std_cxx14::make_unique<std::size_t[]>(max_dim + 1);
  colnums  = std_cxx14::make_unique<size_type[]>(max_vec_len);

  ar &boost::serialization::make_array(rowstart.get(), max_dim + 1);
  ar &boost::serialization::make_array(colnums.get(), max_vec_len);
}



inline bool
SparsityPattern::operator==(const SparsityPattern &sp2) const
{
  // it isn't quite necessary to compare *all* member variables. by only
  // comparing the essential ones, we can say that two sparsity patterns are
  // equal even if one is compressed and the other is not (in which case some
  // of the member variables are not yet set correctly)
  if (rows != sp2.rows || cols != sp2.cols || compressed != sp2.compressed ||
      store_diagonal_first_in_row != sp2.store_diagonal_first_in_row)
    return false;

  for (size_type i = 0; i < rows + 1; ++i)
    if (rowstart[i] != sp2.rowstart[i])
      return false;

  for (size_type i = 0; i < rowstart[rows]; ++i)
    if (colnums[i] != sp2.colnums[i])
      return false;

  return true;
}



namespace internal
{
  namespace SparsityPatternTools
  {
    using size_type = types::global_dof_index;

    inline size_type
    get_column_index_from_iterator(const size_type i)
    {
      return i;
    }



    template <typename value>
    inline size_type
    get_column_index_from_iterator(const std::pair<size_type, value> &i)
    {
      return i.first;
    }



    template <typename value>
    inline size_type
    get_column_index_from_iterator(const std::pair<const size_type, value> &i)
    {
      return i.first;
    }
  } // namespace SparsityPatternTools
} // namespace internal



template <typename ForwardIterator>
void
SparsityPattern::copy_from(const size_type       n_rows,
                           const size_type       n_cols,
                           const ForwardIterator begin,
                           const ForwardIterator end)
{
  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<unsigned int> 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);

  // now enter all the elements into the matrix. note that if the matrix is
  // quadratic, then we already have the diagonal element preallocated
  //
  // for use in the inner loop, we define an alias to the type of the inner
  // iterators
  size_type row = 0;
  using inner_iterator =
    typename std::iterator_traits<ForwardIterator>::value_type::const_iterator;
  for (ForwardIterator i = begin; i != end; ++i, ++row)
    {
      size_type *          cols = &colnums[rowstart[row]] + (is_square ? 1 : 0);
      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, ExcIndexRange(col, 0, n_cols));

          if ((col != row) || !is_square)
            *cols++ = col;
        }
    }

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


#endif // DOXYGEN

DEAL_II_NAMESPACE_CLOSE

#endif