dynamic_sparsity_pattern.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2011 - 2017 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_dynamic_sparsity_pattern_h
#define dealii_dynamic_sparsity_pattern_h


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

#include <deal.II/base/index_set.h>
#include <deal.II/base/subscriptor.h>
#include <deal.II/base/utilities.h>

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

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

DEAL_II_NAMESPACE_OPEN

class DynamicSparsityPattern;


namespace DynamicSparsityPatternIterators
{
  // forward declaration
  class Iterator;

  using size_type = types::global_dof_index;

  class Accessor
  {
  public:
    Accessor(const DynamicSparsityPattern *sparsity_pattern,
             const size_type               row,
             const unsigned int            index_within_row);

    Accessor(const DynamicSparsityPattern *sparsity_pattern);

    size_type
    row() const;

    size_type
    index() const;

    size_type
    column() const;

    bool
    operator==(const Accessor &) const;

    bool
    operator<(const Accessor &) const;

  protected:
    const DynamicSparsityPattern *sparsity_pattern;

    size_type current_row;

    std::vector<size_type>::const_iterator current_entry;

    std::vector<size_type>::const_iterator end_of_row;

    void
    advance();

    friend class Iterator;
  };



  class Iterator
  {
  public:
    Iterator(const DynamicSparsityPattern *sp,
             const size_type               row,
             const unsigned int            index_within_row);

    Iterator(const DynamicSparsityPattern *sp);

    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 DynamicSparsityPatternIterators


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

  using iterator = DynamicSparsityPatternIterators::Iterator;

  using const_iterator = DynamicSparsityPatternIterators::Iterator;

  DynamicSparsityPattern();

  DynamicSparsityPattern(const DynamicSparsityPattern &);

  DynamicSparsityPattern(const size_type m,
                         const size_type n,
                         const IndexSet &rowset = IndexSet());

  DynamicSparsityPattern(const IndexSet &indexset);

  DynamicSparsityPattern(const size_type n);

  DynamicSparsityPattern &
  operator=(const DynamicSparsityPattern &);

  void
  reinit(const size_type m,
         const size_type n,
         const IndexSet &rowset = IndexSet());

  void
  compress();

  bool
  empty() const;

  size_type
  max_entries_per_row() 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_unique_and_sorted = false);

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

  void
  symmetrize();

  template <typename SparsityPatternTypeLeft, typename SparsityPatternTypeRight>
  void
  compute_mmult_pattern(const SparsityPatternTypeLeft & left,
                        const SparsityPatternTypeRight &right);

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

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

  size_type
  n_rows() const;

  size_type
  n_cols() const;

  size_type
  row_length(const size_type row) const;

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

  // @{

  iterator
  begin() const;

  iterator
  end() const;

  iterator
  begin(const size_type r) const;

  iterator
  end(const size_type r) const;

  // @}

  size_type
  bandwidth() const;

  size_type
  n_nonzero_elements() const;

  const IndexSet &
  row_index_set() const;

  static bool
  stores_only_added_elements();

  size_type
  memory_consumption() const;

private:
  bool have_entries;

  size_type rows;

  size_type cols;

  IndexSet rowset;


  struct Line
  {
  public:
    std::vector<size_type> entries;

    void
    add(const size_type col_num);

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

    size_type
    memory_consumption() const;
  };


  std::vector<Line> lines;

  // make the accessor class a friend
  friend class DynamicSparsityPatternIterators::Accessor;
};

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


namespace DynamicSparsityPatternIterators
{
  inline Accessor::Accessor(const DynamicSparsityPattern *sparsity_pattern,
                            const size_type               row,
                            const unsigned int            index_within_row)
    : sparsity_pattern(sparsity_pattern)
    , current_row(row)
    , current_entry(
        ((sparsity_pattern->rowset.size() == 0) ?
           sparsity_pattern->lines[current_row].entries.begin() :
           sparsity_pattern
             ->lines[sparsity_pattern->rowset.index_within_set(current_row)]
             .entries.begin()) +
        index_within_row)
    , end_of_row(
        (sparsity_pattern->rowset.size() == 0) ?
          sparsity_pattern->lines[current_row].entries.end() :
          sparsity_pattern
            ->lines[sparsity_pattern->rowset.index_within_set(current_row)]
            .entries.end())
  {
    AssertIndexRange(current_row, sparsity_pattern->n_rows());
    Assert((sparsity_pattern->rowset.size() == 0) ||
             sparsity_pattern->rowset.is_element(current_row),
           ExcMessage("You can't create an iterator into a "
                      "DynamicSparsityPattern's row that is not "
                      "actually stored by that sparsity pattern "
                      "based on the IndexSet argument to it."));
    AssertIndexRange(
      index_within_row,
      ((sparsity_pattern->rowset.size() == 0) ?
         sparsity_pattern->lines[current_row].entries.size() :
         sparsity_pattern
           ->lines[sparsity_pattern->rowset.index_within_set(current_row)]
           .entries.size()));
  }


