sparse_matrix.h

Go to the documentation of this file.

// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 2002 - 2023 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Part of the source code is dual licensed under Apache-2.0 WITH
// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
// governing the source code and code contributions can be found in
// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
//
// ------------------------------------------------------------------------
 
#ifndef dealii_sparse_matrix_h
#define dealii_sparse_matrix_h
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/enable_observer_pointer.h>
#include <deal.II/base/mpi_stub.h>
#include <deal.II/base/observer_pointer.h>
 
#include <deal.II/lac/exceptions.h>
#include <deal.II/lac/identity_matrix.h>
#include <deal.II/lac/sparsity_pattern.h>
#include <deal.II/lac/vector_operation.h>
 
#include <iterator>
#include <memory>
 
 
DEAL_II_NAMESPACE_OPEN
 
// Forward declarations
#ifndef DOXYGEN
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> &);
  }
} // namespace Utilities
#  endif
 
#  ifdef DEAL_II_WITH_TRILINOS
namespace TrilinosWrappers
{
  class SparseMatrix;
}
#  endif
#endif
 
namespace SparseMatrixIterators
{
  using size_type = types::global_dof_index;
 
  // 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:
    using MatrixType = const SparseMatrix<number>;
 
    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;
 
    // Make iterator class a friend.
    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:
    using MatrixType = SparseMatrix<number>;
 
    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;
 
    // Make iterator class a friend.
    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> &;
 
    using difference_type = size_type;
 
    Iterator(MatrixType *matrix, const std::size_t index_within_matrix);
 
    Iterator(MatrixType *matrix);
 
    Iterator(const SparseMatrixIterators::Iterator<number, false> &i);
 
    const Iterator<number, Constness> &
    operator=(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;
  };
 
} // namespace SparseMatrixIterators
 
DEAL_II_NAMESPACE_CLOSE
 
namespace std
{
  template <typename number, bool Constness>
  struct iterator_traits<
    ::SparseMatrixIterators::Iterator<number, Constness>>
  {
    using iterator_category = forward_iterator_tag;
    using value_type =
      typename ::SparseMatrixIterators::Iterator<number,
                                                       Constness>::value_type;
    using difference_type = typename ::SparseMatrixIterators::
      Iterator<number, Constness>::difference_type;
  };
} // namespace std
 
DEAL_II_NAMESPACE_OPEN
 
// TODO: Add multithreading to the other vmult functions.
 
template <typename number>
class SparseMatrix : public virtual EnableObserverPointer
{
public:
  using size_type = types::global_dof_index;
 
  using value_type = number;
 
  using real_type = typename numbers::NumberTraits<number>::real_type;
 
  using const_iterator = SparseMatrixIterators::Iterator<number, true>;
 
  using iterator = SparseMatrixIterators::Iterator<number, false>;
 
  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() override;
 
  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);
 
  template <typename number2>
  void
  reinit(const SparseMatrix<number2> &sparse_matrix);
 
  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              omega = 1.) const;
 
  template <typename somenumber>
  void
  precondition_TSOR(Vector<somenumber>       &dst,
                    const Vector<somenumber> &src,
                    const number              omega = 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 omega = 1.) const;
 
  template <typename somenumber>
  void
  TSOR(Vector<somenumber> &v, const number omega = 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                  omega = 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                  omega = 1.) const;
 
  template <typename somenumber>
  void
  Jacobi_step(Vector<somenumber>       &v,
              const Vector<somenumber> &b,
              const number              omega = 1.) const;
 
  template <typename somenumber>
  void
  SOR_step(Vector<somenumber>       &v,
           const Vector<somenumber> &b,
           const number              omega = 1.) const;
 
  template <typename somenumber>
  void
  TSOR_step(Vector<somenumber>       &v,
            const Vector<somenumber> &b,
            const number              omega = 1.) const;
 
  template <typename somenumber>
  void
  SSOR_step(Vector<somenumber>       &v,
            const Vector<somenumber> &b,
            const number              omega = 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 <typename 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 char        *separator   = " ") const;
 
  void
  print_pattern(std::ostream &out, const double threshold = 0.) const;
 
  void
  print_as_numpy_arrays(std::ostream      &out,
                        const unsigned int precision = 9) 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 vectors 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:
  ObserverPointer<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;
 
  // To allow it calling private prepare_add() and prepare_set().
  template <typename>
  friend class BlockMatrixBase;
 
  // Also give access to internal details to the iterator/accessor classes.
  template <typename, bool>
  friend class SparseMatrixIterators::Iterator;
  template <typename, bool>
  friend class SparseMatrixIterators::Accessor;
 
#ifdef DEAL_II_WITH_MPI
  // Give access to internal datastructures to perform MPI operations.
  template <typename Number>
  friend void
  Utilities::MPI::sum(const SparseMatrix<Number> &,
                      const MPI_Comm,
                      SparseMatrix<Number> &);
