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 at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#ifndef dealii_dynamic_sparsity_pattern_h
#define dealii_dynamic_sparsity_pattern_h


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

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

DEAL_II_NAMESPACE_OPEN

class DynamicSparsityPattern;


namespace DynamicSparsityPatternIterators
{
  // forward declaration
  class Iterator;

  typedef types::global_dof_index size_type;

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


class DynamicSparsityPattern : public Subscriptor
{
public:
  typedef types::global_dof_index size_type;

  typedef
  DynamicSparsityPatternIterators::Iterator
  iterator;

  typedef
  DynamicSparsityPatternIterators::Iterator
  const_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;
  }
}


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