chunk_sparse_matrix.h

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


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

#  include <deal.II/base/smartpointer.h>
#  include <deal.II/base/subscriptor.h>

#  include <deal.II/lac/chunk_sparsity_pattern.h>
#  include <deal.II/lac/exceptions.h>
#  include <deal.II/lac/identity_matrix.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:
    using MatrixType = const ChunkSparseMatrix<number>;

    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:
    using MatrixType = ChunkSparseMatrix<number>;

    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:
    using MatrixType = typename Accessor<number, Constness>::MatrixType;

    using value_type = const Accessor<number, Constness> &;

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

} // namespace ChunkSparseMatrixIterators



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

  using value_type = number;

  using real_type = typename numbers::NumberTraits<number>::real_type;

  using const_iterator = ChunkSparseMatrixIterators::Iterator<number, true>;

  using iterator = ChunkSparseMatrixIterators::Iterator<number, false>;

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

  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 (std::abs(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(std::abs(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 an alias to the type of the inner
  // iterators
  using inner_iterator =
    typename std::iterator_traits<ForwardIterator>::value_type::const_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;
  }

} // namespace ChunkSparseMatrixIterators



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

DEAL_II_NAMESPACE_CLOSE

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