lapack_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) 2005 - 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_lapack_full_matrix_h
#define dealii_lapack_full_matrix_h
 
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/mutex.h>
#include <deal.II/base/observer_pointer.h>
#include <deal.II/base/table.h>
 
#include <deal.II/lac/lapack_support.h>
#include <deal.II/lac/vector_memory.h>
 
#include <complex>
#include <memory>
#include <vector>
 
DEAL_II_NAMESPACE_OPEN
 
// forward declarations
#ifndef DOXYGEN
template <typename number>
class Vector;
template <typename number>
class BlockVector;
template <typename number>
class FullMatrix;
template <typename number>
class SparseMatrix;
#endif
 
template <typename number>
class LAPACKFullMatrix : public TransposeTable<number>
{
public:
  using size_type = std::make_unsigned_t<types::blas_int>;
 
  explicit LAPACKFullMatrix(const size_type size = 0);
 
  LAPACKFullMatrix(const size_type rows, const size_type cols);
 
  LAPACKFullMatrix(const LAPACKFullMatrix &);
 
  LAPACKFullMatrix<number> &
  operator=(const LAPACKFullMatrix<number> &);
 
  template <typename number2>
  LAPACKFullMatrix<number> &
  operator=(const FullMatrix<number2> &);
 
  template <typename number2>
  LAPACKFullMatrix<number> &
  operator=(const SparseMatrix<number2> &);
 
  LAPACKFullMatrix<number> &
  operator=(const number d);
 
  LAPACKFullMatrix<number> &
  operator*=(const number factor);
 
  LAPACKFullMatrix<number> &
  operator/=(const number factor);
 
  void
  set(const size_type i, const size_type j, const number value);
 
  void
  add(const number a, const LAPACKFullMatrix<number> &B);
 
  void
  rank1_update(const number a, const Vector<number> &v);
 
  void
  apply_givens_rotation(const std::array<number, 3> &csr,
                        const size_type              i,
                        const size_type              k,
                        const bool                   left = true);
 
  template <typename MatrixType>
  void
  copy_from(const MatrixType &);
 
  void
  reinit(const size_type size);
 
  void
  grow_or_shrink(const size_type size);
 
  void
  remove_row_and_column(const size_type row, const size_type col);
 
  void
  reinit(const size_type rows, const size_type cols);
 
  void
  set_property(const LAPACKSupport::Property property);
 
  size_type
  m() const;
 
  size_type
  n() const;
 
  template <typename MatrixType>
  void
  fill(const MatrixType &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,
       const number      factor       = 1.,
       const bool        transpose    = false);
 
 
  template <typename number2>
  void
  vmult(Vector<number2>       &w,
        const Vector<number2> &v,
        const bool             adding = false) const;
 
  void
  vmult(Vector<number>       &w,
        const Vector<number> &v,
        const bool            adding = false) const;
 
  template <typename number2>
  void
  vmult_add(Vector<number2> &w, const Vector<number2> &v) const;
 
  void
  vmult_add(Vector<number> &w, const Vector<number> &v) const;
 
  template <typename number2>
  void
  Tvmult(Vector<number2>       &w,
         const Vector<number2> &v,
         const bool             adding = false) const;
 
  void
  Tvmult(Vector<number>       &w,
         const Vector<number> &v,
         const bool            adding = false) const;
 
  template <typename number2>
  void
  Tvmult_add(Vector<number2> &w, const Vector<number2> &v) const;
 
  void
  Tvmult_add(Vector<number> &w, const Vector<number> &v) const;
 
 
  void
  mmult(LAPACKFullMatrix<number>       &C,
        const LAPACKFullMatrix<number> &B,
        const bool                      adding = false) const;
 
  void
  mmult(FullMatrix<number>             &C,
        const LAPACKFullMatrix<number> &B,
        const bool                      adding = false) const;
 
  void
  Tmmult(LAPACKFullMatrix<number>       &C,
         const LAPACKFullMatrix<number> &B,
         const bool                      adding = false) const;
 
  void
  Tmmult(FullMatrix<number>             &C,
         const LAPACKFullMatrix<number> &B,
         const bool                      adding = false) const;
 
  void
  Tmmult(LAPACKFullMatrix<number>       &C,
         const LAPACKFullMatrix<number> &B,
         const Vector<number>           &V,
         const bool                      adding = false) const;
 
  void
  mTmult(LAPACKFullMatrix<number>       &C,
         const LAPACKFullMatrix<number> &B,
         const bool                      adding = false) const;
 
  void
  mTmult(FullMatrix<number>             &C,
         const LAPACKFullMatrix<number> &B,
         const bool                      adding = false) const;
 
  void
  TmTmult(LAPACKFullMatrix<number>       &C,
          const LAPACKFullMatrix<number> &B,
          const bool                      adding = false) const;
 