  inline Accessor::Accessor(const DynamicSparsityPattern *sparsity_pattern)
    : sparsity_pattern(sparsity_pattern)
    , current_row(numbers::invalid_size_type)
    , current_entry()
    , end_of_row()
  {}



  inline size_type
  Accessor::row() const
  {
    Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());

    return current_row;
  }


  inline size_type
  Accessor::column() const
  {
    Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());

    return *current_entry;
  }


  inline size_type
  Accessor::index() const
  {
    Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());

    return (current_entry -
            ((sparsity_pattern->rowset.size() == 0) ?
               sparsity_pattern->lines[current_row].entries.begin() :
               sparsity_pattern
                 ->lines[sparsity_pattern->rowset.index_within_set(current_row)]
                 .entries.begin()));
  }



  inline bool
  Accessor::operator==(const Accessor &other) const
  {
    // compare the sparsity pattern the iterator points into, the
    // current row, and the location within this row. ignore the
    // latter if the row is past-the-end because in that case the
    // current_entry field may not point to a deterministic location
    return (sparsity_pattern == other.sparsity_pattern &&
            current_row == other.current_row &&
            ((current_row == numbers::invalid_size_type) ||
             (current_entry == other.current_entry)));
  }



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

    // if *this is past-the-end, then it is less than no one
    if (current_row == numbers::invalid_size_type)
      return (false);
    // now *this should be an valid value
    Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());

    // if other is past-the-end
    if (other.current_row == numbers::invalid_size_type)
      return (true);
    // now other should be an valid value
    Assert(other.current_row < sparsity_pattern->n_rows(), ExcInternalError());

    // both iterators are not one-past-the-end
    return ((current_row < other.current_row) ||
            ((current_row == other.current_row) &&
             (current_entry < other.current_entry)));
  }


  inline void
  Accessor::advance()
  {
    Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());

    // move to the next element in this row
    ++current_entry;

    // if this moves us beyond the end of the row, go to the next row
    // if possible, or set the iterator to an invalid state if not.
    //
    // going to the next row is a bit complicated because we may have
    // to skip over empty rows, and because we also have to avoid rows
    // that aren't listed in a possibly passed IndexSet argument of
    // the sparsity pattern. consequently, rather than trying to
    // duplicate code here, just call the begin() function of the
    // sparsity pattern itself
    if (current_entry == end_of_row)
      {
        if (current_row + 1 < sparsity_pattern->n_rows())
          *this = *sparsity_pattern->begin(current_row + 1);
        else
          *this = Accessor(sparsity_pattern); // invalid object
      }
  }



  inline Iterator::Iterator(const DynamicSparsityPattern *sparsity_pattern,
                            const size_type               row,
                            const unsigned int            index_within_row)
    : accessor(sparsity_pattern, row, index_within_row)
  {}



  inline Iterator::Iterator(const DynamicSparsityPattern *sparsity_pattern)
    : accessor(sparsity_pattern)
  {}



  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
  {
    (void)other;
    Assert(accessor.sparsity_pattern == other.accessor.sparsity_pattern,
           ExcInternalError());
    Assert(false, ExcNotImplemented());

    return 0;
  }
} // namespace DynamicSparsityPatternIterators


inline void
DynamicSparsityPattern::Line::add(const size_type j)
{
  // first check the last element (or if line is still empty)
  if ((entries.size() == 0) || (entries.back() < j))
    {
      entries.push_back(j);
      return;
    }

  // do a binary search to find the place where to insert:
  std::vector<size_type>::iterator it =
    Utilities::lower_bound(entries.begin(), entries.end(), j);

  // If this entry is a duplicate, exit immediately
  if (*it == j)
    return;

  // Insert at the right place in the vector. Vector grows automatically to
  // fit elements. Always doubles its size.
  entries.insert(it, j);
}



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



