chunk_sparse_matrix.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2008 - 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_chunk_sparse_matrix_h
#define dealii_chunk_sparse_matrix_h


#include <deal.II/base/config.h>
#include <deal.II/base/subscriptor.h>
#include <deal.II/base/smartpointer.h>
#include <deal.II/lac/chunk_sparsity_pattern.h>
#include <deal.II/lac/identity_matrix.h>
#include <deal.II/lac/exceptions.h>

#include <memory>


DEAL_II_NAMESPACE_OPEN

template <typename number> class Vector;
template <typename number> class FullMatrix;

namespace ChunkSparseMatrixIterators
{
  // forward declaration
  template <typename number, bool Constness>
  class Iterator;

  template <typename number, bool Constness>
  class Accessor : public ChunkSparsityPatternIterators::Accessor
  {
  public:
    number value() const;

    number &value();

    const ChunkSparseMatrix<number> &get_matrix () const;
  };



  template <typename number>
  class Accessor<number,true> : public ChunkSparsityPatternIterators::Accessor
  {
  public:
    typedef const ChunkSparseMatrix<number> MatrixType;

    Accessor (MatrixType         *matrix,
              const unsigned int  row);

    Accessor (MatrixType         *matrix);

    Accessor (const ChunkSparseMatrixIterators::Accessor<number,false> &a);

    number value() const;

    const MatrixType &get_matrix () const;

  private:
    MatrixType *matrix;

    using ChunkSparsityPatternIterators::Accessor::advance;

    template <typename, bool>
    friend class Iterator;
  };


  template <typename number>
  class Accessor<number,false> : public ChunkSparsityPatternIterators::Accessor
  {
  private:
    class Reference
    {
    public:
      Reference (const Accessor *accessor,
                 const bool dummy);

      operator number () const;

      const Reference &operator = (const number n) const;

      const Reference &operator += (const number n) const;

      const Reference &operator -= (const number n) const;

      const Reference &operator *= (const number n) const;

      const Reference &operator /= (const number n) const;

    private:
      const Accessor *accessor;
    };

  public:
    typedef ChunkSparseMatrix<number> MatrixType;

    Accessor (MatrixType         *matrix,
              const unsigned int  row);

    Accessor (MatrixType         *matrix);

    Reference value() const;

    MatrixType &get_matrix () const;

  private:
    MatrixType *matrix;

    using ChunkSparsityPatternIterators::Accessor::advance;

    template <typename, bool>
    friend class Iterator;
  };



  template <typename number, bool Constness>
  class Iterator
  {
  public:
    typedef
    typename Accessor<number,Constness>::MatrixType
    MatrixType;

    typedef
    const Accessor<number,Constness> &value_type;

    Iterator (MatrixType        *matrix,
              const unsigned int row);

    Iterator (MatrixType *matrix);

    Iterator (const ChunkSparseMatrixIterators::Iterator<number,false> &i);

    Iterator &operator++ ();

    Iterator operator++ (int);

    const Accessor<number,Constness> &operator* () const;

    const Accessor<number,Constness> *operator-> () const;

    bool operator == (const Iterator &) const;

    bool operator != (const Iterator &) const;

    bool operator < (const Iterator &) const;

    bool operator > (const Iterator &) const;

    int operator - (const Iterator &p) const;

    Iterator operator + (const unsigned int n) const;

  private:
    Accessor<number,Constness> accessor;
  };

}



template <typename number>
class ChunkSparseMatrix : public virtual Subscriptor
{
public:
  typedef types::global_dof_index size_type;

  typedef number value_type;

  typedef typename numbers::NumberTraits<number>::real_type real_type;

  typedef
  ChunkSparseMatrixIterators::Iterator<number,true>
  const_iterator;

  typedef
  ChunkSparseMatrixIterators::Iterator<number,false>
  iterator;

  struct Traits
  {
    static const bool zero_addition_can_be_elided = true;
  };

  ChunkSparseMatrix ();

