sparse_matrix.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1999 - 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 at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#ifndef dealii_sparse_matrix_h
#define dealii_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/sparsity_pattern.h>
#include <deal.II/lac/identity_matrix.h>
#include <deal.II/lac/exceptions.h>
#include <deal.II/lac/vector_operation.h>
#ifdef DEAL_II_WITH_MPI
#include <mpi.h>
#endif

#include <memory>


DEAL_II_NAMESPACE_OPEN

template <typename number> class Vector;
template <typename number> class FullMatrix;
template <typename Matrix> class BlockMatrixBase;
template <typename number> class SparseILU;
#ifdef DEAL_II_WITH_MPI
namespace Utilities
{
  namespace MPI
  {
    template <typename Number> void sum (const SparseMatrix<Number> &, const MPI_Comm &, SparseMatrix<Number> &);
  }
}
#endif

#ifdef DEAL_II_WITH_TRILINOS
namespace TrilinosWrappers
{
  class SparseMatrix;
}
#endif

namespace SparseMatrixIterators
{
  typedef types::global_dof_index size_type;

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

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

    number &value();

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



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

    Accessor (MatrixType         *matrix,
              const std::size_t   index_within_matrix);

    Accessor (MatrixType         *matrix);

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

    number value() const;

    const MatrixType &get_matrix () const;

  private:
    MatrixType *matrix;

    using SparsityPatternIterators::Accessor::advance;

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


  template <typename number>
  class Accessor<number,false> : public SparsityPatternIterators::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 SparseMatrix<number> MatrixType;

    Accessor (MatrixType         *matrix,
              const std::size_t   index);

    Accessor (MatrixType         *matrix);

    Reference value() const;

    MatrixType &get_matrix () const;

  private:
    MatrixType *matrix;

    using SparsityPatternIterators::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 std::size_t  index_within_matrix);

    Iterator (MatrixType *matrix);

    Iterator (const SparseMatrixIterators::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 size_type n) const;

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

}

//TODO: Add multithreading to the other vmult functions.

template <typename number>
class SparseMatrix : 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
  SparseMatrixIterators::Iterator<number,true>
  const_iterator;

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

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

  SparseMatrix ();

  SparseMatrix (const SparseMatrix &);

  SparseMatrix (SparseMatrix<number> &&m) noexcept;

  explicit SparseMatrix (const SparsityPattern &sparsity);

  SparseMatrix (const SparsityPattern &sparsity,
                const IdentityMatrix  &id);

  virtual ~SparseMatrix ();

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

  SparseMatrix<number> &operator = (SparseMatrix<number> &&m) noexcept;

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

  SparseMatrix &operator = (const double d);

  virtual void reinit (const SparsityPattern &sparsity);

  virtual void clear ();

  bool empty () const;

  size_type m () const;

  size_type n () const;

  size_type get_row_length (const size_type row) const;

  std::size_t n_nonzero_elements () const;

  std::size_t n_actually_nonzero_elements (const double threshold = 0.) const;

  const SparsityPattern &get_sparsity_pattern () const;

  std::size_t memory_consumption () const;

  void compress (::VectorOperation::values);


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

  template <typename number2>
  void set (const std::vector<size_type> &indices,
            const FullMatrix<number2>    &full_matrix,
            const bool                    elide_zero_values = false);

  template <typename number2>
  void set (const std::vector<size_type> &row_indices,
            const std::vector<size_type> &col_indices,
            const FullMatrix<number2>    &full_matrix,
            const bool                    elide_zero_values = false);

  template <typename number2>
  void set (const size_type               row,
            const std::vector<size_type> &col_indices,
            const std::vector<number2>   &values,
            const bool                    elide_zero_values = false);

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

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

  template <typename number2>
  void add (const std::vector<size_type> &indices,
            const FullMatrix<number2>    &full_matrix,
            const bool                    elide_zero_values = true);

  template <typename number2>
  void add (const std::vector<size_type> &row_indices,
            const std::vector<size_type> &col_indices,
            const FullMatrix<number2>    &full_matrix,
            const bool                    elide_zero_values = true);

  template <typename number2>
  void add (const size_type               row,
            const std::vector<size_type> &col_indices,
            const std::vector<number2>   &values,
            const bool                    elide_zero_values = true);

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

  SparseMatrix &operator *= (const number factor);

  SparseMatrix &operator /= (const number factor);

  void symmetrize ();

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

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

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

