full_matrix.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1999 - 2017 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#ifndef dealii__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/lac/exceptions.h>
#include <deal.II/lac/identity_matrix.h>
#include <deal.II/base/tensor.h>

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

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:
  typedef unsigned int size_type;

  typedef number value_type;


  typedef typename numbers::NumberTraits<number>::real_type 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;

    unsigned short 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 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 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 size_type dst_r=0,
          const size_type 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 unsigned int  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);
  DeclException0 (ExcSourceEqualsDestination);
  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
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
{
  Assert (r<m(), ExcIndexRange(r,0,m()));
  return const_iterator(this, r, 0);
}



template <typename number>
inline
typename FullMatrix<number>::const_iterator
FullMatrix<number>::end (const size_type r) const
{
  Assert (r<m(), ExcIndexRange(r,0,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 unsigned int 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 unsigned int old_precision = s.precision (p);
  const unsigned int 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