  ChunkSparseMatrix (const ChunkSparseMatrix &);

  explicit ChunkSparseMatrix (const ChunkSparsityPattern &sparsity);

  ChunkSparseMatrix (const ChunkSparsityPattern &sparsity,
                     const IdentityMatrix  &id);

  virtual ~ChunkSparseMatrix ();

  ChunkSparseMatrix<number> &operator = (const ChunkSparseMatrix<number> &);

  ChunkSparseMatrix<number> &
  operator= (const IdentityMatrix  &id);

  ChunkSparseMatrix &operator = (const double d);

  virtual void reinit (const ChunkSparsityPattern &sparsity);

  virtual void clear ();

  bool empty () const;

  size_type m () const;

  size_type n () const;

  size_type n_nonzero_elements () const;

  size_type n_actually_nonzero_elements () const;

  const ChunkSparsityPattern &get_sparsity_pattern () const;

  std::size_t memory_consumption () const;


  void set (const size_type i,
            const size_type j,
            const number value);

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

  template <typename number2>
  void add (const size_type  row,
            const size_type  n_cols,
            const size_type *col_indices,
            const number2   *values,
            const bool       elide_zero_values = true,
            const bool       col_indices_are_sorted = false);

  ChunkSparseMatrix &operator *= (const number factor);

  ChunkSparseMatrix &operator /= (const number factor);

  void symmetrize ();

  template <typename somenumber>
  ChunkSparseMatrix<number> &
  copy_from (const ChunkSparseMatrix<somenumber> &source);

  template <typename ForwardIterator>
  void copy_from (const ForwardIterator begin,
                  const ForwardIterator end);

  template <typename somenumber>
  void copy_from (const FullMatrix<somenumber> &matrix);

  template <typename somenumber>
  void add (const number factor,
            const ChunkSparseMatrix<somenumber> &matrix);



  number operator () (const size_type i,
                      const size_type j) const;

  number el (const size_type i,
             const size_type j) const;

  number diag_element (const size_type i) const;

  void extract_row_copy (const size_type row,
                         const size_type array_length,
                         size_type      &row_length,
                         size_type      *column_indices,
                         number         *values) const;


  template <class OutVector, class InVector>
  void vmult (OutVector &dst,
              const InVector &src) const;

  template <class OutVector, class InVector>
  void Tvmult (OutVector &dst,
               const InVector &src) const;

  template <class OutVector, class InVector>
  void vmult_add (OutVector &dst,
                  const InVector &src) const;

  template <class OutVector, class InVector>
  void Tvmult_add (OutVector &dst,
                   const InVector &src) const;

  template <typename somenumber>
  somenumber matrix_norm_square (const Vector<somenumber> &v) const;

  template <typename somenumber>
  somenumber matrix_scalar_product (const Vector<somenumber> &u,
                                    const Vector<somenumber> &v) const;
  template <typename somenumber>
  somenumber residual (Vector<somenumber>       &dst,
                       const Vector<somenumber> &x,
                       const Vector<somenumber> &b) const;



  real_type l1_norm () const;

  real_type linfty_norm () const;

  real_type frobenius_norm () const;


  template <typename somenumber>
  void precondition_Jacobi (Vector<somenumber>       &dst,
                            const Vector<somenumber> &src,
                            const number              omega = 1.) const;

  template <typename somenumber>
  void precondition_SSOR (Vector<somenumber>       &dst,
                          const Vector<somenumber> &src,
                          const number              om = 1.) const;

  template <typename somenumber>
  void precondition_SOR (Vector<somenumber>       &dst,
                         const Vector<somenumber> &src,
                         const number              om = 1.) const;

  template <typename somenumber>
  void precondition_TSOR (Vector<somenumber>       &dst,
                          const Vector<somenumber> &src,
                          const number              om = 1.) const;

  template <typename somenumber>
  void SSOR (Vector<somenumber> &v,
             const number        omega = 1.) const;

  template <typename somenumber>
  void SOR (Vector<somenumber> &v,
            const number        om = 1.) const;