#ifdef DEAL_II_WITH_TRILINOS

  SparseMatrix<number> &
  copy_from (const TrilinosWrappers::SparseMatrix &matrix);
#endif

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



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

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

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

  number diag_element (const size_type i) const;

  number &diag_element (const size_type i);


  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;

  template <typename numberB, typename numberC>
  void mmult (SparseMatrix<numberC>       &C,
              const SparseMatrix<numberB> &B,
              const Vector<number>        &V = Vector<number>(),
              const bool                   rebuild_sparsity_pattern = true) const;

  template <typename numberB, typename numberC>
  void Tmmult (SparseMatrix<numberC>       &C,
               const SparseMatrix<numberB> &B,
               const Vector<number>       &V = Vector<number>(),
               const bool                   rebuild_sparsity_pattern = true) 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                    omega = 1.,
                          const std::vector<std::size_t> &pos_right_of_diagonal=std::vector<std::size_t>()) 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 Jacobi_step (Vector<somenumber> &v,
                    const Vector<somenumber> &b,
                    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;

  iterator begin ();

  const_iterator end () const;

  iterator end ();

  const_iterator begin (const size_type r) const;

  iterator begin (const size_type r);

  const_iterator end (const size_type r) const;

  iterator end (const size_type r);


  template <class StreamType>
  void print (StreamType &out,
              const bool  across = false,
              const bool  diagonal_first = true) 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.");
  DeclExceptionMsg (ExcDifferentSparsityPatterns,
                    "When copying one sparse matrix into another, "
                    "or when adding one sparse matrix to another, "
                    "both matrices need to refer to the same "
                    "sparsity pattern.");
  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.");

protected:
  void prepare_add();

  void prepare_set();

private:
  SmartPointer<const SparsityPattern,SparseMatrix<number> > cols;

  std::unique_ptr<number[]> val;

  std::size_t max_len;

  // make all other sparse matrices friends
  template <typename somenumber> friend class SparseMatrix;
  template <typename somenumber> friend class SparseLUDecomposition;
  template <typename> friend class SparseILU;

  template <typename> friend class BlockMatrixBase;

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

#ifdef DEAL_II_WITH_MPI

  template <typename Number> friend void Utilities::MPI::sum (const SparseMatrix<Number> &, const MPI_Comm &, SparseMatrix<Number> &);
#endif
};

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



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


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


// Inline the set() and add() functions, since they will be called frequently.
template <typename number>
inline
void
SparseMatrix<number>::set (const size_type i,
                           const size_type j,
                           const number       value)
{
  AssertIsFinite(value);

  const size_type index = cols->operator()(i, j);

  // 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
  if (index == SparsityPattern::invalid_entry)
    {
      Assert ((index != SparsityPattern::invalid_entry) ||
              (value == number()),
              ExcInvalidIndex(i, j));
      return;
    }

  val[index] = value;
}



template <typename number>
template <typename number2>
inline
void
SparseMatrix<number>::set (const std::vector<size_type> &indices,
                           const FullMatrix<number2>    &values,
                           const bool                    elide_zero_values)
{
  Assert (indices.size() == values.m(),
          ExcDimensionMismatch(indices.size(), values.m()));
  Assert (values.m() == values.n(), ExcNotQuadratic());

  for (size_type i=0; i<indices.size(); ++i)
    set (indices[i], indices.size(), indices.data(), &values(i,0),
         elide_zero_values);
}



template <typename number>
template <typename number2>
inline
void
SparseMatrix<number>::set (const std::vector<size_type> &row_indices,
                           const std::vector<size_type> &col_indices,
                           const FullMatrix<number2>    &values,
                           const bool                    elide_zero_values)
{
  Assert (row_indices.size() == values.m(),
          ExcDimensionMismatch(row_indices.size(), values.m()));
  Assert (col_indices.size() == values.n(),
          ExcDimensionMismatch(col_indices.size(), values.n()));

  for (size_type i=0; i<row_indices.size(); ++i)
    set (row_indices[i], col_indices.size(), col_indices.data(), &values(i,0),
         elide_zero_values);
}