inline types::global_dof_index
DynamicSparsityPattern::n_cols() const
{
  return cols;
}



inline void
DynamicSparsityPattern::add(const size_type i, const size_type j)
{
  Assert(i < rows, ExcIndexRangeType<size_type>(i, 0, rows));
  Assert(j < cols, ExcIndexRangeType<size_type>(j, 0, cols));

  if (rowset.size() > 0 && !rowset.is_element(i))
    return;

  have_entries = true;

  const size_type rowindex =
    rowset.size() == 0 ? i : rowset.index_within_set(i);
  lines[rowindex].add(j);
}



template <typename ForwardIterator>
inline void
DynamicSparsityPattern::add_entries(const size_type row,
                                    ForwardIterator begin,
                                    ForwardIterator end,
                                    const bool      indices_are_sorted)
{
  Assert(row < rows, ExcIndexRangeType<size_type>(row, 0, rows));

  if (rowset.size() > 0 && !rowset.is_element(row))
    return;

  if (!have_entries && begin < end)
    have_entries = true;

  const size_type rowindex =
    rowset.size() == 0 ? row : rowset.index_within_set(row);
  lines[rowindex].add_entries(begin, end, indices_are_sorted);
}



inline types::global_dof_index
DynamicSparsityPattern::row_length(const size_type row) const
{
  Assert(row < n_rows(), ExcIndexRangeType<size_type>(row, 0, n_rows()));

  if (!have_entries)
    return 0;

  if (rowset.size() > 0 && !rowset.is_element(row))
    return 0;

  const size_type rowindex =
    rowset.size() == 0 ? row : rowset.index_within_set(row);
  return lines[rowindex].entries.size();
}



inline types::global_dof_index
DynamicSparsityPattern::column_number(const size_type row,
                                      const size_type index) const
{
  Assert(row < n_rows(), ExcIndexRangeType<size_type>(row, 0, n_rows()));
  Assert(rowset.size() == 0 || rowset.is_element(row), ExcInternalError());

  const size_type local_row =
    rowset.size() ? rowset.index_within_set(row) : row;
  Assert(index < lines[local_row].entries.size(),
         ExcIndexRangeType<size_type>(index,
                                      0,
                                      lines[local_row].entries.size()));
  return lines[local_row].entries[index];
}



inline DynamicSparsityPattern::iterator
DynamicSparsityPattern::begin() const
{
  return begin(0);
}


inline DynamicSparsityPattern::iterator
DynamicSparsityPattern::end() const
{
  return iterator(this);
}



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

  if (!have_entries)
    return iterator(this);

  if (rowset.size() > 0)
    {
      // We have an IndexSet that describes the locally owned set. For
      // performance reasons we need to make sure that we don't do a
      // linear search over 0..n_rows(). Instead, find the first entry
      // >= row r in the locally owned set (this is done in log
      // n_ranges time inside at()). From there, we move forward until
      // we find a non-empty row. By iterating over the IndexSet instead
      // of incrementing the row index, we potentially skip over entries
      // not in the rowset.
      IndexSet::ElementIterator it = rowset.at(r);
      if (it == rowset.end())
        return end(); // we don't own any row between r and the end

      // Instead of using row_length(*it)==0 in the while loop below,
      // which involves an expensive index_within_set() call, we
      // look at the lines vector directly. This works, because we are
      // walking over this vector entry by entry anyways.
      size_type rowindex = rowset.index_within_set(*it);

      while (it != rowset.end() && lines[rowindex].entries.size() == 0)
        {
          ++it;
          ++rowindex;
        }

      if (it == rowset.end())
        return end();
      else
        return iterator(this, *it, 0);
    }

  // Without an index set we have to do a linear search starting at
  // row r until we find a non-empty one. We will check the lines vector
  // directly instead of going through the slower row_length() function
  size_type row = r;

  while (row < n_rows() && lines[row].entries.size() == 0)
    {
      ++row;
    }

  if (row == n_rows())
    return iterator(this);
  else
    return iterator(this, row, 0);
}



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

  unsigned int row = r + 1;
  if (row == n_rows())
    return iterator(this);
  else
    return begin(row);
}



inline const IndexSet &
DynamicSparsityPattern::row_index_set() const
{
  return rowset;
}



inline bool
DynamicSparsityPattern::stores_only_added_elements()
{
  return true;
}


DEAL_II_NAMESPACE_CLOSE

#endif