  template <typename somenumber>
  void TSOR (Vector<somenumber> &v,
             const number        om = 1.) const;

  template <typename somenumber>
  void PSOR (Vector<somenumber> &v,
             const std::vector<size_type> &permutation,
             const std::vector<size_type> &inverse_permutation,
             const number        om = 1.) const;

  template <typename somenumber>
  void TPSOR (Vector<somenumber> &v,
              const std::vector<size_type> &permutation,
              const std::vector<size_type> &inverse_permutation,
              const number        om = 1.) const;

  template <typename somenumber>
  void SOR_step (Vector<somenumber> &v,
                 const Vector<somenumber> &b,
                 const number        om = 1.) const;

  template <typename somenumber>
  void TSOR_step (Vector<somenumber> &v,
                  const Vector<somenumber> &b,
                  const number        om = 1.) const;

  template <typename somenumber>
  void SSOR_step (Vector<somenumber> &v,
                  const Vector<somenumber> &b,
                  const number        om = 1.) const;


  const_iterator begin () const;

  const_iterator end () const;

  iterator begin ();

  iterator end ();

  const_iterator begin (const unsigned int r) const;

  const_iterator end (const unsigned int r) const;

  iterator begin (const unsigned int r);

  iterator end (const unsigned int r);


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

  void print_formatted (std::ostream       &out,
                        const unsigned int  precision   = 3,
                        const bool          scientific  = true,
                        const unsigned int  width       = 0,
                        const char         *zero_string = " ",
                        const double        denominator = 1.) const;

  void print_pattern(std::ostream &out,
                     const double threshold = 0.) const;

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

  void block_read (std::istream &in);

  DeclException2 (ExcInvalidIndex,
                  int, int,
                  << "You are trying to access the matrix entry with index <"
                  << arg1 << ',' << arg2
                  << ">, but this entry does not exist in the sparsity pattern"
                  "of this matrix."
                  "\n\n"
                  "The most common cause for this problem is that you used "
                  "a method to build the sparsity pattern that did not "
                  "(completely) take into account all of the entries you "
                  "will later try to write into. An example would be "
                  "building a sparsity pattern that does not include "
                  "the entries you will write into due to constraints "
                  "on degrees of freedom such as hanging nodes or periodic "
                  "boundary conditions. In such cases, building the "
                  "sparsity pattern will succeed, but you will get errors "
                  "such as the current one at one point or other when "
                  "trying to write into the entries of the matrix.");
  DeclException0 (ExcDifferentChunkSparsityPatterns);
  DeclException2 (ExcIteratorRange,
                  int, int,
                  << "The iterators denote a range of " << arg1
                  << " elements, but the given number of rows was " << arg2);
  DeclExceptionMsg (ExcSourceEqualsDestination,
                    "You are attempting an operation on two matrices that "
                    "are the same object, but the operation requires that the "
                    "two objects are in fact different.");
private:
  SmartPointer<const ChunkSparsityPattern,ChunkSparseMatrix<number> > cols;

  std::unique_ptr<number[]> val;

  size_type max_len;

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

  // make all other sparse matrices friends
  template <typename somenumber> friend class ChunkSparseMatrix;

  template <typename,bool> friend class ChunkSparseMatrixIterators::Iterator;
  template <typename,bool> friend class ChunkSparseMatrixIterators::Accessor;
};

#ifndef DOXYGEN
/*---------------------- Inline functions -----------------------------------*/



template <typename number>
inline
typename ChunkSparseMatrix<number>::size_type
ChunkSparseMatrix<number>::m () const
{
  Assert (cols != nullptr, ExcNotInitialized());
  return cols->rows;
}


template <typename number>
inline
typename ChunkSparseMatrix<number>::size_type
ChunkSparseMatrix<number>::n () const
{
  Assert (cols != nullptr, ExcNotInitialized());
  return cols->cols;
}



