dense_matrix.h

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
 
// This library is free software; you can 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.
 
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
 
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
 
 
#ifndef LIBMESH_DENSE_MATRIX_H
#define LIBMESH_DENSE_MATRIX_H
 
// Local Includes
#include "libmesh/libmesh_common.h"
#include "libmesh/dense_matrix_base.h"
#include "libmesh/int_range.h"
 
// For the definition of PetscBLASInt.
#if (LIBMESH_HAVE_PETSC)
# include "libmesh/petsc_macro.h"
# ifdef I
#  define LIBMESH_SAW_I
# endif
 
#include "libmesh/ignore_warnings.h"
# include <petscsys.h>
#include "libmesh/restore_warnings.h"
 
# ifndef LIBMESH_SAW_I
#  undef I // Avoid complex.h contamination
# endif
#endif
 
// C++ includes
#include <vector>
#include <algorithm>
#include <initializer_list>
 
#ifdef LIBMESH_HAVE_METAPHYSICL
#include "metaphysicl/dualnumber_decl.h"
#include "metaphysicl/raw_type.h"
#endif
 
namespace libMesh
{
 
// Forward Declarations
template <typename T> class DenseVector;
 
template<typename T>
class DenseMatrix : public DenseMatrixBase<T>
{
public:
 
  DenseMatrix(const unsigned int new_m=0,
              const unsigned int new_n=0);
 
  template <typename T2>
  DenseMatrix(unsigned int nrow,
              unsigned int ncol,
              std::initializer_list<T2> init_list);
 
  DenseMatrix (DenseMatrix &&) = default;
  DenseMatrix (const DenseMatrix &) = default;
  DenseMatrix & operator= (const DenseMatrix &) = default;
  DenseMatrix & operator= (DenseMatrix &&) = default;
  virtual ~DenseMatrix() = default;
 
  virtual void zero() override final;
 
  DenseMatrix sub_matrix(unsigned int row_id, unsigned int row_size,
                         unsigned int col_id, unsigned int col_size) const;
 
  T operator() (const unsigned int i,
                const unsigned int j) const;
 
  T & operator() (const unsigned int i,
                  const unsigned int j);
 
  virtual T el(const unsigned int i,
               const unsigned int j) const override final
  { return (*this)(i,j); }
 
  virtual T & el(const unsigned int i,
                 const unsigned int j) override final
  { return (*this)(i,j); }
 
  virtual void left_multiply (const DenseMatrixBase<T> & M2) override final;
 
  template <typename T2>
  void left_multiply (const DenseMatrixBase<T2> & M2);
 
  virtual void right_multiply (const DenseMatrixBase<T> & M2) override final;
 
  template <typename T2>
  void right_multiply (const DenseMatrixBase<T2> & M2);
 
  void vector_mult (DenseVector<T> & dest,
                    const DenseVector<T> & arg) const;
 
  template <typename T2>
  void vector_mult (DenseVector<typename CompareTypes<T,T2>::supertype> & dest,
                    const DenseVector<T2> & arg) const;
 
  void vector_mult_transpose (DenseVector<T> & dest,
                              const DenseVector<T> & arg) const;
 
  template <typename T2>
  void vector_mult_transpose (DenseVector<typename CompareTypes<T,T2>::supertype> & dest,
                              const DenseVector<T2> & arg) const;
 
  void vector_mult_add (DenseVector<T> & dest,
                        const T factor,
                        const DenseVector<T> & arg) const;
 
  template <typename T2, typename T3>
  void vector_mult_add (DenseVector<typename CompareTypes<T, typename CompareTypes<T2,T3>::supertype>::supertype> & dest,
                        const T2 factor,
                        const DenseVector<T3> & arg) const;
 
  void get_principal_submatrix (unsigned int sub_m, unsigned int sub_n, DenseMatrix<T> & dest) const;
 
  void get_principal_submatrix (unsigned int sub_m, DenseMatrix<T> & dest) const;
 
  void outer_product(const DenseVector<T> & a, const DenseVector<T> & b);
 
  template <typename T2>
  DenseMatrix<T> & operator = (const DenseMatrix<T2> & other_matrix);
 
  void swap(DenseMatrix<T> & other_matrix);
 
  void resize(const unsigned int new_m,
              const unsigned int new_n);
 
  void scale (const T factor);
 
  void scale_column (const unsigned int col, const T factor);
 
  DenseMatrix<T> & operator *= (const T factor);
 
  template<typename T2, typename T3>
  typename std::enable_if<
    ScalarTraits<T2>::value, void >::type add (const T2 factor,
                                               const DenseMatrix<T3> & mat);
 
