full_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.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------

#ifndef dealii_full_matrix_h
#define dealii_full_matrix_h


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

#include <deal.II/base/numbers.h>
#include <deal.II/base/table.h>
#include <deal.II/base/tensor.h>

#include <deal.II/differentiation/ad/ad_number_traits.h>

#include <deal.II/lac/exceptions.h>
#include <deal.II/lac/identity_matrix.h>

#include <cstring>
#include <iomanip>
#include <vector>

DEAL_II_NAMESPACE_OPEN


// forward declarations
template <typename number>
class Vector;
template <typename number>
class LAPACKFullMatrix;


template <typename number>
class FullMatrix : public Table<2, number>
{
public:
  // The assertion in full_matrix.templates.h for whether or not a number is
  // finite is not compatible for AD number types.
  static_assert(
    !Differentiation::AD::is_ad_number<number>::value,
    "The FullMatrix class does not support auto-differentiable numbers.");

  using size_type = std::size_t;

  using value_type = number;


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


  class const_iterator;

  class Accessor
  {
  public:
    Accessor(const FullMatrix<number> *matrix,
             const size_type           row,
             const size_type           col);

    size_type
    row() const;

    size_type
    column() const;

    number
    value() const;

  protected:
    const FullMatrix<number> *matrix;

    size_type a_row;

    size_type a_col;

    /*
     * Make enclosing class a friend.
     */
    friend class const_iterator;
  };

  class const_iterator
  {
  public:
    const_iterator(const FullMatrix<number> *matrix,
                   const size_type           row,
                   const size_type           col);

    const_iterator &
    operator++();

    const_iterator
    operator++(int);

    const Accessor &operator*() const;

    const Accessor *operator->() const;

    bool
    operator==(const const_iterator &) const;
    bool
    operator!=(const const_iterator &) const;

    bool
    operator<(const const_iterator &) const;

    bool
    operator>(const const_iterator &) const;

  private:
    Accessor accessor;
  };

  explicit FullMatrix(const size_type n = 0);

  FullMatrix(const size_type rows, const size_type cols);

  FullMatrix(const size_type rows, const size_type cols, const number *entries);

  FullMatrix(const IdentityMatrix &id);
  template <typename number2>
  FullMatrix<number> &
  operator=(const FullMatrix<number2> &);

  FullMatrix<number> &
  operator=(const number d);

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

  template <typename number2>
  FullMatrix<number> &
  operator=(const LAPACKFullMatrix<number2> &);


  template <typename MatrixType>
  void
  copy_from(const MatrixType &);

  template <typename MatrixType>
  void
  copy_transposed(const MatrixType &);

  template <int dim>
  void
  copy_from(const Tensor<2, dim> &T,
            const unsigned int    src_r_i = 0,
            const unsigned int    src_r_j = dim - 1,
            const unsigned int    src_c_i = 0,
            const unsigned int    src_c_j = dim - 1,
            const size_type       dst_r   = 0,
            const size_type       dst_c   = 0);

  template <int dim>
  void copy_to(Tensor<2, dim> &   T,
               const size_type    src_r_i = 0,
               const size_type    src_r_j = dim - 1,
               const size_type    src_c_i = 0,
               const size_type    src_c_j = dim - 1,
               const unsigned int dst_r   = 0,
               const unsigned int dst_c   = 0) const;

  template <typename MatrixType, typename index_type>
  void
  extract_submatrix_from(const MatrixType &             matrix,
                         const std::vector<index_type> &row_index_set,
                         const std::vector<index_type> &column_index_set);

  template <typename MatrixType, typename index_type>
  void
  scatter_matrix_to(const std::vector<index_type> &row_index_set,
                    const std::vector<index_type> &column_index_set,
                    MatrixType &                   matrix) const;

  template <typename number2>
  void
  fill(const FullMatrix<number2> &src,
       const size_type            dst_offset_i = 0,
       const size_type            dst_offset_j = 0,
       const size_type            src_offset_i = 0,
       const size_type            src_offset_j = 0);


  template <typename number2>
  void
  fill(const number2 *);

  template <typename number2>
  void
  fill_permutation(const FullMatrix<number2> &   src,
                   const std::vector<size_type> &p_rows,
                   const std::vector<size_type> &p_cols);

  void
  set(const size_type i, const size_type j, const number value);
  bool
  operator==(const FullMatrix<number> &) const;

  size_type
  m() const;

  size_type
  n() const;

  bool
  all_zero() const;

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

  template <typename number2>
  number2
  matrix_scalar_product(const Vector<number2> &u,
                        const Vector<number2> &v) const;

  real_type
  l1_norm() const;

  real_type
  linfty_norm() const;

