sparsity_pattern.h

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


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

#include <deal.II/base/array_view.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/linear_index_iterator.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 SparsityPatternBase;
class DynamicSparsityPattern;
class ChunkSparsityPattern;
template <typename number>
class FullMatrix;
template <typename number>
class SparseMatrix;
template <typename number>
class SparseLUDecomposition;
template <typename number>
class SparseILU;

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

    Accessor(const SparsityPatternBase *matrix, const std::size_t linear_index);

    Accessor(const SparsityPatternBase *matrix);

    Accessor();

    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:
    DeclExceptionMsg(DummyAccessor,
                     "The instance of this class was initialized"
                     " without SparsityPattern object, which"
                     " means that it is a dummy accessor that can"
                     " not do any operations.");

    const SparsityPatternBase *container;

    std::size_t linear_index;

    void
    advance();

    friend class LinearIndexIterator<Iterator, Accessor>;

    friend class ChunkSparsityPatternIterators::Accessor;
  };



  class Iterator : public LinearIndexIterator<Iterator, Accessor>
  {
  public:
    using size_type = types::global_dof_index;

    using container_pointer_type = SparsityPatternBase *;

    Iterator(const SparsityPatternBase *sp, const std::size_t linear_index);

    Iterator(const Accessor &accessor);
  };
} // namespace SparsityPatternIterators



class SparsityPatternBase : 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;

  // @{
  SparsityPatternBase();

  ~SparsityPatternBase() override = default;

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

  virtual void
  reinit(const size_type                      m,
         const size_type                      n,
         const ArrayView<const unsigned int> &row_lengths) = 0;

  void
  symmetrize();

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

  // @}

  // @{

  iterator
  begin() const;

  iterator
  end() const;

  iterator
  begin(const size_type r) const;

  iterator
  end(const size_type r) const;


  // @}

  // @{

  bool
  operator==(const SparsityPatternBase &) const;

  bool
  empty() const;

  bool
  exists(const size_type i, const size_type j) 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;

  std::size_t
  memory_consumption() const;

  // @}

  // @{

  size_type
  column_number(const size_type row, const unsigned int index) const;

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

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

  // @}

  // @{

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

  // @}

  
  DeclExceptionMsg(
    ExcNotCompressed,
    "The operation you attempted is only allowed after the SparsityPattern "
    "has been set up and compress() was called.");

  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(
    ExcMatrixIsCompressed,
    "The operation you attempted changes the structure of the SparsityPattern "
    "and is not possible after compress() has been called.");

  // @}


protected:
  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;

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

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

class SparsityPattern : public SparsityPatternBase
{
public:
  using size_type = SparsityPatternBase::size_type;

  using const_iterator = SparsityPatternBase::const_iterator;

  using iterator = SparsityPatternBase::iterator;

  using SparsityPatternBase::reinit;

  // @{
  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 &);

  virtual void
  reinit(const size_type                      m,
         const size_type                      n,
         const ArrayView<const unsigned int> &row_lengths) override;

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

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

  // @}


  // @{
  bool
  operator==(const SparsityPattern &) 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;

  // @}
  // @{

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

  void
  block_read(std::istream &in);

  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(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:
  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 DynamicSparsityPattern;

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


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

#ifndef DOXYGEN


namespace SparsityPatternIterators
{
  inline Accessor::Accessor(const SparsityPatternBase *sparsity_pattern,
                            const std::size_t          i)
    : container(sparsity_pattern)
    , linear_index(i)
  {}



  inline Accessor::Accessor(const SparsityPatternBase *sparsity_pattern)
    : container(sparsity_pattern)
    , linear_index(container->rowstart[container->rows])
  {}



  inline Accessor::Accessor()
    : container(nullptr)
    , linear_index(numbers::invalid_size_type)
  {}



  inline bool
  Accessor::is_valid_entry() const
  {
    Assert(container != nullptr, DummyAccessor());
    return (linear_index < container->rowstart[container->rows] &&
            container->colnums[linear_index] != SparsityPattern::invalid_entry);
  }



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

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



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

    return (container->colnums[linear_index]);
  }



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

    return linear_index - container->rowstart[row()];
  }



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

    return linear_index;
  }



  inline bool
  Accessor::operator==(const Accessor &other) const
  {
    Assert(container != nullptr, DummyAccessor());
    Assert(other.container != nullptr, DummyAccessor());
    return (container == other.container && linear_index == other.linear_index);
  }



  inline bool
  Accessor::operator<(const Accessor &other) const
  {
    Assert(container != nullptr, DummyAccessor());
    Assert(other.container != nullptr, DummyAccessor());
    Assert(container == other.container, ExcInternalError());

    return linear_index < other.linear_index;
  }



  inline void
  Accessor::advance()
  {
    Assert(container != nullptr, DummyAccessor());
    Assert(linear_index < container->rowstart[container->rows],
           ExcIteratorPastEnd());
    ++linear_index;
  }


  inline Iterator::Iterator(const SparsityPatternBase *sp,
                            const std::size_t          linear_index)
    : LinearIndexIterator<Iterator, Accessor>(Accessor(sp, linear_index))
  {}


  inline Iterator::Iterator(const Accessor &accessor)
    : LinearIndexIterator<Iterator, Accessor>(accessor)
  {}


} // namespace SparsityPatternIterators



inline SparsityPatternBase::iterator
SparsityPatternBase::begin() const
{
  if (n_rows() > 0)
    return {this, rowstart[0]};
  else
    return end();
}



inline SparsityPatternBase::iterator
SparsityPatternBase::end() const
{
  if (n_rows() > 0)
    return {this, rowstart[rows]};
  else
    return {nullptr, 0};
}



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

  return {this, rowstart[r]};
}



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

  return {this, rowstart[r + 1]};
}



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



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



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



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



inline unsigned int
SparsityPatternBase::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
SparsityPatternBase::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
SparsityPatternBase::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
SparsityPatternBase::save(Archive &ar, const unsigned int) const
{
  // forward to serialization function in the base class.
  ar &boost::serialization::base_object<const Subscriptor>(*this);

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

  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
SparsityPatternBase::load(Archive &ar, const unsigned int)
{
  // forward to serialization function in the base class.
  ar &boost::serialization::base_object<Subscriptor>(*this);

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

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



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



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



inline bool
SparsityPatternBase::operator==(const SparsityPatternBase &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)
    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;
}



inline bool
SparsityPattern::operator==(const SparsityPattern &sp2) const
{
  return (static_cast<const SparsityPatternBase &>(*this) == sp2) &&
         (store_diagonal_first_in_row == sp2.store_diagonal_first_in_row);
}



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