  bool operator== (const DenseMatrix<T> & mat) const;
 
  bool operator!= (const DenseMatrix<T> & mat) const;
 
  DenseMatrix<T> & operator+= (const DenseMatrix<T> & mat);
 
  DenseMatrix<T> & operator-= (const DenseMatrix<T> & mat);
 
  auto min () const -> decltype(libmesh_real(T(0)));
 
  auto max () const -> decltype(libmesh_real(T(0)));
 
  auto l1_norm () const;
 
  auto linfty_norm () const;
 
  void left_multiply_transpose (const DenseMatrix<T> & A);
 
  template <typename T2>
  void left_multiply_transpose (const DenseMatrix<T2> & A);
 
 
  void right_multiply_transpose (const DenseMatrix<T> & A);
 
  template <typename T2>
  void right_multiply_transpose (const DenseMatrix<T2> & A);
 
  T transpose (const unsigned int i,
               const unsigned int j) const;
 
  void get_transpose(DenseMatrix<T> & dest) const;
 
  std::vector<T> & get_values() { return _val; }
 
  const std::vector<T> & get_values() const { return _val; }
 
  void condense(const unsigned int i,
                const unsigned int j,
                const T val,
                DenseVector<T> & rhs)
  { DenseMatrixBase<T>::condense (i, j, val, rhs); }
 
  void lu_solve (const DenseVector<T> & b,
                 DenseVector<T> & x);
 
  template <typename T2>
  void cholesky_solve(const DenseVector<T2> & b,
                      DenseVector<T2> & x);
 
  void svd(DenseVector<Real> & sigma);
 
  void svd(DenseVector<Real> & sigma,
           DenseMatrix<Number> & U,
           DenseMatrix<Number> & VT);
 
  void svd_solve(const DenseVector<T> & rhs,
                 DenseVector<T> & x,
                 Real rcond=std::numeric_limits<Real>::epsilon()) const;
 
  void evd(DenseVector<T> & lambda_real,
           DenseVector<T> & lambda_imag);
 
  void evd_left(DenseVector<T> & lambda_real,
                DenseVector<T> & lambda_imag,
                DenseMatrix<T> & VL);
 
  void evd_right(DenseVector<T> & lambda_real,
                 DenseVector<T> & lambda_imag,
                 DenseMatrix<T> & VR);
 
  void evd_left_and_right(DenseVector<T> & lambda_real,
                          DenseVector<T> & lambda_imag,
                          DenseMatrix<T> & VL,
                          DenseMatrix<T> & VR);
 
  T det();
 
  // void inverse();
 
  bool use_blas_lapack;
 
  struct UseBlasLapack
  {
    static const bool value = false;
  };
 
private:
 
  std::vector<T> _val;
 
  void _lu_decompose ();
 
  void _lu_back_substitute (const DenseVector<T> & b,
                            DenseVector<T> & x) const;
 
  void _cholesky_decompose();
 
  template <typename T2>
  void _cholesky_back_substitute(const DenseVector<T2> & b,
                                 DenseVector<T2> & x) const;
 
  enum DecompositionType {LU=0, CHOLESKY=1, LU_BLAS_LAPACK, NONE};
 
  DecompositionType _decomposition_type;
 
  enum _BLAS_Multiply_Flag {
    LEFT_MULTIPLY = 0,
    RIGHT_MULTIPLY,
    LEFT_MULTIPLY_TRANSPOSE,
    RIGHT_MULTIPLY_TRANSPOSE
  };
 
  void _multiply_blas(const DenseMatrixBase<T> & other,
                      _BLAS_Multiply_Flag flag);
 
  void _lu_decompose_lapack();
 
  void _svd_lapack(DenseVector<Real> & sigma);
 
  void _svd_lapack(DenseVector<Real> & sigma,
                   DenseMatrix<Number> & U,
                   DenseMatrix<Number> & VT);
 
  void _svd_solve_lapack(const DenseVector<T> & rhs,
                         DenseVector<T> & x,
                         Real rcond) const;
 
  void _svd_helper (char JOBU,
                    char JOBVT,
                    std::vector<Real> & sigma_val,
                    std::vector<Number> & U_val,
                    std::vector<Number> & VT_val);
 
  void _evd_lapack(DenseVector<T> & lambda_real,
                   DenseVector<T> & lambda_imag,
                   DenseMatrix<T> * VL = nullptr,
                   DenseMatrix<T> * VR = nullptr);
 
#if (LIBMESH_HAVE_PETSC && LIBMESH_USE_REAL_NUMBERS)
  typedef PetscBLASInt pivot_index_t;
#else
  typedef int pivot_index_t;
#endif
  std::vector<pivot_index_t> _pivots;
 