template <typename number>
inline
const ChunkSparsityPattern &
ChunkSparseMatrix<number>::get_sparsity_pattern () const
{
  Assert (cols != nullptr, ExcNotInitialized());
  return *cols;
}



template <typename number>
inline
typename ChunkSparseMatrix<number>::size_type
ChunkSparseMatrix<number>::compute_location (const size_type i,
                                             const size_type j) const
{
  const size_type chunk_size = cols->get_chunk_size();
  const size_type chunk_index
    = cols->sparsity_pattern(i/chunk_size, j/chunk_size);

  if (chunk_index == ChunkSparsityPattern::invalid_entry)
    return ChunkSparsityPattern::invalid_entry;
  else
    {
      return (chunk_index * chunk_size * chunk_size
              +
              (i % chunk_size) * chunk_size
              +
              (j % chunk_size));
    }
}


template <typename number>
inline
void ChunkSparseMatrix<number>::set (const size_type i,
                                     const size_type j,
                                     const number value)
{

  AssertIsFinite(value);

  Assert (cols != nullptr, ExcNotInitialized());
  // it is allowed to set elements of the matrix that are not part of the
  // sparsity pattern, if the value to which we set it is zero
  const size_type index = compute_location(i,j);
  Assert ((index != SparsityPattern::invalid_entry) ||
          (value == 0.),
          ExcInvalidIndex(i,j));

  if (index != SparsityPattern::invalid_entry)
    val[index] = value;
}



template <typename number>
inline
void ChunkSparseMatrix<number>::add (const size_type i,
                                     const size_type j,
                                     const number value)
{

  AssertIsFinite(value);

  Assert (cols != nullptr, ExcNotInitialized());

  if (value != 0.)
    {
      const size_type index = compute_location(i,j);
      Assert ((index != ChunkSparsityPattern::invalid_entry),
              ExcInvalidIndex(i,j));

      val[index] += value;
    }
}



template <typename number>
template <typename number2>
inline
void ChunkSparseMatrix<number>::add (const size_type  row,
                                     const size_type  n_cols,
                                     const size_type *col_indices,
                                     const number2   *values,
                                     const bool       /*elide_zero_values*/,
                                     const bool       /*col_indices_are_sorted*/)
{
  // TODO: could be done more efficiently...
  for (size_type col=0; col<n_cols; ++col)
    add(row, col_indices[col], static_cast<number>(values[col]));
}



template <typename number>
inline
ChunkSparseMatrix<number> &
ChunkSparseMatrix<number>::operator *= (const number factor)
{
  Assert (cols != nullptr, ExcNotInitialized());
  Assert (val != nullptr, ExcNotInitialized());

  const size_type chunk_size = cols->get_chunk_size();

  // multiply all elements of the matrix with the given factor. this includes
  // the padding elements in chunks that overlap the boundaries of the actual
  // matrix -- but since multiplication with a number does not violate the
  // invariant of keeping these elements at zero nothing can happen
  number             *val_ptr    = val.get();
  const number *const end_ptr    = val.get() +
                                   cols->sparsity_pattern.n_nonzero_elements()
                                   *
                                   chunk_size * chunk_size;
  while (val_ptr != end_ptr)
    *val_ptr++ *= factor;

  return *this;
}



template <typename number>
inline
ChunkSparseMatrix<number> &
ChunkSparseMatrix<number>::operator /= (const number factor)
{
  Assert (cols != nullptr, ExcNotInitialized());
  Assert (val != nullptr, ExcNotInitialized());
  Assert (factor !=0, ExcDivideByZero());

  const number factor_inv = 1. / factor;

  const size_type chunk_size = cols->get_chunk_size();

  // multiply all elements of the matrix with the given factor. this includes
  // the padding elements in chunks that overlap the boundaries of the actual
  // matrix -- but since multiplication with a number does not violate the
  // invariant of keeping these elements at zero nothing can happen
  number             *val_ptr    = val.get();
  const number *const end_ptr    = val.get() +
                                   cols->sparsity_pattern.n_nonzero_elements()
                                   *
                                   chunk_size * chunk_size;

  while (val_ptr != end_ptr)
    *val_ptr++ *= factor_inv;

  return *this;
}