  void
  TmTmult(FullMatrix<number>             &C,
          const LAPACKFullMatrix<number> &B,
          const bool                      adding = false) const;
 
  void
  transpose(LAPACKFullMatrix<number> &B) const;
 
  void
  scale_rows(const Vector<number> &V);
 
  void
  compute_lu_factorization();
 
  void
  compute_cholesky_factorization();
 
  number
  reciprocal_condition_number(const number l1_norm) const;
 
  number
  reciprocal_condition_number() const;
 
  number
  determinant() const;
 
  number
  l1_norm() const;
 
  number
  linfty_norm() const;
 
  number
  frobenius_norm() const;
 
  number
  trace() const;
 
  void
  invert();
 
  void
  solve(Vector<number> &v, const bool transposed = false) const;
 
  void
  solve(LAPACKFullMatrix<number> &B, const bool transposed = false) const;
 
  void
  compute_eigenvalues(const bool right_eigenvectors = false,
                      const bool left_eigenvectors  = false);
 
  void
  compute_eigenvalues_symmetric(const number        lower_bound,
                                const number        upper_bound,
                                const number        abs_accuracy,
                                Vector<number>     &eigenvalues,
                                FullMatrix<number> &eigenvectors);
 
  void
  compute_generalized_eigenvalues_symmetric(
    LAPACKFullMatrix<number>    &B,
    const number                 lower_bound,
    const number                 upper_bound,
    const number                 abs_accuracy,
    Vector<number>              &eigenvalues,
    std::vector<Vector<number>> &eigenvectors,
    const types::blas_int        itype = 1);
 
  void
  compute_generalized_eigenvalues_symmetric(
    LAPACKFullMatrix<number>    &B,
    std::vector<Vector<number>> &eigenvectors,
    const types::blas_int        itype = 1);
 
  void
  compute_svd();
 
  void
  compute_inverse_svd(const double threshold = 0.);
 
  void
  compute_inverse_svd_with_kernel(const unsigned int kernel_size);
 
  std::complex<typename numbers::NumberTraits<number>::real_type>
  eigenvalue(const size_type i) const;
 
  FullMatrix<std::complex<typename numbers::NumberTraits<number>::real_type>>
  get_right_eigenvectors() const;
 
  FullMatrix<std::complex<typename numbers::NumberTraits<number>::real_type>>
  get_left_eigenvectors() const;
 
  number
  singular_value(const size_type i) const;
 
  inline const LAPACKFullMatrix<number> &
  get_svd_u() const;
 
  inline const LAPACKFullMatrix<number> &
  get_svd_vt() 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;
 
  LAPACKSupport::State
  get_state() const;
 
private:
  number
  norm(const char type) const;
 
  LAPACKSupport::State state;
 
  LAPACKSupport::Property property;
 
  mutable std::vector<number> work;
 
  mutable std::vector<types::blas_int> iwork;
 
  std::vector<types::blas_int> ipiv;
 
  std::vector<number> inv_work;
 
  std::vector<typename numbers::NumberTraits<number>::real_type> wr;
 
  std::vector<number> wi;
 
  std::vector<number> vl;
 
  std::vector<number> vr;
 
  std::unique_ptr<LAPACKFullMatrix<number>> svd_u;
 
  std::unique_ptr<LAPACKFullMatrix<number>> svd_vt;
 
  mutable Threads::Mutex mutex;
};
 
 
 
template <typename number>
class PreconditionLU : public EnableObserverPointer
{
public:
  void
  initialize(const LAPACKFullMatrix<number> &);
  void
  initialize(const LAPACKFullMatrix<number> &, VectorMemory<Vector<number>> &);
  void
  vmult(Vector<number> &, const Vector<number> &) const;
  void
  Tvmult(Vector<number> &, const Vector<number> &) const;
  void
  vmult(BlockVector<number> &, const BlockVector<number> &) const;
  void
  Tvmult(BlockVector<number> &, const BlockVector<number> &) const;
 
private:
  ObserverPointer<const LAPACKFullMatrix<number>, PreconditionLU<number>>
                                                                        matrix;
  ObserverPointer<VectorMemory<Vector<number>>, PreconditionLU<number>> mem;
};
 
/*---------------------- Inline functions -----------------------------------*/
 
template <typename number>
inline void
LAPACKFullMatrix<number>::set(const size_type i,
                              const size_type j,
                              const number    value)
{
  (*this)(i, j) = value;
}
 
 
 
template <typename number>
inline typename LAPACKFullMatrix<number>::size_type
LAPACKFullMatrix<number>::m() const
{
  return static_cast<size_type>(this->n_rows());
}
 
 
 
template <typename number>
inline typename LAPACKFullMatrix<number>::size_type
LAPACKFullMatrix<number>::n() const
{
  return static_cast<size_type>(this->n_cols());
}
 
 
 