  void _lu_back_substitute_lapack (const DenseVector<T> & b,
                                   DenseVector<T> & x);
 
  void _matvec_blas(T alpha, T beta,
                    DenseVector<T> & dest,
                    const DenseVector<T> & arg,
                    bool trans=false) const;
 
  template <typename T2>
  void _left_multiply_transpose (const DenseMatrix<T2> & A);
 
  template <typename T2>
  void _right_multiply_transpose (const DenseMatrix<T2> & A);
};
 
 
 
 
 
// ------------------------------------------------------------
namespace DenseMatrices
{
 
typedef DenseMatrix<Real> RealDenseMatrix;
 
typedef DenseMatrix<Complex> ComplexDenseMatrix;
 
}
 
 
 
using namespace DenseMatrices;
 
// The PETSc Lapack wrappers are only for PetscScalar, therefore we
// can't e.g. get a Lapack version of DenseMatrix<Real>::lu_solve()
// when libmesh/PETSc are compiled with complex numbers.
#if defined(LIBMESH_HAVE_PETSC) && \
  defined(LIBMESH_USE_REAL_NUMBERS) && \
  defined(LIBMESH_DEFAULT_DOUBLE_PRECISION)
template <>
struct DenseMatrix<double>::UseBlasLapack
{
  static const bool value = true;
};
#endif
 
 
// ------------------------------------------------------------
// Dense Matrix member functions
template<typename T>
inline
DenseMatrix<T>::DenseMatrix(const unsigned int new_m,
                            const unsigned int new_n) :
  DenseMatrixBase<T>(new_m,new_n),
  use_blas_lapack(DenseMatrix<T>::UseBlasLapack::value),
  _val(),
  _decomposition_type(NONE)
{
  this->resize(new_m,new_n);
}
 
template <typename T>
template <typename T2>
DenseMatrix<T>::DenseMatrix(unsigned int nrow,
                            unsigned int ncol,
                            std::initializer_list<T2> init_list) :
  DenseMatrixBase<T>(nrow, ncol),
  use_blas_lapack(DenseMatrix<T>::UseBlasLapack::value),
  _val(init_list.begin(), init_list.end()),
  _decomposition_type(NONE)
{
  // Make sure the user passed us an amount of data which is
  // consistent with the size of the matrix.
  libmesh_assert_equal_to(nrow * ncol, init_list.size());
}
 
 
 
template<typename T>
inline
void DenseMatrix<T>::swap(DenseMatrix<T> & other_matrix)
{
  std::swap(this->_m, other_matrix._m);
  std::swap(this->_n, other_matrix._n);
  _val.swap(other_matrix._val);
  DecompositionType _temp = _decomposition_type;
  _decomposition_type = other_matrix._decomposition_type;
  other_matrix._decomposition_type = _temp;
}
 
 
template <typename T>
template <typename T2>
inline
DenseMatrix<T> &
DenseMatrix<T>::operator=(const DenseMatrix<T2> & mat)
{
  unsigned int mat_m = mat.m(), mat_n = mat.n();
  this->resize(mat_m, mat_n);
  for (unsigned int i=0; i<mat_m; i++)
    for (unsigned int j=0; j<mat_n; j++)
      (*this)(i,j) = mat(i,j);
 
  return *this;
}
 
 
 
template<typename T>
inline
void DenseMatrix<T>::resize(const unsigned int new_m,
                            const unsigned int new_n)
{
  _val.resize(new_m*new_n);
 
  this->_m = new_m;
  this->_n = new_n;
 
  // zero and set decomposition_type to NONE
  this->zero();
}
 
 
 
template<typename T>
inline
void DenseMatrix<T>::zero()
{
  _decomposition_type = NONE;
 
  std::fill (_val.begin(), _val.end(), static_cast<T>(0));
}
 
 
 
template<typename T>
inline
DenseMatrix<T> DenseMatrix<T>::sub_matrix(unsigned int row_id, unsigned int row_size,
                                          unsigned int col_id, unsigned int col_size) const
{
  libmesh_assert_less (row_id + row_size - 1, this->_m);
  libmesh_assert_less (col_id + col_size - 1, this->_n);
 
  DenseMatrix<T> sub;
  sub._m = row_size;
  sub._n = col_size;
  sub._val.resize(row_size * col_size);
 
  unsigned int end_col = this->_n - col_size - col_id;
  unsigned int p = row_id * this->_n;
  unsigned int q = 0;
  for (unsigned int i=0; i<row_size; i++)
  {
    // skip the beginning columns
    p += col_id;
    for (unsigned int j=0; j<col_size; j++)
      sub._val[q++] = _val[p++];
    // skip the rest columns
    p += end_col;
  }
 
  return sub;
}
 
 
 
template<typename T>
inline
T DenseMatrix<T>::operator () (const unsigned int i,
                               const unsigned int j) const
{
  libmesh_assert_less (i*j, _val.size());
  libmesh_assert_less (i, this->_m);
  libmesh_assert_less (j, this->_n);
 
 
  //  return _val[(i) + (this->_m)*(j)]; // col-major
  return _val[(i)*(this->_n) + (j)]; // row-major
}
 
 
 
template<typename T>
inline
T & DenseMatrix<T>::operator () (const unsigned int i,
                                 const unsigned int j)
{
  libmesh_assert_less (i*j, _val.size());
  libmesh_assert_less (i, this->_m);
  libmesh_assert_less (j, this->_n);
 
  //return _val[(i) + (this->_m)*(j)]; // col-major
  return _val[(i)*(this->_n) + (j)]; // row-major
}
 
 
 
 
 
template<typename T>
inline
void DenseMatrix<T>::scale (const T factor)
{
  for (auto & v : _val)
    v *= factor;
}
 
 
template<typename T>
inline
void DenseMatrix<T>::scale_column (const unsigned int col, const T factor)
{
  for (auto i : make_range(this->m()))
    (*this)(i, col) *= factor;
}
 
 
 
template<typename T>
inline
DenseMatrix<T> & DenseMatrix<T>::operator *= (const T factor)
{
  this->scale(factor);
  return *this;
}
 
 
 
template<typename T>
template<typename T2, typename T3>
inline
typename std::enable_if<
  ScalarTraits<T2>::value, void >::type
DenseMatrix<T>::add (const T2 factor,
                     const DenseMatrix<T3> & mat)
{
  libmesh_assert_equal_to (this->m(), mat.m());
  libmesh_assert_equal_to (this->n(), mat.n());
 
  for (auto i : make_range(this->m()))
    for (auto j : make_range(this->n()))
      (*this)(i,j) += factor * mat(i,j);
}
 
 
 
template<typename T>
inline
bool DenseMatrix<T>::operator == (const DenseMatrix<T> & mat) const
{
  for (auto i : index_range(_val))
    if (_val[i] != mat._val[i])
      return false;
 
  return true;
}
 
 
 
template<typename T>
inline
bool DenseMatrix<T>::operator != (const DenseMatrix<T> & mat) const
{
  for (auto i : index_range(_val))
    if (_val[i] != mat._val[i])
      return true;
 
  return false;
}
 
 
 
template<typename T>
inline
DenseMatrix<T> & DenseMatrix<T>::operator += (const DenseMatrix<T> & mat)
{
  for (auto i : index_range(_val))
    _val[i] += mat._val[i];
 
  return *this;
}
 
 
 
template<typename T>
inline
DenseMatrix<T> & DenseMatrix<T>::operator -= (const DenseMatrix<T> & mat)
{
  for (auto i : index_range(_val))
    _val[i] -= mat._val[i];
 
  return *this;
}
 
 
 
template<typename T>
inline
auto DenseMatrix<T>::min () const -> decltype(libmesh_real(T(0)))
{
  libmesh_assert (this->_m);
  libmesh_assert (this->_n);
  auto my_min = libmesh_real((*this)(0,0));
 
  for (unsigned int i=0; i!=this->_m; i++)
    {
      for (unsigned int j=0; j!=this->_n; j++)
        {
          auto current = libmesh_real((*this)(i,j));
          my_min = (my_min < current? my_min : current);
        }
    }
  return my_min;
}
 
 
 
template<typename T>
inline
auto DenseMatrix<T>::max () const -> decltype(libmesh_real(T(0)))
{
  libmesh_assert (this->_m);
  libmesh_assert (this->_n);
  auto my_max = libmesh_real((*this)(0,0));
 
  for (unsigned int i=0; i!=this->_m; i++)
    {
      for (unsigned int j=0; j!=this->_n; j++)
        {
          auto current = libmesh_real((*this)(i,j));
          my_max = (my_max > current? my_max : current);
        }
    }
  return my_max;
}
 
 
 
template<typename T>
inline
auto DenseMatrix<T>::l1_norm () const
{
  libmesh_assert (this->_m);
  libmesh_assert (this->_n);
 
  using std::abs;
  auto columnsum = abs(T(0));
  for (unsigned int i=0; i!=this->_m; i++)
    {
      columnsum += abs((*this)(i,0));
    }
  auto my_max = columnsum;
  for (unsigned int j=1; j!=this->_n; j++)
    {
      columnsum = 0.;
      for (unsigned int i=0; i!=this->_m; i++)
        {
          columnsum += abs((*this)(i,j));
        }
      my_max = (my_max > columnsum? my_max : columnsum);
    }
  return my_max;
}
 
 
 
template<typename T>
inline
auto DenseMatrix<T>::linfty_norm () const
{
  libmesh_assert (this->_m);
  libmesh_assert (this->_n);
  using std::abs;
 
  auto rowsum = abs(T(0));
  for (unsigned int j=0; j!=this->_n; j++)
    {
      rowsum += abs((*this)(0,j));
    }
  auto my_max = rowsum;
  for (unsigned int i=1; i!=this->_m; i++)
    {
      rowsum = 0.;
      for (unsigned int j=0; j!=this->_n; j++)
        {
          rowsum += abs((*this)(i,j));
        }
      my_max = (my_max > rowsum? my_max : rowsum);
    }
  return my_max;
}
 
 
 
template<typename T>
inline
T DenseMatrix<T>::transpose (const unsigned int i,
                             const unsigned int j) const
{
  // Implement in terms of operator()
  return (*this)(j,i);
}
 
 
 
 
 
// template<typename T>
// inline
// void DenseMatrix<T>::condense(const unsigned int iv,
//       const unsigned int jv,
//       const T val,
//       DenseVector<T> & rhs)
// {
//   libmesh_assert_equal_to (this->_m, rhs.size());
//   libmesh_assert_equal_to (iv, jv);
 
 
//   // move the known value into the RHS
//   // and zero the column
//   for (auto i : make_range(this->m()))
//     {
//       rhs(i) -= ((*this)(i,jv))*val;
//       (*this)(i,jv) = 0.;
//     }
 
//   // zero the row
//   for (auto j : make_range(this->n()))
//     (*this)(iv,j) = 0.;
 
//   (*this)(iv,jv) = 1.;
//   rhs(iv) = val;
 
// }
 
 
// A matrix is finite iff every component is
template <typename T>
bool isfinite (const DenseMatrix<T> & var)
{
  using std::isfinite;
  using libMesh::isfinite; // for T==complex
  const auto m = var.m(), n = var.n();
  for (auto i : make_range(m))
    for (auto j : make_range(n))
      if (!isfinite(var(i,j)))
        return false;
  return true;
}
 
 
// A matrix is infinite iff some component is infinite but no
// component is NaN.
//
// This is arguably inconsistent with our std::complex overload (and
// the C99 Annex G recommendations for _Complex, and C++ std::complex
// arithmetic), which treats mixed (inf,NaN) pairs as infinite, but
// this is probably safer for users.
template <typename T>
bool isinf (const DenseMatrix<T> & var)
{
  using std::isinf;
  using libMesh::isinf; // for T==complex
  using std::isnan;
  using libMesh::isnan;
  const auto m = var.m(), n = var.n();
  bool has_inf = false;
  for (auto i : make_range(m))
    for (auto j : make_range(n))
      {
        // NaN anywhere makes us NaN, not inf
        if (isnan(var(i,j)))
          return false;
        has_inf = has_inf || isinf(var(i,j));
      }
  return has_inf;
}
 
 
 
// A matrix is NaN iff some component is NaN
//
// This is arguably inconsistent with our std::complex overload (and
// the C99 Annex G recommendations for _Complex, and C++ std::complex
// arithmetic), which treats mixed (inf,NaN) pairs as infinite, but
// this is probably safer for users.
template <typename T>
bool isnan (const DenseMatrix<T> & var)
{
  using std::isnan;
  using libMesh::isnan; // for T==complex
  const auto m = var.m(), n = var.n();
  for (auto i : make_range(m))
    for (auto j : make_range(n))
      if (isnan(var(i,j)))
        return true;
  return false;
}
 
 
} // namespace libMesh
 
#ifdef LIBMESH_HAVE_METAPHYSICL
namespace MetaPhysicL
{
template <typename T>
struct RawType<libMesh::DenseMatrix<T>>
{
  typedef libMesh::DenseMatrix<typename RawType<T>::value_type> value_type;
 
  static value_type value (const libMesh::DenseMatrix<T> & in)
    {
      const auto m = in.m(), n = in.n();
      value_type ret(m, n);
      for (unsigned int i = 0; i < m; ++i)
        for (unsigned int j = 0; j < n; ++j)
          ret(i,j) = raw_value(in(i,j));
 
      return ret;
    }
};
}
#endif
 
 
#endif // LIBMESH_DENSE_MATRIX_H