template <typename number>
inline
number ChunkSparseMatrix<number>::operator () (const size_type i,
                                               const size_type j) const
{
  Assert (cols != nullptr, ExcNotInitialized());
  AssertThrow (compute_location(i,j) != SparsityPattern::invalid_entry,
               ExcInvalidIndex(i,j));
  return val[compute_location(i,j)];
}



template <typename number>
inline
number ChunkSparseMatrix<number>::el (const size_type i,
                                      const size_type j) const
{
  Assert (cols != nullptr, ExcNotInitialized());
  const size_type index = compute_location(i,j);

  if (index != ChunkSparsityPattern::invalid_entry)
    return val[index];
  else
    return 0;
}



template <typename number>
inline
number ChunkSparseMatrix<number>::diag_element (const size_type i) const
{
  Assert (cols != nullptr, ExcNotInitialized());
  Assert (m() == n(),  ExcNotQuadratic());
  AssertIndexRange(i, m());

  // Use that the first element in each row of a quadratic matrix is the main
  // diagonal of the chunk sparsity pattern
  const size_type chunk_size = cols->get_chunk_size();
  return val[cols->sparsity_pattern.rowstart[i/chunk_size]
             *
             chunk_size * chunk_size
             +
             (i % chunk_size) * chunk_size
             +
             (i % chunk_size)];
}



template <typename number>
template <typename ForwardIterator>
inline
void
ChunkSparseMatrix<number>::copy_from (const ForwardIterator begin,
                                      const ForwardIterator end)
{
  Assert (static_cast<size_type >(std::distance (begin, end)) == m(),
          ExcIteratorRange (std::distance (begin, end), m()));

  // for use in the inner loop, we define a typedef to the type of the inner
  // iterators
  typedef typename std::iterator_traits<ForwardIterator>::value_type::const_iterator inner_iterator;
  size_type row=0;
  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)
        // write entries
        set (row, j->first, j->second);
    }
}



//---------------------------------------------------------------------------


