dense_matrix.h

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2024 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 boostcopy::enable_if_c<
    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 -> decltype(std::abs(T(0)));

  auto linfty_norm () const -> decltype(std::abs(T(0)));

  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 boostcopy::enable_if_c<
  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 -> decltype(std::abs(T(0)))
{
  libmesh_assert (this->_m);
  libmesh_assert (this->_n);

  auto columnsum = std::abs(T(0));
  for (unsigned int i=0; i!=this->_m; i++)
    {
      columnsum += std::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 += std::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 -> decltype(std::abs(T(0)))
{
  libmesh_assert (this->_m);
  libmesh_assert (this->_n);

  auto rowsum = std::abs(T(0));
  for (unsigned int j=0; j!=this->_n; j++)
    {
      rowsum += std::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 += std::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;

// }




} // 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