sparsity_pattern.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2000 - 2016 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__sparsity_pattern_h
#define dealii__sparsity_pattern_h


#include <deal.II/base/config.h>
#include <deal.II/base/exceptions.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 <vector>
#include <iostream>
#include <algorithm>

DEAL_II_NAMESPACE_OPEN

class SparsityPattern;
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
  {
    typedef types::global_dof_index size_type;

    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 SparsityPatternIterators
{
  // forward declaration
  class Iterator;

  typedef types::global_dof_index size_type;

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



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

  typedef
  SparsityPatternIterators::Iterator
  const_iterator;

  typedef
  SparsityPatternIterators::Iterator
  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 ();

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

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


  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::size_t *rowstart;

  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,
                       sparsity_pattern->rowstart + sparsity_pattern->rows + 1,
                       index_within_sparsity);
    return insert_point - sparsity_pattern->rowstart - 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
  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;
  }
}



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 ((rowstart!=0) && (colnums!=0), ExcEmptyObject());
  Assert (compressed, ExcNotCompressed());
  return rowstart[rows]-rowstart[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, max_dim + 1);
  ar &boost::serialization::make_array(colnums, 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;

  if (rowstart != 0)
    delete[] rowstart;
  rowstart = new std::size_t[max_dim + 1];

  if (colnums != 0)
    delete[] colnums;
  colnums = new size_type[max_vec_len];

  ar &boost::serialization::make_array(rowstart, max_dim + 1);
  ar &boost::serialization::make_array(colnums, 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
  {
    typedef types::global_dof_index size_type;

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



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 a typedef to the type of the inner
  // iterators
  size_type row = 0;
  typedef typename std::iterator_traits<ForwardIterator>::value_type::const_iterator inner_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