namespace ChunkSparseMatrixIterators
{
  template <typename number>
  inline
  Accessor<number,true>::
  Accessor (const MatrixType   *matrix,
            const unsigned int  row)
    :
    ChunkSparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern(),
                                             row),
    matrix (matrix)
  {}



  template <typename number>
  inline
  Accessor<number,true>::
  Accessor (const MatrixType *matrix)
    :
    ChunkSparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern()),
    matrix (matrix)
  {}



  template <typename number>
  inline
  Accessor<number,true>::
  Accessor (const ChunkSparseMatrixIterators::Accessor<number,false> &a)
    :
    ChunkSparsityPatternIterators::Accessor (a),
    matrix (&a.get_matrix())
  {}



  template <typename number>
  inline
  number
  Accessor<number, true>::value () const
  {
    const unsigned int chunk_size = matrix->get_sparsity_pattern().get_chunk_size();
    return matrix->val[reduced_index() * chunk_size * chunk_size
                       +
                       chunk_row * chunk_size
                       +
                       chunk_col];
  }



  template <typename number>
  inline
  const typename Accessor<number, true>::MatrixType &
  Accessor<number, true>::get_matrix () const
  {
    return *matrix;
  }



  template <typename number>
  inline
  Accessor<number, false>::Reference::Reference (
    const Accessor *accessor,
    const bool)
    :
    accessor (accessor)
  {}


  template <typename number>
  inline
  Accessor<number, false>::Reference::operator number() const
  {
    const unsigned int chunk_size = accessor->matrix->get_sparsity_pattern().get_chunk_size();
    return accessor->matrix->val[accessor->reduced_index() * chunk_size * chunk_size
                                 +
                                 accessor->chunk_row * chunk_size
                                 +
                                 accessor->chunk_col];
  }



  template <typename number>
  inline
  const typename Accessor<number, false>::Reference &
  Accessor<number, false>::Reference::operator = (const number n) const
  {
    const unsigned int chunk_size = accessor->matrix->get_sparsity_pattern().get_chunk_size();
    accessor->matrix->val[accessor->reduced_index() * chunk_size * chunk_size
                          +
                          accessor->chunk_row * chunk_size
                          +
                          accessor->chunk_col] = n;
    return *this;
  }



  template <typename number>
  inline
  const typename Accessor<number, false>::Reference &
  Accessor<number, false>::Reference::operator += (const number n) const
  {
    const unsigned int chunk_size = accessor->matrix->get_sparsity_pattern().get_chunk_size();
    accessor->matrix->val[accessor->reduced_index() * chunk_size * chunk_size
                          +
                          accessor->chunk_row * chunk_size
                          +
                          accessor->chunk_col] += n;
    return *this;
  }



  template <typename number>
  inline
  const typename Accessor<number, false>::Reference &
  Accessor<number, false>::Reference::operator -= (const number n) const
  {
    const unsigned int chunk_size = accessor->matrix->get_sparsity_pattern().get_chunk_size();
    accessor->matrix->val[accessor->reduced_index() * chunk_size * chunk_size
                          +
                          accessor->chunk_row * chunk_size
                          +
                          accessor->chunk_col] -= n;
    return *this;
  }



  template <typename number>
  inline
  const typename Accessor<number, false>::Reference &
  Accessor<number, false>::Reference::operator *= (const number n) const
  {
    const unsigned int chunk_size = accessor->matrix->get_sparsity_pattern().get_chunk_size();
    accessor->matrix->val[accessor->reduced_index() * chunk_size * chunk_size
                          +
                          accessor->chunk_row * chunk_size
                          +
                          accessor->chunk_col] *= n;
    return *this;
  }



  template <typename number>
  inline
  const typename Accessor<number, false>::Reference &
  Accessor<number, false>::Reference::operator /= (const number n) const
  {
    const unsigned int chunk_size = accessor->matrix->get_sparsity_pattern().get_chunk_size();
    accessor->matrix->val[accessor->reduced_index() * chunk_size * chunk_size
                          +
                          accessor->chunk_row * chunk_size
                          +
                          accessor->chunk_col] /= n;
    return *this;
  }



  template <typename number>
  inline
  Accessor<number,false>::
  Accessor (MatrixType         *matrix,
            const unsigned int  row)
    :
    ChunkSparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern(),
                                             row),
    matrix (matrix)
  {}



  template <typename number>
  inline
  Accessor<number,false>::
  Accessor (MatrixType         *matrix)
    :
    ChunkSparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern()),
    matrix (matrix)
  {}



  template <typename number>
  inline
  typename Accessor<number, false>::Reference
  Accessor<number, false>::value() const
  {
    return Reference(this,true);
  }




  template <typename number>
  inline
  typename Accessor<number, false>::MatrixType &
  Accessor<number, false>::get_matrix () const
  {
    return *matrix;
  }



  template <typename number, bool Constness>
  inline
  Iterator<number, Constness>::
  Iterator (MatrixType        *matrix,
            const unsigned int row)
    :
    accessor(matrix, row)
  {}



  template <typename number, bool Constness>
  inline
  Iterator<number, Constness>::
  Iterator (MatrixType *matrix)
    :
    accessor(matrix)
  {}



  template <typename number, bool Constness>
  inline
  Iterator<number, Constness>::
  Iterator (const ChunkSparseMatrixIterators::Iterator<number,false> &i)
    :
    accessor(*i)
  {}



  template <typename number, bool Constness>
  inline
  Iterator<number, Constness> &
  Iterator<number,Constness>::operator++ ()
  {
    accessor.advance ();
    return *this;
  }


  template <typename number, bool Constness>
  inline
  Iterator<number,Constness>
  Iterator<number,Constness>::operator++ (int)
  {
    const Iterator iter = *this;
    accessor.advance ();
    return iter;
  }


  template <typename number, bool Constness>
  inline
  const Accessor<number,Constness> &
  Iterator<number,Constness>::operator* () const
  {
    return accessor;
  }


  template <typename number, bool Constness>
  inline
  const Accessor<number,Constness> *
  Iterator<number,Constness>::operator-> () const
  {
    return &accessor;
  }


  template <typename number, bool Constness>
  inline
  bool
  Iterator<number,Constness>::
  operator == (const Iterator &other) const
  {
    return (accessor == other.accessor);
  }


  template <typename number, bool Constness>
  inline
  bool
  Iterator<number,Constness>::
  operator != (const Iterator &other) const
  {
    return ! (*this == other);
  }


  template <typename number, bool Constness>
  inline
  bool
  Iterator<number,Constness>::
  operator < (const Iterator &other) const
  {
    Assert (&accessor.get_matrix() == &other.accessor.get_matrix(),
            ExcInternalError());

    return (accessor < other.accessor);
  }


  template <typename number, bool Constness>
  inline
  bool
  Iterator<number,Constness>::
  operator > (const Iterator &other) const
  {
    return (other < *this);
  }


  template <typename number, bool Constness>
  inline
  int
  Iterator<number,Constness>::
  operator - (const Iterator &other) const
  {
    Assert (&accessor.get_matrix() == &other.accessor.get_matrix(),
            ExcInternalError());

    // TODO: can be optimized
    int difference = 0;
    if (*this < other)
      {
        Iterator copy = *this;
        while (copy != other)
          {
            ++copy;
            --difference;
          }
      }
    else
      {
        Iterator copy = other;
        while (copy != *this)
          {
            ++copy;
            ++difference;
          }
      }
    return difference;
  }



  template <typename number, bool Constness>
  inline
  Iterator<number,Constness>
  Iterator<number,Constness>::
  operator + (const unsigned int n) const
  {
    Iterator x = *this;
    for (unsigned int i=0; i<n; ++i)
      ++x;

    return x;
  }

}



