full_matrix.h

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
// Copyright (C) 1999 - 2025 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Detailed license information governing the source code and contributions
// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
//
// -----------------------------------------------------------------------------
 
#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/lac/exceptions.h>
#include <deal.II/lac/identity_matrix.h>
#include <deal.II/lac/vector_operation.h>
 
#include <cstring>
#include <iomanip>
#include <vector>
 
DEAL_II_NAMESPACE_OPEN
 
 
// forward declarations
#ifndef DOXYGEN
template <typename number>
class Vector;
template <typename number>
class LAPACKFullMatrix;
#endif
 
template <typename number>
class FullMatrix : public Table<2, number>
{
public:
  static_assert(
    std::is_arithmetic_v<typename numbers::NumberTraits<number>::real_type>,
    "The FullMatrix class only supports basic numeric types. In particular, it "
    "does not support automatically differentiated numbers.");
 
 
  using size_type = std::size_t;
 
  using value_type = number;
 
  using iterator = typename Table<2, number>::iterator;
 
  using const_iterator = typename Table<2, number>::const_iterator;
 
  using Table<2, number>::begin;
 
  using Table<2, number>::end;
 
  using real_type = typename numbers::NumberTraits<number>::real_type;
 
  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 <typename 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 char        *separator   = " ") const;
 
  std::size_t
  memory_consumption() const;
 
  iterator
  begin(const size_type r);
 
  iterator
  end(const size_type r);
 
  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
  permute(const std::vector<unsigned int> &row_perm,
          const std::vector<unsigned int> &col_perm);
 
  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.));
 
 
  void
  kronecker_product(const FullMatrix<number> &A,
                    const FullMatrix<number> &B,
                    const bool                adding);
 
 
  template <typename... Rest>
  void
  kronecker_product(const FullMatrix<number> &A,
                    const FullMatrix<number> &B,
                    const Rest &...rest);
 
 
  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;
 
  void compress(VectorOperation::values);
 
  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 vectors that "
                   "are the same object, but the operation requires that the "
                   "two objects are in fact different.");
  DeclException0(ExcMatrixNotPositiveDefinite);
};
 
#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());
 
  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 <int dim>
void
FullMatrix<number>::copy_from(const Tensor<2, dim> &T,
                              const unsigned int    src_r_i,
                              const unsigned int    src_r_j,
                              const unsigned int    src_c_i,
                              const unsigned int    src_c_j,
                              const size_type       dst_r,
                              const size_type       dst_c)
{
  Assert(!this->empty(), ExcEmptyMatrix());
  AssertIndexRange(src_r_j - src_r_i, this->m() - dst_r);
  AssertIndexRange(src_c_j - src_c_i, this->n() - dst_c);
  AssertIndexRange(src_r_j, dim);
  AssertIndexRange(src_c_j, dim);
  AssertIndexRange(src_r_i, src_r_j + 1);
  AssertIndexRange(src_c_i, src_c_j + 1);
 
  for (size_type i = 0; i < src_r_j - src_r_i + 1; ++i)
    for (size_type j = 0; j < src_c_j - src_c_i + 1; ++j)
      {
        const unsigned int src_r_index = static_cast<unsigned int>(i + src_r_i);
        const unsigned int src_c_index = static_cast<unsigned int>(j + src_c_i);
        (*this)(i + dst_r, j + dst_c)  = number(T[src_r_index][src_c_index]);
      }
}
 
 
 
template <typename number>
template <int dim>
void
FullMatrix<number>::copy_to(Tensor<2, dim>    &T,
                            const size_type    src_r_i,
                            const size_type    src_r_j,
                            const size_type    src_c_i,
                            const size_type    src_c_j,
                            const unsigned int dst_r,
                            const unsigned int dst_c) const
{
  Assert(!this->empty(), ExcEmptyMatrix());
  AssertIndexRange(src_r_j - src_r_i, dim - dst_r);
  AssertIndexRange(src_c_j - src_c_i, dim - dst_c);
  AssertIndexRange(src_r_j, this->m());
  AssertIndexRange(src_r_j, this->n());
  AssertIndexRange(src_r_i, src_r_j + 1);
  AssertIndexRange(src_c_j, src_c_j + 1);
 
  for (size_type i = 0; i < src_r_j - src_r_i + 1; ++i)
    for (size_type j = 0; j < src_c_j - src_c_i + 1; ++j)
      {
        const unsigned int dst_r_index = static_cast<unsigned int>(i + dst_r);
        const unsigned int dst_c_index = static_cast<unsigned int>(j + dst_c);
        T[dst_r_index][dst_c_index] = double((*this)(i + src_r_i, j + src_c_i));
      }
}
 
 
 
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 typename FullMatrix<number>::iterator
FullMatrix<number>::begin(const size_type r)
{
  AssertIndexRange(r, m());
  return iterator(this, r, 0);
}
 
 
 
template <typename number>
inline typename FullMatrix<number>::iterator
FullMatrix<number>::end(const size_type r)
{
  AssertIndexRange(r, m());
  return iterator(this, r + 1, 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 <typename 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);
}
 
 
template <typename number>
template <typename... Rest>
void
FullMatrix<number>::kronecker_product(const FullMatrix<number> &A,
                                      const FullMatrix<number> &B,
                                      const Rest &...rest)
{
  this->kronecker_product(A, B, false); // tmp := A ⊗ B
 
  static_assert(
    (std::is_same_v<std::remove_cv_t<std::remove_reference_t<Rest>>,
                    FullMatrix<number>> &&
     ...),
    "kronecker_product: all variadic arguments must be ::FullMatrix<number>");
 
  if constexpr (sizeof...(rest) != 0)
    {
      // copy the current object into tmp and proceed with the next factor
      FullMatrix<number> tmp = *this;
      this->kronecker_product(tmp, rest...);
    }
}
 
 
 
#endif // DOXYGEN
 
DEAL_II_NAMESPACE_CLOSE
 
#endif