chunk_sparsity_pattern.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2008 - 2015 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__chunk_sparsity_pattern_h
#define dealii__chunk_sparsity_pattern_h


#include <deal.II/base/config.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/subscriptor.h>
#include <deal.II/base/vector_slice.h>

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

#include <vector>
#include <iostream>

DEAL_II_NAMESPACE_OPEN


template <typename> class ChunkSparseMatrix;


namespace ChunkSparsityPatternIterators
{
  // forward declaration
  class Iterator;

  class Accessor
  {
  public:
    Accessor (const ChunkSparsityPattern *matrix,
              const unsigned int          row);

    Accessor (const ChunkSparsityPattern *matrix);

    unsigned int row () const;

    std::size_t reduced_index() const;

    unsigned int column () const;

    bool is_valid_entry () const;


    bool operator == (const Accessor &) const;


    bool operator < (const Accessor &) const;

  protected:
    const ChunkSparsityPattern *sparsity_pattern;

    SparsityPatternIterators::Accessor reduced_accessor;

    unsigned int chunk_row;

    unsigned int chunk_col;

    void advance ();

    friend class Iterator;
  };



  class Iterator
  {
  public:
    Iterator (const ChunkSparsityPattern *sp,
              const unsigned int          row);

    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;

  private:
    Accessor accessor;
  };
}



class ChunkSparsityPattern : public Subscriptor
{
public:
  typedef types::global_dof_index size_type;
  typedef ChunkSparsityPatternIterators::Iterator const_iterator;

  typedef ChunkSparsityPatternIterators::Iterator iterator;

  static const size_type invalid_entry = SparsityPattern::invalid_entry;

  ChunkSparsityPattern ();

  ChunkSparsityPattern (const ChunkSparsityPattern &);

  ChunkSparsityPattern (const size_type m,
                        const size_type n,
                        const size_type max_chunks_per_row,
                        const size_type chunk_size);

  ChunkSparsityPattern (const size_type m,
                        const size_type n,
                        const std::vector<size_type> &row_lengths,
                        const size_type chunk_size);

  ChunkSparsityPattern (const size_type n,
                        const size_type max_per_row,
                        const size_type chunk_size);

  ChunkSparsityPattern (const size_type                m,
                        const std::vector<size_type>  &row_lengths,
                        const size_type                chunk_size);

  ~ChunkSparsityPattern ();

  ChunkSparsityPattern &operator = (const ChunkSparsityPattern &);

  void reinit (const size_type m,
               const size_type n,
               const size_type max_per_row,
               const size_type chunk_size);

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

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

  void compress ();

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

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

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

  template <typename Sparsity>
  void create_from (const unsigned int  m,
                    const unsigned int  n,
                    const Sparsity     &sparsity_pattern_for_chunks,
                    const unsigned int  chunk_size,
                    const bool          optimize_diagonal = true);

  bool empty () const;

  size_type get_chunk_size () const;

  size_type max_entries_per_row () const;

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

  void symmetrize ();

  inline size_type n_rows () const;

  inline size_type n_cols () const;

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

  size_type row_length (const size_type row) const;

  size_type bandwidth () const;

  size_type n_nonzero_elements () const;

  bool is_compressed () const;

  bool stores_only_added_elements () const;

  iterator begin () const;

  iterator end () const;

  iterator begin (const unsigned int r) const;

  iterator end (const unsigned int r) 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;

  std::size_t memory_consumption () const;

  DeclException1 (ExcInvalidNumber,
                  int,
                  << "The provided number is invalid here: " << arg1);
  DeclException2 (ExcInvalidIndex,
                  int, int,
                  << "The given index " << arg1
                  << " should be less than " << arg2 << ".");
  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?)");
  DeclException0 (ExcNotCompressed);
  DeclException0 (ExcMatrixIsCompressed);
  DeclException0 (ExcEmptyObject);
  DeclException0 (ExcInvalidConstructorCall);
  DeclException2 (ExcIteratorRange,
                  int, int,
                  << "The iterators denote a range of " << arg1
                  << " elements, but the given number of rows was " << arg2);
  DeclException0 (ExcMETISNotInstalled);
  DeclException1 (ExcInvalidNumberOfPartitions,
                  int,
                  << "The number of partitions you gave is " << arg1
                  << ", but must be greater than zero.");
  DeclException2 (ExcInvalidArraySize,
                  int, int,
                  << "The array has size " << arg1 << " but should have size "
                  << arg2);
private:
  size_type rows;

  size_type cols;

  size_type chunk_size;

  SparsityPattern sparsity_pattern;

  template <typename> friend class ChunkSparseMatrix;

  friend class ChunkSparsityPatternIterators::Accessor;
};


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

#ifndef DOXYGEN

namespace ChunkSparsityPatternIterators
{
  inline
  Accessor::
  Accessor (const ChunkSparsityPattern *sparsity_pattern,
            const unsigned int          row)
    :
    sparsity_pattern(sparsity_pattern),
    reduced_accessor(row==sparsity_pattern->n_rows() ?
                     *sparsity_pattern->sparsity_pattern.end() :
                     *sparsity_pattern->sparsity_pattern.
                     begin(row/sparsity_pattern->get_chunk_size())),
    chunk_row (row==sparsity_pattern->n_rows() ? 0 :
               row%sparsity_pattern->get_chunk_size()),
    chunk_col (0)
  {}



  inline
  Accessor::
  Accessor (const ChunkSparsityPattern *sparsity_pattern)
    :
    sparsity_pattern(sparsity_pattern),
    reduced_accessor(*sparsity_pattern->sparsity_pattern.end()),
    chunk_row (0),
    chunk_col (0)
  {}