  real_type
  frobenius_norm() const;

  real_type
  relative_symmetry_norm2() const;

  number
  determinant() const;

  number
  trace() const;

  template <class StreamType>
  void
  print(StreamType &       s,
        const unsigned int width     = 5,
        const unsigned int precision = 2) 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 double       threshold   = 0.) const;

  std::size_t
  memory_consumption() const;



  const_iterator
  begin() const;

  const_iterator
  end() const;

  const_iterator
  begin(const size_type r) const;

  const_iterator
  end(const size_type r) const;



  FullMatrix &
  operator*=(const number factor);

  FullMatrix &
  operator/=(const number factor);

  template <typename number2>
  void
  add(const number a, const FullMatrix<number2> &A);

  template <typename number2>
  void
  add(const number               a,
      const FullMatrix<number2> &A,
      const number               b,
      const FullMatrix<number2> &B);

  template <typename number2>
  void
  add(const number               a,
      const FullMatrix<number2> &A,
      const number               b,
      const FullMatrix<number2> &B,
      const number               c,
      const FullMatrix<number2> &C);

  template <typename number2>
  void
  add(const FullMatrix<number2> &src,
      const number               factor,
      const size_type            dst_offset_i = 0,
      const size_type            dst_offset_j = 0,
      const size_type            src_offset_i = 0,
      const size_type            src_offset_j = 0);

  template <typename number2>
  void
  Tadd(const number s, const FullMatrix<number2> &B);

  template <typename number2>
  void
  Tadd(const FullMatrix<number2> &src,
       const number               factor,
       const size_type            dst_offset_i = 0,
       const size_type            dst_offset_j = 0,
       const size_type            src_offset_i = 0,
       const size_type            src_offset_j = 0);

  void
  add(const size_type row, const size_type column, const number value);

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

  void
  add_row(const size_type i, const number s, const size_type j);

  void
  add_row(const size_type i,
          const number    s,
          const size_type j,
          const number    t,
          const size_type k);

  void
  add_col(const size_type i, const number s, const size_type j);

  void
  add_col(const size_type i,
          const number    s,
          const size_type j,
          const number    t,
          const size_type k);

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

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

  void
  diagadd(const number s);

  template <typename number2>
  void
  equ(const number a, const FullMatrix<number2> &A);

  template <typename number2>
  void
  equ(const number               a,
      const FullMatrix<number2> &A,
      const number               b,
      const FullMatrix<number2> &B);

  template <typename number2>
  void
  equ(const number               a,
      const FullMatrix<number2> &A,
      const number               b,
      const FullMatrix<number2> &B,
      const number               c,
      const FullMatrix<number2> &C);

  void
  symmetrize();

  void
  gauss_jordan();

  template <typename number2>
  void
  invert(const FullMatrix<number2> &M);

  template <typename number2>
  void
  cholesky(const FullMatrix<number2> &A);

  template <typename number2>
  void
  outer_product(const Vector<number2> &V, const Vector<number2> &W);

  template <typename number2>
  void
  left_invert(const FullMatrix<number2> &M);

  template <typename number2>
  void
  right_invert(const FullMatrix<number2> &M);



  template <typename number2>
  void
  mmult(FullMatrix<number2> &      C,
        const FullMatrix<number2> &B,
        const bool                 adding = false) const;

  template <typename number2>
  void
  Tmmult(FullMatrix<number2> &      C,
         const FullMatrix<number2> &B,
         const bool                 adding = false) const;

  template <typename number2>
  void
  mTmult(FullMatrix<number2> &      C,
         const FullMatrix<number2> &B,
         const bool                 adding = false) const;

  template <typename number2>
  void
  TmTmult(FullMatrix<number2> &      C,
          const FullMatrix<number2> &B,
          const bool                 adding = false) const;

  void
  triple_product(const FullMatrix<number> &A,
                 const FullMatrix<number> &B,
                 const FullMatrix<number> &D,
                 const bool                transpose_B = false,
                 const bool                transpose_D = false,
                 const number              scaling     = number(1.));

  template <typename number2>
  void
  vmult(Vector<number2> &      w,
        const Vector<number2> &v,
        const bool             adding = false) const;

  template <typename number2>
  void
  vmult_add(Vector<number2> &w, const Vector<number2> &v) const;

  template <typename number2>
  void
  Tvmult(Vector<number2> &      w,
         const Vector<number2> &v,
         const bool             adding = false) const;

  template <typename number2>
  void
  Tvmult_add(Vector<number2> &w, const Vector<number2> &v) const;

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

  template <typename number2, typename number3>
  number
  residual(Vector<number2> &      dst,
           const Vector<number2> &x,
           const Vector<number3> &b) const;

  template <typename number2>
  void
  forward(Vector<number2> &dst, const Vector<number2> &src) const;

  template <typename number2>
  void
  backward(Vector<number2> &dst, const Vector<number2> &src) const;


  DeclException0(ExcEmptyMatrix);

  DeclException1(
    ExcNotRegular,
    number,
    << "The maximal pivot is " << arg1
    << ", which is below the threshold. The matrix may be singular.");
  DeclException3(ExcInvalidDestination,
                 size_type,
                 size_type,
                 size_type,
                 << "Target region not in matrix: size in this direction="
                 << arg1 << ", size of new matrix=" << arg2
                 << ", offset=" << arg3);
  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.");
  DeclException0(ExcMatrixNotPositiveDefinite);

  friend class Accessor;
};

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