template <typename number>
inline
typename ChunkSparseMatrix<number>::const_iterator
ChunkSparseMatrix<number>::begin () const
{
  return const_iterator(this, 0);
}


template <typename number>
inline
typename ChunkSparseMatrix<number>::const_iterator
ChunkSparseMatrix<number>::end () const
{
  return const_iterator(this);
}


template <typename number>
inline
typename ChunkSparseMatrix<number>::iterator
ChunkSparseMatrix<number>::begin ()
{
  return iterator(this, 0);
}


template <typename number>
inline
typename ChunkSparseMatrix<number>::iterator
ChunkSparseMatrix<number>::end ()
{
  return iterator(this);
}


template <typename number>
inline
typename ChunkSparseMatrix<number>::const_iterator
ChunkSparseMatrix<number>::begin (const unsigned int r) const
{
  Assert (r<m(), ExcIndexRange(r,0,m()));
  return const_iterator(this, r);
}



template <typename number>
inline
typename ChunkSparseMatrix<number>::const_iterator
ChunkSparseMatrix<number>::end (const unsigned int r) const
{
  Assert (r<m(), ExcIndexRange(r,0,m()));
  return const_iterator(this, r+1);
}



template <typename number>
inline
typename ChunkSparseMatrix<number>::iterator
ChunkSparseMatrix<number>::begin (const unsigned int r)
{
  Assert (r<m(), ExcIndexRange(r,0,m()));
  return iterator(this, r);
}



template <typename number>
inline
typename ChunkSparseMatrix<number>::iterator
ChunkSparseMatrix<number>::end (const unsigned int r)
{
  Assert (r<m(), ExcIndexRange(r,0,m()));
  return iterator(this, r+1);
}




#endif // DOXYGEN


/*----------------------------   chunk_sparse_matrix.h     ---------------------------*/

DEAL_II_NAMESPACE_CLOSE

#endif
/*----------------------------   chunk_sparse_matrix.h     ---------------------------*/