#endif
};
 
#ifndef DOXYGEN
/*---------------------- Inline functions -----------------------------------*/
 
 
 
template <typename number>
template <typename number2>
void
SparseMatrix<number>::reinit(const SparseMatrix<number2> &sparse_matrix)
{
  this->reinit(sparse_matrix.get_sparsity_pattern());
}
 
 
 
template <typename number>
inline typename SparseMatrix<number>::size_type
SparseMatrix<number>::m() const
{
  Assert(cols != nullptr, ExcNeedsSparsityPattern());
  return cols->rows;
}
 
 
template <typename number>
inline typename SparseMatrix<number>::size_type
SparseMatrix<number>::n() const
{
  Assert(cols != nullptr, ExcNeedsSparsityPattern());
  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,
      std::is_sorted(col_indices.begin(), col_indices.end()));
}
 
 
 
template <typename number>
inline SparseMatrix<number> &
SparseMatrix<number>::operator*=(const number factor)
{
  Assert(cols != nullptr, ExcNeedsSparsityPattern());
  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, ExcNeedsSparsityPattern());
  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, ExcNeedsSparsityPattern());
  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, ExcNeedsSparsityPattern());
  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, ExcNeedsSparsityPattern());
  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, ExcNeedsSparsityPattern());
  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, ExcNeedsSparsityPattern());
  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 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 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(linear_index, matrix->n_nonzero_elements());
    return matrix->val[linear_index];
  }
 
 
 
  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->linear_index,
                     accessor->matrix->n_nonzero_elements());
    return accessor->matrix->val[accessor->linear_index];
  }
 
 
 
  template <typename number>
  inline const typename Accessor<number, false>::Reference &
  Accessor<number, false>::Reference::operator=(const number n) const
  {
    AssertIndexRange(accessor->linear_index,
                     accessor->matrix->n_nonzero_elements());
    accessor->matrix->val[accessor->linear_index] = n;
    return *this;
  }
 
 
 
  template <typename number>
  inline const typename Accessor<number, false>::Reference &
  Accessor<number, false>::Reference::operator+=(const number n) const
  {
    AssertIndexRange(accessor->linear_index,
                     accessor->matrix->n_nonzero_elements());
    accessor->matrix->val[accessor->linear_index] += n;
    return *this;
  }
 
 
 
  template <typename number>
  inline const typename Accessor<number, false>::Reference &
  Accessor<number, false>::Reference::operator-=(const number n) const
  {
    AssertIndexRange(accessor->linear_index,
                     accessor->matrix->n_nonzero_elements());
    accessor->matrix->val[accessor->linear_index] -= n;
    return *this;
  }
 
 
 
  template <typename number>
  inline const typename Accessor<number, false>::Reference &
  Accessor<number, false>::Reference::operator*=(const number n) const
  {
    AssertIndexRange(accessor->linear_index,
                     accessor->matrix->n_nonzero_elements());
    accessor->matrix->val[accessor->linear_index] *= n;
    return *this;
  }
 
 
 
  template <typename number>
  inline const typename Accessor<number, false>::Reference &
  Accessor<number, false>::Reference::operator/=(const number n) const
  {
    AssertIndexRange(accessor->linear_index,
                     accessor->matrix->n_nonzero_elements());
    accessor->matrix->val[accessor->linear_index] /= 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 const Iterator<number, Constness> &
  Iterator<number, Constness>::operator=(
    const SparseMatrixIterators::Iterator<number, false> &i)
  {
    accessor = *i;
    return *this;
  }
 
 
 
  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)->linear_index - other->linear_index;
  }
 
 
 
  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;
  }
 
} // namespace SparseMatrixIterators
 
 
 
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
{
  AssertIndexRange(r, 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
{
  AssertIndexRange(r, m());
 
  return const_iterator(this, cols->rowstart[r + 1]);
}
 
 
 
template <typename number>
inline typename SparseMatrix<number>::iterator
SparseMatrix<number>::begin(const size_type r)
{
  AssertIndexRange(r, m());
 
  return iterator(this, cols->rowstart[r]);
}
 
 
 
template <typename number>
inline typename SparseMatrix<number>::iterator
SparseMatrix<number>::end(const size_type r)
{
  AssertIndexRange(r, m());
 
  return iterator(this, cols->rowstart[r + 1]);
}
 
 
 
template <typename number>
template <typename StreamType>
inline void
SparseMatrix<number>::print(StreamType &out,
                            const bool  across,
                            const bool  diagonal_first) const
{
  Assert(cols != nullptr, ExcNeedsSparsityPattern());
  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
 
 
DEAL_II_NAMESPACE_CLOSE
 
#endif