template <typename number>
template <typename number2>
inline
void
SparseMatrix<number>::set (const size_type               row,
                           const std::vector<size_type> &col_indices,
                           const std::vector<number2>   &values,
                           const bool                    elide_zero_values)
{
  Assert (col_indices.size() == values.size(),
          ExcDimensionMismatch(col_indices.size(), values.size()));

  set (row, col_indices.size(), col_indices.data(), values.data(),
       elide_zero_values);
}



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

  if (value == number())
    return;

  const size_type index = cols->operator()(i, j);

  // it is allowed to add elements to the matrix that are not part of the
  // sparsity pattern, if the value to which we set it is zero
  if (index == SparsityPattern::invalid_entry)
    {
      Assert ((index != SparsityPattern::invalid_entry) ||
              (value == number()),
              ExcInvalidIndex(i, j));
      return;
    }

  val[index] += value;
}



template <typename number>
template <typename number2>
inline
void
SparseMatrix<number>::add (const std::vector<size_type> &indices,
                           const FullMatrix<number2>    &values,
                           const bool                    elide_zero_values)
{
  Assert (indices.size() == values.m(),
          ExcDimensionMismatch(indices.size(), values.m()));
  Assert (values.m() == values.n(), ExcNotQuadratic());

  for (size_type i=0; i<indices.size(); ++i)
    add (indices[i], indices.size(), indices.data(), &values(i,0),
         elide_zero_values);
}



template <typename number>
template <typename number2>
inline
void
SparseMatrix<number>::add (const std::vector<size_type> &row_indices,
                           const std::vector<size_type> &col_indices,
                           const FullMatrix<number2>    &values,
                           const bool                    elide_zero_values)
{
  Assert (row_indices.size() == values.m(),
          ExcDimensionMismatch(row_indices.size(), values.m()));
  Assert (col_indices.size() == values.n(),
          ExcDimensionMismatch(col_indices.size(), values.n()));

  for (size_type i=0; i<row_indices.size(); ++i)
    add (row_indices[i], col_indices.size(), col_indices.data(), &values(i,0),
         elide_zero_values);
}



template <typename number>
template <typename number2>
inline
void
SparseMatrix<number>::add (const size_type               row,
                           const std::vector<size_type> &col_indices,
                           const std::vector<number2>   &values,
                           const bool                    elide_zero_values)
{
  Assert (col_indices.size() == values.size(),
          ExcDimensionMismatch(col_indices.size(), values.size()));

  add (row, col_indices.size(), col_indices.data(), values.data(),
       elide_zero_values);
}



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

  number             *val_ptr    = val.get();
  const number *const end_ptr    = val.get() + cols->n_nonzero_elements();

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

  return *this;
}



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

  const number factor_inv = number(1.) / factor;

  number             *val_ptr    = val.get();
  const number *const end_ptr    = val.get() + cols->n_nonzero_elements();

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

  return *this;
}



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



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



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

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



template <typename number>
inline
number SparseMatrix<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
  return val[cols->rowstart[i]];
}



template <typename number>
inline
number &SparseMatrix<number>::diag_element (const size_type i)
{
  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
  return val[cols->rowstart[i]];
}



template <typename number>
template <typename ForwardIterator>
void
SparseMatrix<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 SparseMatrixIterators
{
  template <typename number>
  inline
  Accessor<number,true>::
  Accessor (const MatrixType   *matrix,
            const std::size_t   index_within_matrix)
    :
    SparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern(),
                                        index_within_matrix),
    matrix (matrix)
  {}



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



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



  template <typename number>
  inline
  number
  Accessor<number, true>::value () const
  {
    AssertIndexRange(index_within_sparsity, matrix->n_nonzero_elements());
    return matrix->val[index_within_sparsity];
  }



  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
  {
    AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
    return accessor->matrix->val[accessor->index_within_sparsity];
  }



  template <typename number>
  inline
  const typename Accessor<number, false>::Reference &
  Accessor<number, false>::Reference::operator = (const number n) const
  {
    AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
    accessor->matrix->val[accessor->index_within_sparsity] = n;
    return *this;
  }



  template <typename number>
  inline
  const typename Accessor<number, false>::Reference &
  Accessor<number, false>::Reference::operator += (const number n) const
  {
    AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
    accessor->matrix->val[accessor->index_within_sparsity] += n;
    return *this;
  }



  template <typename number>
  inline
  const typename Accessor<number, false>::Reference &
  Accessor<number, false>::Reference::operator -= (const number n) const
  {
    AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
    accessor->matrix->val[accessor->index_within_sparsity] -= n;
    return *this;
  }



  template <typename number>
  inline
  const typename Accessor<number, false>::Reference &
  Accessor<number, false>::Reference::operator *= (const number n) const
  {
    AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
    accessor->matrix->val[accessor->index_within_sparsity] *= n;
    return *this;
  }



  template <typename number>
  inline
  const typename Accessor<number, false>::Reference &
  Accessor<number, false>::Reference::operator /= (const number n) const
  {
    AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
    accessor->matrix->val[accessor->index_within_sparsity] /= n;
    return *this;
  }



  template <typename number>
  inline
  Accessor<number,false>::
  Accessor (MatrixType         *matrix,
            const std::size_t   index)
    :
    SparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern(),
                                        index),
    matrix (matrix)
  {}



  template <typename number>
  inline
  Accessor<number,false>::
  Accessor (MatrixType         *matrix)
    :
    SparsityPatternIterators::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 std::size_t  index)
    :
    accessor(matrix, index)
  {}



  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 SparseMatrixIterators::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());

    return (*this)->index_within_sparsity - other->index_within_sparsity;
  }



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

    return x;
  }

}



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


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


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