template <typename number>
inline typename FullMatrix<number>::size_type
FullMatrix<number>::m() const
{
  return this->n_rows();
}



template <typename number>
inline typename FullMatrix<number>::size_type
FullMatrix<number>::n() const
{
  return this->n_cols();
}



template <typename number>
FullMatrix<number> &
FullMatrix<number>::operator=(const number d)
{
  Assert(d == number(0), ExcScalarAssignmentOnlyForZeroValue());
  (void)d; // removes -Wunused-parameter warning in optimized mode

  if (this->n_elements() != 0)
    this->reset_values();

  return *this;
}



template <typename number>
template <typename number2>
inline void
FullMatrix<number>::fill(const number2 *src)
{
  Table<2, number>::fill(src);
}



template <typename number>
template <typename MatrixType>
void
FullMatrix<number>::copy_from(const MatrixType &M)
{
  this->reinit(M.m(), M.n());

  // loop over the elements of the argument matrix row by row, as suggested
  // in the documentation of the sparse matrix iterator class, and
  // copy them into the current object
  for (size_type row = 0; row < M.m(); ++row)
    {
      const typename MatrixType::const_iterator end_row = M.end(row);
      for (typename MatrixType::const_iterator entry = M.begin(row);
           entry != end_row;
           ++entry)
        this->el(row, entry->column()) = entry->value();
    }
}



template <typename number>
template <typename MatrixType>
void
FullMatrix<number>::copy_transposed(const MatrixType &M)
{
  this->reinit(M.n(), M.m());

  // loop over the elements of the argument matrix row by row, as suggested
  // in the documentation of the sparse matrix iterator class, and
  // copy them into the current object
  for (size_type row = 0; row < M.m(); ++row)
    {
      const typename MatrixType::const_iterator end_row = M.end(row);
      for (typename MatrixType::const_iterator entry = M.begin(row);
           entry != end_row;
           ++entry)
        this->el(entry->column(), row) = entry->value();
    }
}



template <typename number>
template <typename MatrixType, typename index_type>
inline void
FullMatrix<number>::extract_submatrix_from(
  const MatrixType &             matrix,
  const std::vector<index_type> &row_index_set,
  const std::vector<index_type> &column_index_set)
{
  AssertDimension(row_index_set.size(), this->n_rows());
  AssertDimension(column_index_set.size(), this->n_cols());

  const size_type n_rows_submatrix = row_index_set.size();
  const size_type n_cols_submatrix = column_index_set.size();

  for (size_type sub_row = 0; sub_row < n_rows_submatrix; ++sub_row)
    for (size_type sub_col = 0; sub_col < n_cols_submatrix; ++sub_col)
      (*this)(sub_row, sub_col) =
        matrix.el(row_index_set[sub_row], column_index_set[sub_col]);
}



template <typename number>
template <typename MatrixType, typename index_type>
inline void
FullMatrix<number>::scatter_matrix_to(
  const std::vector<index_type> &row_index_set,
  const std::vector<index_type> &column_index_set,
  MatrixType &                   matrix) const
{
  AssertDimension(row_index_set.size(), this->n_rows());
  AssertDimension(column_index_set.size(), this->n_cols());

  const size_type n_rows_submatrix = row_index_set.size();
  const size_type n_cols_submatrix = column_index_set.size();

  for (size_type sub_row = 0; sub_row < n_rows_submatrix; ++sub_row)
    for (size_type sub_col = 0; sub_col < n_cols_submatrix; ++sub_col)
      matrix.set(row_index_set[sub_row],
                 column_index_set[sub_col],
                 (*this)(sub_row, sub_col));
}


template <typename number>
inline void
FullMatrix<number>::set(const size_type i,
                        const size_type j,
                        const number    value)
{
  (*this)(i, j) = value;
}



template <typename number>
template <typename number2>
void
FullMatrix<number>::vmult_add(Vector<number2> &      w,
                              const Vector<number2> &v) const
{
  vmult(w, v, true);
}