template <typename number>
template <typename MatrixType>
inline void
LAPACKFullMatrix<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();
    }
 
  state = LAPACKSupport::matrix;
}
 
 
 
template <typename number>
template <typename MatrixType>
inline void
LAPACKFullMatrix<number>::fill(const MatrixType &M,
                               const size_type   dst_offset_i,
                               const size_type   dst_offset_j,
                               const size_type   src_offset_i,
                               const size_type   src_offset_j,
                               const number      factor,
                               const bool        transpose)
{
  // loop over the elements of the argument matrix row by row, as suggested
  // in the documentation of the sparse matrix iterator class
  for (size_type row = src_offset_i; 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)
        {
          const size_type i = transpose ? entry->column() : row;
          const size_type j = transpose ? row : entry->column();
 
          const size_type dst_i = dst_offset_i + i - src_offset_i;
          const size_type dst_j = dst_offset_j + j - src_offset_j;
          if (dst_i < this->n_rows() && dst_j < this->n_cols())
            (*this)(dst_i, dst_j) = factor * entry->value();
        }
    }
 
  state = LAPACKSupport::matrix;
}
 
 
 
template <typename number>
template <typename number2>
void
LAPACKFullMatrix<number>::vmult(Vector<number2> &,
                                const Vector<number2> &,
                                const bool) const
{
  Assert(false,
         ExcMessage("LAPACKFullMatrix<number>::vmult must be called with a "
                    "matching Vector<double> vector type."));
}
 
 
 
template <typename number>
template <typename number2>
void
LAPACKFullMatrix<number>::vmult_add(Vector<number2> &,
                                    const Vector<number2> &) const
{
  Assert(false,
         ExcMessage("LAPACKFullMatrix<number>::vmult_add must be called with a "
                    "matching Vector<double> vector type."));
}
 
 
 
template <typename number>
template <typename number2>
void
LAPACKFullMatrix<number>::Tvmult(Vector<number2> &,
                                 const Vector<number2> &,
                                 const bool) const
{
  Assert(false,
         ExcMessage("LAPACKFullMatrix<number>::Tvmult must be called with a "
                    "matching Vector<double> vector type."));
}
 
 
 
template <typename number>
template <typename number2>
void
LAPACKFullMatrix<number>::Tvmult_add(Vector<number2> &,
                                     const Vector<number2> &) const
{
  Assert(false,
         ExcMessage("LAPACKFullMatrix<number>::Tvmult_add must be called "
                    "with a matching Vector<double> vector type."));
}
 
 
 
namespace internal
{
  namespace LAPACKFullMatrixImplementation
  {
    template <typename RealNumber>
    std::complex<RealNumber>
    pack_complex(const RealNumber &real_part, const RealNumber &imaginary_part)
    {
      return std::complex<RealNumber>(real_part, imaginary_part);
    }
 
    // The eigenvalues in LAPACKFullMatrix with complex-valued matrices are
    // contained in the 'wi' array, ignoring the 'wr' array.
    template <typename Number>
    std::complex<Number>
    pack_complex(const Number &, const std::complex<Number> &complex_number)
    {
      return complex_number;
    }
  } // namespace LAPACKFullMatrixImplementation
} // namespace internal
 
 
 
template <typename number>
inline std::complex<typename numbers::NumberTraits<number>::real_type>
LAPACKFullMatrix<number>::eigenvalue(const size_type i) const
{
  Assert(state & LAPACKSupport::eigenvalues, ExcInvalidState());
  Assert(wr.size() == this->n_rows(), ExcInternalError());
  Assert(wi.size() == this->n_rows(), ExcInternalError());
  AssertIndexRange(i, this->n_rows());
 
  return internal::LAPACKFullMatrixImplementation::pack_complex(wr[i], wi[i]);
}
 
 
 
template <typename number>
inline number
LAPACKFullMatrix<number>::singular_value(const size_type i) const
{
  Assert(state == LAPACKSupport::svd || state == LAPACKSupport::inverse_svd,
         LAPACKSupport::ExcState(state));
  AssertIndexRange(i, wr.size());
 
  return wr[i];
}
 
 
 
template <typename number>
inline const LAPACKFullMatrix<number> &
LAPACKFullMatrix<number>::get_svd_u() const
{
  Assert(state == LAPACKSupport::svd || state == LAPACKSupport::inverse_svd,
         LAPACKSupport::ExcState(state));
 
  return *svd_u;
}
 
 
 
template <typename number>
inline const LAPACKFullMatrix<number> &
LAPACKFullMatrix<number>::get_svd_vt() const
{
  Assert(state == LAPACKSupport::svd || state == LAPACKSupport::inverse_svd,
         LAPACKSupport::ExcState(state));
 
  return *svd_vt;
}
 
 
 
DEAL_II_NAMESPACE_CLOSE
 
#endif