template <typename number>
inline
typename SparseMatrix<number>::iterator
SparseMatrix<number>::end ()
{
  return iterator(this, cols->rowstart[cols->rows]);
}


template <typename number>
inline
typename SparseMatrix<number>::const_iterator
SparseMatrix<number>::begin (const size_type r) const
{
  Assert (r<m(), ExcIndexRange(r,0,m()));

  return const_iterator(this, cols->rowstart[r]);
}



template <typename number>
inline
typename SparseMatrix<number>::const_iterator
SparseMatrix<number>::end (const size_type r) const
{
  Assert (r<m(), ExcIndexRange(r,0,m()));

  return const_iterator(this, cols->rowstart[r+1]);
}



template <typename number>
inline
typename SparseMatrix<number>::iterator
SparseMatrix<number>::begin (const size_type r)
{
  Assert (r<m(), ExcIndexRange(r,0,m()));

  return iterator(this, cols->rowstart[r]);
}



template <typename number>
inline
typename SparseMatrix<number>::iterator
SparseMatrix<number>::end (const size_type r)
{
  Assert (r<m(), ExcIndexRange(r,0,m()));

  return iterator(this, cols->rowstart[r+1]);
}



template <typename number>
template <class StreamType>
inline
void SparseMatrix<number>::print (StreamType &out,
                                  const bool  across,
                                  const bool  diagonal_first) const
{
  Assert (cols != nullptr, ExcNotInitialized());
  Assert (val != nullptr, ExcNotInitialized());

  bool hanging_diagonal = false;
  number diagonal = number();

  for (size_type i=0; i<cols->rows; ++i)
    {
      for (size_type j=cols->rowstart[i]; j<cols->rowstart[i+1]; ++j)
        {
          if (!diagonal_first && i == cols->colnums[j])
            {
              diagonal = val[j];
              hanging_diagonal = true;
            }
          else
            {
              if (hanging_diagonal && cols->colnums[j]>i)
                {
                  if (across)
                    out << ' ' << i << ',' << i << ':' << diagonal;
                  else
                    out << '(' << i << ',' << i << ") " << diagonal << std::endl;
                  hanging_diagonal = false;
                }
              if (across)
                out << ' ' << i << ',' << cols->colnums[j] << ':' << val[j];
              else
                out << "(" << i << "," << cols->colnums[j] << ") " << val[j] << std::endl;
            }
        }
      if (hanging_diagonal)
        {
          if (across)
            out << ' ' << i << ',' << i << ':' << diagonal;
          else
            out << '(' << i << ',' << i << ") " << diagonal << std::endl;
          hanging_diagonal = false;
        }
    }
  if (across)
    out << std::endl;
}


template <typename number>
inline
void
SparseMatrix<number>::
prepare_add()
{
  //nothing to do here
}



template <typename number>
inline
void
SparseMatrix<number>::
prepare_set()
{
  //nothing to do here
}

#endif // DOXYGEN


/*----------------------------   sparse_matrix.h     ---------------------------*/

DEAL_II_NAMESPACE_CLOSE

#endif
/*----------------------------   sparse_matrix.h     ---------------------------*/