  inline
  bool
  Accessor::is_valid_entry () const
  {
    return reduced_accessor.is_valid_entry()
           &&
           sparsity_pattern->get_chunk_size()*reduced_accessor.row()+chunk_row <
           sparsity_pattern->n_rows()
           &&
           sparsity_pattern->get_chunk_size()*reduced_accessor.column()+chunk_col <
           sparsity_pattern->n_cols();
  }



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

    return sparsity_pattern->get_chunk_size()*reduced_accessor.row() +
           chunk_row;
  }



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

    return sparsity_pattern->get_chunk_size()*reduced_accessor.column() +
           chunk_col;
  }



  inline
  std::size_t
  Accessor::reduced_index() const
  {
    Assert (is_valid_entry() == true, ExcInvalidIterator());

    return reduced_accessor.index_within_sparsity;
  }




  inline
  bool
  Accessor::operator == (const Accessor &other) const
  {
    // no need to check for equality of sparsity patterns as this is done in
    // the reduced case already and every ChunkSparsityPattern has its own
    // reduced sparsity pattern
    return (reduced_accessor == other.reduced_accessor &&
            chunk_row == other.chunk_row &&
            chunk_col == other.chunk_col);
  }



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

    if (chunk_row != other.chunk_row)
      {
        if (reduced_accessor.index_within_sparsity ==
            reduced_accessor.sparsity_pattern->n_nonzero_elements())
          return false;
        if (other.reduced_accessor.index_within_sparsity ==
            reduced_accessor.sparsity_pattern->n_nonzero_elements())
          return true;

        const unsigned int
        global_row = sparsity_pattern->get_chunk_size()*reduced_accessor.row()
                     +chunk_row,
                     other_global_row = sparsity_pattern->get_chunk_size()*
                                        other.reduced_accessor.row()+other.chunk_row;
        if (global_row < other_global_row)
          return true;
        else if (global_row > other_global_row)
          return false;
      }

    return (reduced_accessor.index_within_sparsity <
            other.reduced_accessor.index_within_sparsity ||
            (reduced_accessor.index_within_sparsity ==
             other.reduced_accessor.index_within_sparsity &&
             chunk_col < other.chunk_col));
  }


  inline
  void
  Accessor::advance ()
  {
    const unsigned int chunk_size = sparsity_pattern->get_chunk_size();
    Assert (chunk_row < chunk_size && chunk_col < chunk_size,
            ExcIteratorPastEnd());
    Assert (reduced_accessor.row() * chunk_size + chunk_row <
            sparsity_pattern->n_rows()
            &&
            reduced_accessor.column() * chunk_size + chunk_col <
            sparsity_pattern->n_cols(),
            ExcIteratorPastEnd());
    if (chunk_size == 1)
      {
        reduced_accessor.advance();
        return;
      }

    ++chunk_col;

    // end of chunk
    if (chunk_col == chunk_size
        ||
        reduced_accessor.column() * chunk_size + chunk_col ==
        sparsity_pattern->n_cols())
      {
        const unsigned int reduced_row = reduced_accessor.row();
        // end of row
        if (reduced_accessor.index_within_sparsity + 1 ==
            reduced_accessor.sparsity_pattern->rowstart[reduced_row+1])
          {
            ++chunk_row;

            chunk_col = 0;

            // end of chunk rows or end of matrix
            if (chunk_row == chunk_size ||
                (reduced_row * chunk_size + chunk_row ==
                 sparsity_pattern->n_rows()))
              {
                chunk_row = 0;
                reduced_accessor.advance();
              }
            // go back to the beginning of the same reduced row but with
            // chunk_row increased by one
            else
              reduced_accessor.index_within_sparsity =
                reduced_accessor.sparsity_pattern->rowstart[reduced_row];
          }
        // advance within chunk
        else
          {
            reduced_accessor.advance();
            chunk_col = 0;
          }
      }
  }



  inline
  Iterator::Iterator (const ChunkSparsityPattern *sparsity_pattern,
                      const unsigned int          row)
    :
    accessor(sparsity_pattern, row)
  {}



  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 ! (accessor == other.accessor);
  }


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

}



inline
ChunkSparsityPattern::iterator
ChunkSparsityPattern::begin () const
{
  return iterator(this, 0);
}


inline
ChunkSparsityPattern::iterator
ChunkSparsityPattern::end () const
{
  return iterator(this, n_rows());
}



inline
ChunkSparsityPattern::iterator
ChunkSparsityPattern::begin (const unsigned int r) const
{
  Assert (r<n_rows(), ExcIndexRange(r,0,n_rows()));
  return iterator(this, r);
}



inline
ChunkSparsityPattern::iterator
ChunkSparsityPattern::end (const unsigned int r) const
{
  Assert (r<n_rows(), ExcIndexRange(r,0,n_rows()))
  return iterator(this, r+1);
}



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


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



inline
ChunkSparsityPattern::size_type
ChunkSparsityPattern::get_chunk_size () const
{
  return chunk_size;
}



inline
bool
ChunkSparsityPattern::is_compressed () const
{
  return sparsity_pattern.compressed;
}



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

  // now enter all the elements into the matrix
  size_type row = 0;
  typedef typename std::iterator_traits<ForwardIterator>::value_type::const_iterator inner_iterator;
  for (ForwardIterator i=begin; i!=end; ++i, ++row)
    {
      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, ExcInvalidIndex(col,n_cols));

          add (row, col);
        }
    }

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


#endif // DOXYGEN

DEAL_II_NAMESPACE_CLOSE

#endif