template <typename number>
template <typename number2>
void
FullMatrix<number>::Tvmult_add(Vector<number2> &      w,
                               const Vector<number2> &v) const
{
  Tvmult(w, v, true);
}


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


template <typename number>
inline FullMatrix<number>::Accessor::Accessor(const FullMatrix<number> *matrix,
                                              const size_type           r,
                                              const size_type           c)
  : matrix(matrix)
  , a_row(r)
  , a_col(c)
{}


template <typename number>
inline typename FullMatrix<number>::size_type
FullMatrix<number>::Accessor::row() const
{
  return a_row;
}


template <typename number>
inline typename FullMatrix<number>::size_type
FullMatrix<number>::Accessor::column() const
{
  return a_col;
}


template <typename number>
inline number
FullMatrix<number>::Accessor::value() const
{
  AssertIsFinite(matrix->el(a_row, a_col));
  return matrix->el(a_row, a_col);
}


template <typename number>
inline FullMatrix<number>::const_iterator::const_iterator(
  const FullMatrix<number> *matrix,
  const size_type           r,
  const size_type           c)
  : accessor(matrix, r, c)
{}


template <typename number>
inline typename FullMatrix<number>::const_iterator &
FullMatrix<number>::const_iterator::operator++()
{
  Assert(accessor.a_row < accessor.matrix->m(), ExcIteratorPastEnd());

  ++accessor.a_col;
  if (accessor.a_col >= accessor.matrix->n())
    {
      accessor.a_col = 0;
      accessor.a_row++;
    }
  return *this;
}


template <typename number>
inline typename FullMatrix<number>::const_iterator
FullMatrix<number>::const_iterator::operator++(int)
{
  const typename FullMatrix<number>::const_iterator current = *this;
  ++(*this);

  return current;
}


template <typename number>
inline const typename FullMatrix<number>::Accessor &
  FullMatrix<number>::const_iterator::operator*() const
{
  return accessor;
}


template <typename number>
inline const typename FullMatrix<number>::Accessor *
  FullMatrix<number>::const_iterator::operator->() const
{
  return &accessor;
}


template <typename number>
inline bool
FullMatrix<number>::const_iterator::
operator==(const const_iterator &other) const
{
  return (accessor.row() == other.accessor.row() &&
          accessor.column() == other.accessor.column());
}


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


template <typename number>
inline bool
FullMatrix<number>::const_iterator::operator<(const const_iterator &other) const
{
  return (accessor.row() < other.accessor.row() ||
          (accessor.row() == other.accessor.row() &&
           accessor.column() < other.accessor.column()));
}


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


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


template <typename number>
inline typename FullMatrix<number>::const_iterator
FullMatrix<number>::end() const
{
  return const_iterator(this, m(), 0);
}


template <typename number>
inline typename FullMatrix<number>::const_iterator
FullMatrix<number>::begin(const size_type r) const
{
  AssertIndexRange(r, m());
  return const_iterator(this, r, 0);
}



template <typename number>
inline typename FullMatrix<number>::const_iterator
FullMatrix<number>::end(const size_type r) const
{
  AssertIndexRange(r, m());
  return const_iterator(this, r + 1, 0);
}



template <typename number>
inline void
FullMatrix<number>::add(const size_type r, const size_type c, const number v)
{
  AssertIndexRange(r, this->m());
  AssertIndexRange(c, this->n());

  this->operator()(r, c) += v;
}



template <typename number>
template <typename number2, typename index_type>
inline void
FullMatrix<number>::add(const size_type   row,
                        const size_type   n_cols,
                        const index_type *col_indices,
                        const number2 *   values,
                        const bool,
                        const bool)
{
  AssertIndexRange(row, this->m());
  for (size_type col = 0; col < n_cols; ++col)
    {
      AssertIndexRange(col_indices[col], this->n());
      this->operator()(row, col_indices[col]) += values[col];
    }
}


template <typename number>
template <class StreamType>
inline void
FullMatrix<number>::print(StreamType &       s,
                          const unsigned int w,
                          const unsigned int p) const
{
  Assert(!this->empty(), ExcEmptyMatrix());

  // save the state of out stream
  const std::streamsize old_precision = s.precision(p);
  const std::streamsize old_width     = s.width(w);

  for (size_type i = 0; i < this->m(); ++i)
    {
      for (size_type j = 0; j < this->n(); ++j)
        {
          s.width(w);
          s.precision(p);
          s << this->el(i, j);
        }
      s << std::endl;
    }

  // reset output format
  s.precision(old_precision);
  s.width(old_width);
}


#endif // DOXYGEN

DEAL_II_NAMESPACE_CLOSE

#endif