libMesh::DenseMatrix< T > Class Template Reference

Defines a dense matrix for use in Finite Element-type computations. More...

#include <dense_matrix.h>

Inheritance diagram for libMesh::DenseMatrix< T >:

Classes
struct	UseBlasLapack
	Helper structure for determining whether to use blas_lapack. More...

Public Member Functions
	DenseMatrix (const unsigned int new_m=0, const unsigned int new_n=0)
	Constructor.
template<typename T2 >
	DenseMatrix (unsigned int nrow, unsigned int ncol, std::initializer_list< T2 > init_list)
	Constructor taking the number of rows, columns, and an initializer_list, which must be of length nrow * ncol, of row-major values to initialize the DenseMatrix with.
	DenseMatrix (DenseMatrix &&)=default
	The 5 special functions can be defaulted for this class, as it does not manage any memory itself.
	DenseMatrix (const DenseMatrix &)=default
DenseMatrix &	operator= (const DenseMatrix &)=default
DenseMatrix &	operator= (DenseMatrix &&)=default
virtual	~DenseMatrix ()=default
virtual void	zero () override final
	Sets all elements of the matrix to 0 and resets any decomposition flag which may have been previously set.
DenseMatrix	sub_matrix (unsigned int row_id, unsigned int row_size, unsigned int col_id, unsigned int col_size) const
	Get submatrix with the smallest row and column indices and the submatrix size.
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
virtual T &	el (const unsigned int i, const unsigned int j) override final
virtual void	left_multiply (const DenseMatrixBase< T > &M2) override final
	Performs the operation: (*this) <- M2 * (*this)
template<typename T2 >
void	left_multiply (const DenseMatrixBase< T2 > &M2)
	Left multiplies by the matrix M2 of different type.
virtual void	right_multiply (const DenseMatrixBase< T > &M2) override final
	Performs the operation: (*this) <- (*this) * M3.
template<typename T2 >
void	right_multiply (const DenseMatrixBase< T2 > &M2)
	Right multiplies by the matrix M2 of different type.
void	vector_mult (DenseVector< T > &dest, const DenseVector< T > &arg) const
	Performs the matrix-vector multiplication, dest := (*this) * arg.
template<typename T2 >
void	vector_mult (DenseVector< typename CompareTypes< T, T2 >::supertype > &dest, const DenseVector< T2 > &arg) const
	Performs the matrix-vector multiplication, dest := (*this) * arg on mixed types.
void	vector_mult_transpose (DenseVector< T > &dest, const DenseVector< T > &arg) const
	Performs the matrix-vector multiplication, dest := (*this)^T * arg.
template<typename T2 >
void	vector_mult_transpose (DenseVector< typename CompareTypes< T, T2 >::supertype > &dest, const DenseVector< T2 > &arg) const
	Performs the matrix-vector multiplication, dest := (*this)^T * arg.
void	vector_mult_add (DenseVector< T > &dest, const T factor, const DenseVector< T > &arg) const
	Performs the scaled matrix-vector multiplication, dest += factor * (*this) * arg.
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
	Performs the scaled matrix-vector multiplication, dest += factor * (*this) * arg.
void	get_principal_submatrix (unsigned int sub_m, unsigned int sub_n, DenseMatrix< T > &dest) const
	Put the sub_m x sub_n principal submatrix into dest.
void	get_principal_submatrix (unsigned int sub_m, DenseMatrix< T > &dest) const
	Put the sub_m x sub_m principal submatrix into dest.
void	outer_product (const DenseVector< T > &a, const DenseVector< T > &b)
	Computes the outer (dyadic) product of two vectors and stores in (*this).
template<typename T2 >
DenseMatrix< T > &	operator= (const DenseMatrix< T2 > &other_matrix)
	Assignment-from-other-matrix-type operator.
void	swap (DenseMatrix< T > &other_matrix)
	STL-like swap method.
void	resize (const unsigned int new_m, const unsigned int new_n)
	Resizes the matrix to the specified size and calls zero().
void	scale (const T factor)
	Multiplies every element in the matrix by factor.
void	scale_column (const unsigned int col, const T factor)
	Multiplies every element in the column col matrix by factor.
DenseMatrix< T > &	operator*= (const T factor)
	Multiplies every element in the matrix by factor.
template<typename T2 , typename T3 >
std::enable_if< ScalarTraits< T2 >::value, void >::type	add (const T2 factor, const DenseMatrix< T3 > &mat)
	Adds factor times mat to this matrix.
bool	operator== (const DenseMatrix< T > &mat) const
bool	operator!= (const DenseMatrix< T > &mat) const
DenseMatrix< T > &	operator+= (const DenseMatrix< T > &mat)
	Adds mat to this matrix.
DenseMatrix< T > &	operator-= (const DenseMatrix< T > &mat)
	Subtracts mat from this matrix.
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)
	Left multiplies by the transpose of the matrix A.
template<typename T2 >
void	left_multiply_transpose (const DenseMatrix< T2 > &A)
	Left multiplies by the transpose of the matrix A which contains a different numerical type.
void	right_multiply_transpose (const DenseMatrix< T > &A)
	Right multiplies by the transpose of the matrix A.
template<typename T2 >
void	right_multiply_transpose (const DenseMatrix< T2 > &A)
	Right multiplies by the transpose of the matrix A which contains a different numerical type.
T	transpose (const unsigned int i, const unsigned int j) const
void	get_transpose (DenseMatrix< T > &dest) const
	Put the tranposed matrix into dest.
std::vector< T > &	get_values ()
const std::vector< T > &	get_values () const
void	condense (const unsigned int i, const unsigned int j, const T val, DenseVector< T > &rhs)
	Condense-out the (i,j) entry of the matrix, forcing it to take on the value val.
void	lu_solve (const DenseVector< T > &b, DenseVector< T > &x)
	Solve the system Ax=b given the input vector b.
template<typename T2 >
void	cholesky_solve (const DenseVector< T2 > &b, DenseVector< T2 > &x)
	For symmetric positive definite (SPD) matrices.
void	svd (DenseVector< Real > &sigma)
	Compute the singular value decomposition of the matrix.
void	svd (DenseVector< Real > &sigma, DenseMatrix< Number > &U, DenseMatrix< Number > &VT)
	Compute the "reduced" singular value decomposition of the matrix.
void	svd_solve (const DenseVector< T > &rhs, DenseVector< T > &x, Real rcond=std::numeric_limits< Real >::epsilon()) const
	Solve the system of equations \( A x = rhs \) for \( x \) in the least-squares sense.
void	evd (DenseVector< T > &lambda_real, DenseVector< T > &lambda_imag)
	Compute the eigenvalues (both real and imaginary parts) of a general matrix.
void	evd_left (DenseVector< T > &lambda_real, DenseVector< T > &lambda_imag, DenseMatrix< T > &VL)
	Compute the eigenvalues (both real and imaginary parts) and left eigenvectors of a general matrix, \( A \).
void	evd_right (DenseVector< T > &lambda_real, DenseVector< T > &lambda_imag, DenseMatrix< T > &VR)
	Compute the eigenvalues (both real and imaginary parts) and right eigenvectors of a general matrix, \( A \).
void	evd_left_and_right (DenseVector< T > &lambda_real, DenseVector< T > &lambda_imag, DenseMatrix< T > &VL, DenseMatrix< T > &VR)
	Compute the eigenvalues (both real and imaginary parts) as well as the left and right eigenvectors of a general matrix.
T	det ()
unsigned int	m () const
unsigned int	n () const
void	print (std::ostream &os=libMesh::out) const
	Pretty-print the matrix, by default to libMesh::out.
void	print_scientific (std::ostream &os, unsigned precision=8) const
	Prints the matrix entries with more decimal places in scientific notation.
template<typename T2 , typename T3 >
std::enable_if< ScalarTraits< T2 >::value, void >::type	add (const T2 factor, const DenseMatrixBase< T3 > &mat)
	Adds factor to every element in the matrix.
DenseVector< T >	diagonal () const
	Return the matrix diagonal.

Public Attributes
bool	use_blas_lapack
	Computes the inverse of the dense matrix (assuming it is invertible) by first computing the LU decomposition and then performing multiple back substitution steps.

Protected Member Functions
void	condense (const unsigned int i, const unsigned int j, const T val, DenseVectorBase< T > &rhs)
	Condense-out the (i,j) entry of the matrix, forcing it to take on the value val.

Static Protected Member Functions
static void	multiply (DenseMatrixBase< T > &M1, const DenseMatrixBase< T > &M2, const DenseMatrixBase< T > &M3)
	Helper function - Performs the computation M1 = M2 * M3 where: M1 = (m x n) M2 = (m x p) M3 = (p x n)

Protected Attributes
unsigned int	_m
	The row dimension.
unsigned int	_n
	The column dimension.

Private Types
enum	DecompositionType { LU =0 , CHOLESKY =1 , LU_BLAS_LAPACK , NONE }
	The decomposition schemes above change the entries of the matrix A. More...
enum	_BLAS_Multiply_Flag { LEFT_MULTIPLY = 0 , RIGHT_MULTIPLY , LEFT_MULTIPLY_TRANSPOSE , RIGHT_MULTIPLY_TRANSPOSE }
	Enumeration used to determine the behavior of the _multiply_blas function. More...
typedef PetscBLASInt	pivot_index_t
	Array used to store pivot indices.
typedef int	pivot_index_t

Private Member Functions
void	_lu_decompose ()
	Form the LU decomposition of the matrix.
void	_lu_back_substitute (const DenseVector< T > &b, DenseVector< T > &x) const
	Solves the system Ax=b through back substitution.
void	_cholesky_decompose ()
	Decomposes a symmetric positive definite matrix into a product of two lower triangular matrices according to A = LL^T.
template<typename T2 >
void	_cholesky_back_substitute (const DenseVector< T2 > &b, DenseVector< T2 > &x) const
	Solves the equation Ax=b for the unknown value x and rhs b based on the Cholesky factorization of A.
void	_multiply_blas (const DenseMatrixBase< T > &other, _BLAS_Multiply_Flag flag)
	The _multiply_blas function computes A <- op(A) * op(B) using BLAS gemm function.
void	_lu_decompose_lapack ()
	Computes an LU factorization of the matrix using the Lapack routine "getrf".
void	_svd_lapack (DenseVector< Real > &sigma)
	Computes an SVD of the matrix using the Lapack routine "getsvd".
void	_svd_lapack (DenseVector< Real > &sigma, DenseMatrix< Number > &U, DenseMatrix< Number > &VT)
	Computes a "reduced" SVD of the matrix using the Lapack routine "getsvd".
void	_svd_solve_lapack (const DenseVector< T > &rhs, DenseVector< T > &x, Real rcond) const
	Called by svd_solve(rhs).
void	_svd_helper (char JOBU, char JOBVT, std::vector< Real > &sigma_val, std::vector< Number > &U_val, std::vector< Number > &VT_val)
	Helper function that actually performs the SVD.
void	_evd_lapack (DenseVector< T > &lambda_real, DenseVector< T > &lambda_imag, DenseMatrix< T > *VL=nullptr, DenseMatrix< T > *VR=nullptr)
	Computes the eigenvalues of the matrix using the Lapack routine "DGEEV".
void	_lu_back_substitute_lapack (const DenseVector< T > &b, DenseVector< T > &x)
	Companion function to _lu_decompose_lapack().
void	_matvec_blas (T alpha, T beta, DenseVector< T > &dest, const DenseVector< T > &arg, bool trans=false) const
	Uses the BLAS GEMV function (through PETSc) to compute.
template<typename T2 >
void	_left_multiply_transpose (const DenseMatrix< T2 > &A)
	Left multiplies by the transpose of the matrix A which may contain a different numerical type.
template<typename T2 >
void	_right_multiply_transpose (const DenseMatrix< T2 > &A)
	Right multiplies by the transpose of the matrix A which may contain a different numerical type.

Private Attributes
std::vector< T >	_val
	The actual data values, stored as a 1D array.
DecompositionType	_decomposition_type
	This flag keeps track of which type of decomposition has been performed on the matrix.
std::vector< pivot_index_t >	_pivots

Detailed Description

template<typename T>
class libMesh::DenseMatrix< T >

Defines a dense matrix for use in Finite Element-type computations.

Useful for storing element stiffness matrices before summation into a global matrix. All overridden virtual functions are documented in dense_matrix_base.h.

Author

Benjamin S. Kirk

Date

2002

A matrix object used for finite element assembly and numerics.

Definition at line 71 of file dense_matrix.h.

Member Typedef Documentation

pivot_index_t [1/2]

template<typename T >

typedef PetscBLASInt libMesh::DenseMatrix< T >::pivot_index_t

private

Array used to store pivot indices.

May be used by whatever factorization is currently active, clients of the class should not rely on it for any reason.

Definition at line 736 of file dense_matrix.h.

pivot_index_t [2/2]

template<typename T >

typedef int libMesh::DenseMatrix< T >::pivot_index_t

private

Definition at line 738 of file dense_matrix.h.

Member Enumeration Documentation

_BLAS_Multiply_Flag

template<typename T >

enum libMesh::DenseMatrix::_BLAS_Multiply_Flag

private

Enumeration used to determine the behavior of the _multiply_blas function.

Enumerator
LEFT_MULTIPLY
RIGHT_MULTIPLY
LEFT_MULTIPLY_TRANSPOSE
RIGHT_MULTIPLY_TRANSPOSE

Definition at line 656 of file dense_matrix.h.

                           {
    LEFT_MULTIPLY = 0,
    RIGHT_MULTIPLY,
    LEFT_MULTIPLY_TRANSPOSE,
    RIGHT_MULTIPLY_TRANSPOSE
  };

DecompositionType

template<typename T >

enum libMesh::DenseMatrix::DecompositionType

private

The decomposition schemes above change the entries of the matrix A.

It is therefore an error to call A.lu_solve() and subsequently call A.cholesky_solve() since the result will probably not match any desired outcome. This typedef keeps track of which decomposition has been called for this matrix.

Enumerator
LU
CHOLESKY
LU_BLAS_LAPACK
NONE

Definition at line 644 of file dense_matrix.h.

{LU=0, CHOLESKY=1, LU_BLAS_LAPACK, NONE};

Constructor & Destructor Documentation

DenseMatrix() [1/4]

template<typename T >

libMesh::DenseMatrix< T >::DenseMatrix ( const unsigned int  new_m = 0, const unsigned int  new_n = 0  )

inline

Constructor.

Creates a dense matrix of dimension m by n.

Definition at line 836 of file dense_matrix.h.

                                                      :
  DenseMatrixBase<T>(new_m,new_n),
  use_blas_lapack(DenseMatrix<T>::UseBlasLapack::value),
  _val(),
  _decomposition_type(NONE)
{
  this->resize(new_m,new_n);
}

References libMesh::DenseMatrix< T >::resize().

DenseMatrix() [2/4]

template<typename T >

template<typename T2 >

libMesh::DenseMatrix< T >::DenseMatrix	(	unsigned int	nrow,
		unsigned int	ncol,
		std::initializer_list< T2 >	init_list
	)

Constructor taking the number of rows, columns, and an initializer_list, which must be of length nrow * ncol, of row-major values to initialize the DenseMatrix with.

Definition at line 848 of file dense_matrix.h.

                                                               :
  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());
}

DenseMatrix() [3/4]

template<typename T >

libMesh::DenseMatrix< T >::DenseMatrix ( DenseMatrix< T > &&  )

default

The 5 special functions can be defaulted for this class, as it does not manage any memory itself.

DenseMatrix() [4/4]

template<typename T >

libMesh::DenseMatrix< T >::DenseMatrix ( const DenseMatrix< T > &  )

default

~DenseMatrix()

template<typename T >

virtual libMesh::DenseMatrix< T >::~DenseMatrix ( )

virtualdefault

Member Function Documentation

_cholesky_back_substitute()

template<typename T >

template<typename T2 >

template LIBMESH_EXPORT void libMesh::DenseMatrix< T >::_cholesky_back_substitute ( const DenseVector< T2 > &  b, DenseVector< T2 > &  x  ) const

private

Solves the equation Ax=b for the unknown value x and rhs b based on the Cholesky factorization of A.

Note

This method may be used when A is real-valued and b and x are complex-valued.

Definition at line 924 of file dense_matrix_impl.h.

{
  // Shorthand notation for number of rows and columns.
  const unsigned int
    n_rows = this->m(),
    n_cols = this->n();
 
  // Just to be really sure...
  libmesh_assert_equal_to (n_rows, n_cols);
 
  // A convenient reference to *this
  const DenseMatrix<T> & A = *this;
 
  // Now compute the solution to Ax =b using the factorization.
  x.resize(n_rows);
 
  // Solve for Ly=b
  for (unsigned int i=0; i<n_cols; ++i)
    {
      T2 temp = b(i);
 
      for (unsigned int k=0; k<i; ++k)
        temp -= A(i,k)*x(k);
 
      x(i) = temp / A(i,i);
    }
 
  // Solve for L^T x = y
  for (unsigned int i=0; i<n_cols; ++i)
    {
      const unsigned int ib = (n_cols-1)-i;
 
      for (unsigned int k=(ib+1); k<n_cols; ++k)
        x(ib) -= A(k,ib) * x(k);
 
      x(ib) /= A(ib,ib);
    }
}

References b, and libMesh::DenseVector< T >::resize().

_cholesky_decompose()

template<typename T >

void libMesh::DenseMatrix< T >::_cholesky_decompose ( )

private

Decomposes a symmetric positive definite matrix into a product of two lower triangular matrices according to A = LL^T.

Note

This program generates an error if the matrix is not SPD.

Definition at line 876 of file dense_matrix_impl.h.

{
  using std::sqrt;
 
  // If we called this function, there better not be any
  // previous decomposition of the matrix.
  libmesh_assert_equal_to (this->_decomposition_type, NONE);
 
  // Shorthand notation for number of rows and columns.
  const unsigned int
    n_rows = this->m(),
    n_cols = this->n();
 
  // Just to be really sure...
  libmesh_assert_equal_to (n_rows, n_cols);
 
  // A convenient reference to *this
  DenseMatrix<T> & A = *this;
 
  for (unsigned int i=0; i<n_rows; ++i)
    {
      for (unsigned int j=i; j<n_cols; ++j)
        {
          for (unsigned int k=0; k<i; ++k)
            A(i,j) -= A(i,k) * A(j,k);
 
          if (i == j)
            {
#ifndef LIBMESH_USE_COMPLEX_NUMBERS
              libmesh_error_msg_if(A(i,j) <= 0.0,
                                   "Error! Can only use Cholesky decomposition with symmetric positive definite matrices.");
#endif
 
              A(i,i) = sqrt(A(i,j));
            }
          else
            A(j,i) = A(i,j) / A(i,i);
        }
    }
 
  // Set the flag for CHOLESKY decomposition
  this->_decomposition_type = CHOLESKY;
}

_evd_lapack()

template<typename T >

template LIBMESH_EXPORT void libMesh::DenseMatrix< T >::_evd_lapack ( DenseVector< T > &  lambda_real, DenseVector< T > &  lambda_imag, DenseMatrix< T > *  VL = nullptr, DenseMatrix< T > *  VR = nullptr  )

private

Computes the eigenvalues of the matrix using the Lapack routine "DGEEV".

If VR and/or VL are not nullptr, then the matrix of right and/or left eigenvectors is also computed and returned by this function.

[ Implementation in dense_matrix_blas_lapack.C ]

Definition at line 695 of file dense_matrix_blas_lapack.C.

{
  // This algorithm only works for square matrices, so verify this up front.
  libmesh_error_msg_if(this->m() != this->n(), "Can only compute eigen-decompositions for square matrices.");
 
  // If the user requests left or right eigenvectors, we have to make
  // sure and pass the transpose of this matrix to Lapack, otherwise
  // it will compute the inverse transpose of what we are
  // after... since we know the matrix is square, we can just swap
  // entries in place.  If the user does not request eigenvectors, we
  // can skip this extra step, since the eigenvalues for A and A^T are
  // the same.
  if (VL || VR)
    {
      for (auto i : make_range(this->_m))
        for (auto j : make_range(i))
          std::swap((*this)(i,j), (*this)(j,i));
    }
 
  // The calling sequence for dgeev is:
  // DGEEV (character  JOBVL,
  //        character  JOBVR,
  //        integer N,
  //        double precision, dimension( lda, * )  A,
  //        integer  LDA,
  //        double precision, dimension( * )  WR,
  //        double precision, dimension( * )  WI,
  //        double precision, dimension( ldvl, * )  VL,
  //        integer  LDVL,
  //        double precision, dimension( ldvr, * )  VR,
  //        integer  LDVR,
  //        double precision, dimension( * )  WORK,
  //        integer  LWORK,
  //        integer  INFO)
 
  // JOBVL (input)
  //       = 'N': left eigenvectors of A are not computed;
  //       = 'V': left eigenvectors of A are computed.
  char JOBVL = VL ? 'V' : 'N';
 
  // JOBVR (input)
  //       = 'N': right eigenvectors of A are not computed;
  //       = 'V': right eigenvectors of A are computed.
  char JOBVR = VR ? 'V' : 'N';
 
  // N (input)
  //   The number of rows/cols of the matrix A.  N >= 0.
  PetscBLASInt N = this->m();
 
  // A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
  //   On entry, the N-by-N matrix A.
  //   On exit, A has been overwritten.
 
  // LDA (input)
  //     The leading dimension of the array A.  LDA >= max(1,N).
  PetscBLASInt LDA = N;
 
  // WR (output) double precision array, dimension (N)
  // WI (output) double precision array, dimension (N)
  //    WR and WI contain the real and imaginary parts,
  //    respectively, of the computed eigenvalues.  Complex
  //    conjugate pairs of eigenvalues appear consecutively
  //    with the eigenvalue having the positive imaginary part
  //    first.
  lambda_real.resize(N);
  lambda_imag.resize(N);
 
  // VL (output) double precision array, dimension (LDVL,N)
  //    If JOBVL = 'V', the left eigenvectors u(j) are stored one
  //    after another in the columns of VL, in the same order
  //    as their eigenvalues.
  //    If JOBVL = 'N', VL is not referenced.
  //    If the j-th eigenvalue is real, then u(j) = VL(:,j),
  //    the j-th column of VL.
  //    If the j-th and (j+1)-st eigenvalues form a complex
  //    conjugate pair, then u(j) = VL(:,j) + i*VL(:,j+1) and
  //    u(j+1) = VL(:,j) - i*VL(:,j+1).
  // Will be set below if needed.
 
  // LDVL (input)
  //      The leading dimension of the array VL.  LDVL >= 1; if
  //      JOBVL = 'V', LDVL >= N.
  PetscBLASInt LDVL = VL ? N : 1;
 
  // VR (output) DOUBLE PRECISION array, dimension (LDVR,N)
  //    If JOBVR = 'V', the right eigenvectors v(j) are stored one
  //    after another in the columns of VR, in the same order
  //    as their eigenvalues.
  //    If JOBVR = 'N', VR is not referenced.
  //    If the j-th eigenvalue is real, then v(j) = VR(:,j),
  //    the j-th column of VR.
  //    If the j-th and (j+1)-st eigenvalues form a complex
  //    conjugate pair, then v(j) = VR(:,j) + i*VR(:,j+1) and
  //    v(j+1) = VR(:,j) - i*VR(:,j+1).
  // Will be set below if needed.
 
  // LDVR (input)
  //      The leading dimension of the array VR.  LDVR >= 1; if
  //      JOBVR = 'V', LDVR >= N.
  PetscBLASInt LDVR = VR ? N : 1;
 
  // WORK (workspace/output) double precision array, dimension (MAX(1,LWORK))
  //      On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  //
  // LWORK (input)
  //       The dimension of the array WORK.  LWORK >= max(1,3*N), and
  //       if JOBVL = 'V' or JOBVR = 'V', LWORK >= 4*N.  For good
  //       performance, LWORK must generally be larger.
  //
  //       If LWORK = -1, then a workspace query is assumed; the routine
  //       only calculates the optimal size of the WORK array, returns
  //       this value as the first entry of the WORK array, and no error
  //       message related to LWORK is issued by XERBLA.
  PetscBLASInt LWORK = (VR || VL) ? 4*N : 3*N;
  std::vector<T> WORK(LWORK);
 
  // INFO (output)
  //      = 0:  successful exit
  //      < 0:  if INFO = -i, the i-th argument had an illegal value.
  //      > 0:  if INFO = i, the QR algorithm failed to compute all the
  //            eigenvalues, and no eigenvectors or condition numbers
  //            have been computed; elements 1:ILO-1 and i+1:N of WR
  //            and WI contain eigenvalues which have converged.
  PetscBLASInt INFO = 0;
 
  // Get references to raw data
  std::vector<T> & lambda_real_val = lambda_real.get_values();
  std::vector<T> & lambda_imag_val = lambda_imag.get_values();
 
  // Set up eigenvector storage if necessary.
  T * VR_ptr = nullptr;
  if (VR)
    {
      VR->resize(N, N);
      VR_ptr = VR->get_values().data();
    }
 
  T * VL_ptr = nullptr;
  if (VL)
    {
      VL->resize(N, N);
      VL_ptr = VL->get_values().data();
    }
 
  // Ready to call the Lapack routine through PETSc's interface
  LAPACKgeev_(&JOBVL,
              &JOBVR,
              &N,
              pPS(_val.data()),
              &LDA,
              pPS(lambda_real_val.data()),
              pPS(lambda_imag_val.data()),
              pPS(VL_ptr),
              &LDVL,
              pPS(VR_ptr),
              &LDVR,
              pPS(WORK.data()),
              &LWORK,
              &INFO);
 
  // Check return value for errors
  libmesh_error_msg_if(INFO != 0, "INFO=" << INFO << ", Error during Lapack eigenvalue calculation!");
 
  // If the user requested either right or left eigenvectors, LAPACK
  // has now computed the transpose of the desired matrix, i.e. V^T
  // instead of V.  We could leave this up to user code to handle, but
  // rather than risking getting very unexpected results, we'll just
  // transpose it in place before handing it back.
  if (VR)
    {
      for (auto i : make_range(N))
        for (auto j : make_range(i))
          std::swap((*VR)(i,j), (*VR)(j,i));
    }
 
  if (VL)
    {
      for (auto i : make_range(N))
        for (auto j : make_range(i))
          std::swap((*VL)(i,j), (*VL)(j,i));
    }
}

_left_multiply_transpose()

template<typename T >

template<typename T2 >

void libMesh::DenseMatrix< T >::_left_multiply_transpose ( const DenseMatrix< T2 > &  A )

private

Left multiplies by the transpose of the matrix A which may contain a different numerical type.

Definition at line 112 of file dense_matrix_impl.h.

{
  //Check to see if we are doing (A^T)*A
  if (this == &A)
    {
      //libmesh_here();
      DenseMatrix<T> B(*this);
 
      // Simple but inefficient way
      // return this->left_multiply_transpose(B);
 
      // More efficient, but more code way
      // If A is mxn, the result will be a square matrix of Size n x n.
      const unsigned int n_rows = A.m();
      const unsigned int n_cols = A.n();
 
      // resize() *this and also zero out all entries.
      this->resize(n_cols,n_cols);
 
      // Compute the lower-triangular part
      for (unsigned int i=0; i<n_cols; ++i)
        for (unsigned int j=0; j<=i; ++j)
          for (unsigned int k=0; k<n_rows; ++k) // inner products are over n_rows
            (*this)(i,j) += B(k,i)*B(k,j);
 
      // Copy lower-triangular part into upper-triangular part
      for (unsigned int i=0; i<n_cols; ++i)
        for (unsigned int j=i+1; j<n_cols; ++j)
          (*this)(i,j) = (*this)(j,i);
    }
 
  else
    {
      DenseMatrix<T> B(*this);
 
      this->resize (A.n(), B.n());
 
      libmesh_assert_equal_to (A.m(), B.m());
      libmesh_assert_equal_to (this->m(), A.n());
      libmesh_assert_equal_to (this->n(), B.n());
 
      const unsigned int m_s = A.n();
      const unsigned int p_s = A.m();
      const unsigned int n_s = this->n();
 
      // Do it this way because there is a
      // decent chance (at least for constraint matrices)
      // that A.transpose(i,k) = 0.
      for (unsigned int i=0; i<m_s; i++)
        for (unsigned int k=0; k<p_s; k++)
          if (A.transpose(i,k) != 0.)
            for (unsigned int j=0; j<n_s; j++)
              (*this)(i,j) += A.transpose(i,k)*B(k,j);
    }
}

References libMesh::DenseMatrixBase< T >::m(), and libMesh::DenseMatrixBase< T >::n().

_lu_back_substitute()

template<typename T >

void libMesh::DenseMatrix< T >::_lu_back_substitute ( const DenseVector< T > &  b, DenseVector< T > &  x  ) const

private

Solves the system Ax=b through back substitution.

This function is private since it is only called as part of the implementation of the lu_solve(...) function.

Definition at line 580 of file dense_matrix_impl.h.

{
  const unsigned int
    n_cols = this->n();
 
  libmesh_assert_equal_to (this->m(), n_cols);
  libmesh_assert_equal_to (this->m(), b.size());
 
  x.resize (n_cols);
 
  // A convenient reference to *this
  const DenseMatrix<T> & A = *this;
 
  // Temporary vector storage.  We use this instead of
  // modifying the RHS.
  DenseVector<T> z = b;
 
  // Lower-triangular "top to bottom" solve step, taking into account pivots
  for (unsigned int i=0; i<n_cols; ++i)
    {
      // Swap
      if (_pivots[i] != static_cast<pivot_index_t>(i))
        std::swap( z(i), z(_pivots[i]) );
 
      x(i) = z(i);
 
      for (unsigned int j=0; j<i; ++j)
        x(i) -= A(i,j)*x(j);
 
      x(i) /= A(i,i);
    }
 
  // Upper-triangular "bottom to top" solve step
  const unsigned int last_row = n_cols-1;
 
  for (int i=last_row; i>=0; --i)
    {
      for (int j=i+1; j<static_cast<int>(n_cols); ++j)
        x(i) -= A(i,j)*x(j);
    }
}

References b, and libMesh::DenseVector< T >::resize().

_lu_back_substitute_lapack()

template<typename T >

template LIBMESH_EXPORT void libMesh::DenseMatrix< T >::_lu_back_substitute_lapack ( const DenseVector< T > &  b, DenseVector< T > &  x  )

private

Companion function to _lu_decompose_lapack().

Do not use directly, called through the public lu_solve() interface. This function is logically const in that it does not modify the matrix, but since we are just calling LAPACK routines, it's less const_cast hassle to just declare the function non-const. [ Implementation in dense_matrix_blas_lapack.C ]

Definition at line 901 of file dense_matrix_blas_lapack.C.

{
  // The calling sequence for getrs is:
  // dgetrs(TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
 
  // trans (input)
  //       'n' for no transpose, 't' for transpose
  char TRANS[] = "t";
 
  // N (input)
  //   The order of the matrix A.  N >= 0.
  PetscBLASInt N = this->m();
 
 
  // NRHS (input)
  //      The number of right hand sides, i.e., the number of columns
  //      of the matrix B.  NRHS >= 0.
  PetscBLASInt NRHS = 1;
 
  // A (input) double precision array, dimension (LDA,N)
  //   The factors L and U from the factorization A = P*L*U
  //   as computed by dgetrf.
 
  // LDA (input)
  //     The leading dimension of the array A.  LDA >= max(1,N).
  PetscBLASInt LDA = N;
 
  // ipiv (input) int array, dimension (N)
  //      The pivot indices from DGETRF; for 1<=i<=N, row i of the
  //      matrix was interchanged with row IPIV(i).
  // Here, we pass _pivots.data() which was computed in _lu_decompose_lapack
 
  // B (input/output) double precision array, dimension (LDB,NRHS)
  //   On entry, the right hand side matrix B.
  //   On exit, the solution matrix X.
  // Here, we pass a copy of the rhs vector's data array in x, so that the
  // passed right-hand side b is unmodified.  I don't see a way around this
  // copy if we want to maintain an unmodified rhs in LibMesh.
  x = b;
  std::vector<T> & x_vec = x.get_values();
 
  // We can avoid the copy if we don't care about overwriting the RHS: just
  // pass b to the Lapack routine and then swap with x before exiting
  // std::vector<T> & x_vec = b.get_values();
 
  // LDB (input)
  //     The leading dimension of the array B.  LDB >= max(1,N).
  PetscBLASInt LDB = N;
 
  // INFO (output)
  //      = 0:  successful exit
  //      < 0:  if INFO = -i, the i-th argument had an illegal value
  PetscBLASInt INFO = 0;
 
  // Finally, ready to call the Lapack getrs function
  LAPACKgetrs_(TRANS, &N, &NRHS, pPS(_val.data()), &LDA,
               _pivots.data(), pPS(x_vec.data()), &LDB, &INFO);
 
  // Check return value for errors
  libmesh_error_msg_if(INFO != 0, "INFO=" << INFO << ", Error during Lapack LU solve!");
 
  // Don't do this if you already made a copy of b above
  // Swap b and x.  The solution will then be in x, and whatever was originally
  // in x, maybe garbage, maybe nothing, will be in b.
  // FIXME: Rewrite the LU and Cholesky solves to just take one input, and overwrite
  // the input.  This *should* make user code simpler, as they don't have to create
  // an extra vector just to pass it in to the solve function!
  // b.swap(x);
}

References b, libMesh::DenseVector< T >::get_values(), and libMesh::pPS().

_lu_decompose()

template<typename T >

void libMesh::DenseMatrix< T >::_lu_decompose ( )

private

Form the LU decomposition of the matrix.

This function is private since it is only called as part of the implementation of the lu_solve(...) function.

Definition at line 631 of file dense_matrix_impl.h.

{
  using std::abs;
 
  // If this function was called, there better not be any
  // previous decomposition of the matrix.
  libmesh_assert_equal_to (this->_decomposition_type, NONE);
 
  // Get the matrix size and make sure it is square
  const unsigned int
    n_rows = this->m();
 
  // A convenient reference to *this
  DenseMatrix<T> & A = *this;
 
  _pivots.resize(n_rows);
 
  for (unsigned int i=0; i<n_rows; ++i)
    {
      // Find the pivot row by searching down the i'th column
      _pivots[i] = i;
 
      // abs(complex) must return a Real!
      auto the_max = abs( A(i,i) );
      for (unsigned int j=i+1; j<n_rows; ++j)
        {
          auto candidate_max = abs( A(j,i) );
          if (the_max < candidate_max)
            {
              the_max = candidate_max;
              _pivots[i] = j;
            }
        }
 
      // libMesh::out << "the_max=" << the_max << " found at row " << _pivots[i] << std::endl;
 
      // If the max was found in a different row, interchange rows.
      // Here we interchange the *entire* row, in Gaussian elimination
      // you would only interchange the subrows A(i,j) and A(p(i),j), for j>i
      if (_pivots[i] != static_cast<pivot_index_t>(i))
        {
          for (unsigned int j=0; j<n_rows; ++j)
            std::swap( A(i,j), A(_pivots[i], j) );
        }
 
      // If the max abs entry found is zero, the matrix is singular
      if (A(i,i) == libMesh::zero)
        libmesh_degenerate_mapping_msg ("Matrix A is singular!");
 
      // Scale upper triangle entries of row i by the diagonal entry
      // Note: don't scale the diagonal entry itself!
      const T diag_inv = 1. / A(i,i);
      for (unsigned int j=i+1; j<n_rows; ++j)
        A(i,j) *= diag_inv;
 
      // Update the remaining sub-matrix A[i+1:m][i+1:m]
      // by subtracting off (the diagonal-scaled)
      // upper-triangular part of row i, scaled by the
      // i'th column entry of each row.  In terms of
      // row operations, this is:
      // for each r > i
      //   SubRow(r) = SubRow(r) - A(r,i)*SubRow(i)
      //
      // If we were scaling the i'th column as well, like
      // in Gaussian elimination, this would 'zero' the
      // entry in the i'th column.
      for (unsigned int row=i+1; row<n_rows; ++row)
        for (unsigned int col=i+1; col<n_rows; ++col)
          A(row,col) -= A(row,i) * A(i,col);
 
    } // end i loop
 
  // Set the flag for LU decomposition
  this->_decomposition_type = LU;
}

_lu_decompose_lapack()

template<typename T >

template LIBMESH_EXPORT void libMesh::DenseMatrix< T >::_lu_decompose_lapack ( )

private

Computes an LU factorization of the matrix using the Lapack routine "getrf".

This routine should only be used by the "use_blas_lapack" branch of the lu_solve() function. After the call to this function, the matrix is replaced by its factorized version, and the DecompositionType is set to LU_BLAS_LAPACK. [ Implementation in dense_matrix_blas_lapack.C ]

Definition at line 210 of file dense_matrix_blas_lapack.C.

{
  // If this function was called, there better not be any
  // previous decomposition of the matrix.
  libmesh_assert_equal_to (this->_decomposition_type, NONE);
 
  // The calling sequence for dgetrf is:
  // dgetrf(M, N, A, lda, ipiv, info)
 
  // M (input)
  //   The number of rows of the matrix A.  M >= 0.
  // In C/C++, pass the number of *cols* of A
  PetscBLASInt M = this->n();
 
  // N (input)
  //   The number of columns of the matrix A.  N >= 0.
  // In C/C++, pass the number of *rows* of A
  PetscBLASInt N = this->m();
 
  // A (input/output) double precision array, dimension (LDA,N)
  //   On entry, the M-by-N matrix to be factored.
  //   On exit, the factors L and U from the factorization
  //   A = P*L*U; the unit diagonal elements of L are not stored.
 
  // LDA (input)
  //     The leading dimension of the array A.  LDA >= max(1,M).
  PetscBLASInt LDA = M;
 
  // ipiv (output) integer array, dimension (min(m,n))
  //      The pivot indices; for 1 <= i <= min(m,n), row i of the
  //      matrix was interchanged with row IPIV(i).
  // Here, we pass _pivots.data(), a private class member used to store pivots
  this->_pivots.resize( std::min(M,N) );
 
  // info (output)
  //      = 0:  successful exit
  //      < 0:  if INFO = -i, the i-th argument had an illegal value
  //      > 0:  if INFO = i, U(i,i) is exactly zero. The factorization
  //            has been completed, but the factor U is exactly
  //            singular, and division by zero will occur if it is used
  //            to solve a system of equations.
  PetscBLASInt INFO = 0;
 
  // Ready to call the actual factorization routine through PETSc's interface
  LAPACKgetrf_(&M, &N, pPS(this->_val.data()), &LDA, _pivots.data(),
               &INFO);
 
  // Check return value for errors
  libmesh_error_msg_if(INFO != 0, "INFO=" << INFO << ", Error during Lapack LU factorization!");
 
  // Set the flag for LU decomposition
  this->_decomposition_type = LU_BLAS_LAPACK;
}

References libMesh::pPS().

_matvec_blas()

template<typename T >

template LIBMESH_EXPORT void libMesh::DenseMatrix< T >::_matvec_blas ( T  alpha, T  beta, DenseVector< T > &  dest, const DenseVector< T > &  arg, bool  trans = false  ) const

private

Uses the BLAS GEMV function (through PETSc) to compute.

dest := alpha*A*arg + beta*dest

where alpha and beta are scalars, A is this matrix, and arg and dest are input vectors of appropriate size. If trans is true, the transpose matvec is computed instead. By default, trans==false.

[ Implementation in dense_matrix_blas_lapack.C ]

Definition at line 990 of file dense_matrix_blas_lapack.C.

{
  // Ensure that dest and arg sizes are compatible
  if (!trans)
    {
      // dest  ~ A     * arg
      // (mx1)   (mxn) * (nx1)
      libmesh_error_msg_if((dest.size() != this->m()) || (arg.size() != this->n()),
                           "Improper input argument sizes!");
    }
 
  else // trans == true
    {
      // Ensure that dest and arg are proper size
      // dest  ~ A^T   * arg
      // (nx1)   (nxm) * (mx1)
      libmesh_error_msg_if((dest.size() != this->n()) || (arg.size() != this->m()),
                           "Improper input argument sizes!");
    }
 
  // Calling sequence for dgemv:
  //
  // dgemv(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
 
  // TRANS (input)
  //       't' for transpose, 'n' for non-transpose multiply
  // We store everything in row-major order, so pass the transpose flag for
  // non-transposed matvecs and the 'n' flag for transposed matvecs
  char TRANS[] = "t";
  if (trans)
    TRANS[0] = 'n';
 
  // M (input)
  //   On entry, M specifies the number of rows of the matrix A.
  // In C/C++, pass the number of *cols* of A
  PetscBLASInt M = this->n();
 
  // N (input)
  //   On entry, N specifies the number of columns of the matrix A.
  // In C/C++, pass the number of *rows* of A
  PetscBLASInt N = this->m();
 
  // ALPHA (input)
  // The scalar constant passed to this function
 
  // A (input) double precision array of DIMENSION ( LDA, n ).
  //   Before entry, the leading m by n part of the array A must
  //   contain the matrix of coefficients.
  // The matrix, *this.  Note that _matvec_blas is called from
  // a const function, vector_mult(), and so we have made this function const
  // as well.  Since BLAS knows nothing about const, we have to cast it away
  // now.
  DenseMatrix<T> & a_ref = const_cast<DenseMatrix<T> &> ( *this );
  std::vector<T> & a = a_ref.get_values();
 
  // LDA (input)
  //     On entry, LDA specifies the first dimension of A as declared
  //     in the calling (sub) program. LDA must be at least
  //     max( 1, m ).
  PetscBLASInt LDA = M;
 
  // X (input) double precision array of DIMENSION at least
  //   ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
  //   and at least
  //   ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
  //   Before entry, the incremented array X must contain the
  //   vector x.
  // Here, we must cast away the const-ness of "arg" since BLAS knows
  // nothing about const
  DenseVector<T> & x_ref = const_cast<DenseVector<T> &> ( arg );
  std::vector<T> & x = x_ref.get_values();
 
  // INCX (input)
  //      On entry, INCX specifies the increment for the elements of
  //      X. INCX must not be zero.
  PetscBLASInt INCX = 1;
 
  // BETA (input)
  //      On entry, BETA specifies the scalar beta. When BETA is
  //      supplied as zero then Y need not be set on input.
  // The second scalar constant passed to this function
 
  // Y (input) double precision array of DIMENSION at least
  //   ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
  //   and at least
  //   ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
  //   Before entry with BETA non-zero, the incremented array Y
  //   must contain the vector y. On exit, Y is overwritten by the
  //   updated vector y.
  // The input vector "dest"
  std::vector<T> & y = dest.get_values();
 
  // INCY (input)
  //      On entry, INCY specifies the increment for the elements of
  //      Y. INCY must not be zero.
  PetscBLASInt INCY = 1;
 
  // Finally, ready to call the BLAS function
  BLASgemv_(TRANS, &M, &N, pPS(&alpha), pPS(a.data()), &LDA,
            pPS(x.data()), &INCX, pPS(&beta), pPS(y.data()), &INCY);
}

References libMesh::DenseMatrix< T >::get_values(), libMesh::DenseVector< T >::get_values(), libMesh::pPS(), and libMesh::DenseVector< T >::size().

_multiply_blas()

template<typename T >

template LIBMESH_EXPORT void libMesh::DenseMatrix< T >::_multiply_blas ( const DenseMatrixBase< T > &  other, _BLAS_Multiply_Flag  flag  )

private

The _multiply_blas function computes A <- op(A) * op(B) using BLAS gemm function.

Used in the right_multiply(), left_multiply(), right_multiply_transpose(), and left_multiply_transpose() routines. [ Implementation in dense_matrix_blas_lapack.C ]

Definition at line 37 of file dense_matrix_blas_lapack.C.

{
  int result_size = 0;
 
  // For each case, determine the size of the final result make sure
  // that the inner dimensions match
  switch (flag)
    {
    case LEFT_MULTIPLY:
      {
        result_size = other.m() * this->n();
        if (other.n() == this->m())
          break;
      }
      libmesh_fallthrough();
    case RIGHT_MULTIPLY:
      {
        result_size = other.n() * this->m();
        if (other.m() == this->n())
          break;
      }
      libmesh_fallthrough();
    case LEFT_MULTIPLY_TRANSPOSE:
      {
        result_size = other.n() * this->n();
        if (other.m() == this->m())
          break;
      }
      libmesh_fallthrough();
    case RIGHT_MULTIPLY_TRANSPOSE:
      {
        result_size = other.m() * this->m();
        if (other.n() == this->n())
          break;
      }
      libmesh_fallthrough();
    default:
      libmesh_error_msg("Unknown flag selected or matrices are incompatible for multiplication.");
    }
 
  // For this to work, the passed arg. must actually be a DenseMatrix<T>
  const DenseMatrix<T> * const_that = cast_ptr<const DenseMatrix<T> *>(&other);
 
  // Also, although 'that' is logically const in this BLAS routine,
  // the PETSc BLAS interface does not specify that any of the inputs are
  // const.  To use it, I must cast away const-ness.
  DenseMatrix<T> * that = const_cast<DenseMatrix<T> *> (const_that);
 
  // Initialize A, B pointers for LEFT_MULTIPLY* cases
  DenseMatrix<T> * A = this;
  DenseMatrix<T> * B = that;
 
  // For RIGHT_MULTIPLY* cases, swap the meaning of A and B.
  // Here is a full table of combinations we can pass to BLASgemm, and what the answer is when finished:
  // pass A B   -> (Fortran) -> A^T B^T -> (C++) -> (A^T B^T)^T -> (identity) -> B A   "lt multiply"
  // pass B A   -> (Fortran) -> B^T A^T -> (C++) -> (B^T A^T)^T -> (identity) -> A B   "rt multiply"
  // pass A B^T -> (Fortran) -> A^T B   -> (C++) -> (A^T B)^T   -> (identity) -> B^T A "lt multiply t"
  // pass B^T A -> (Fortran) -> B A^T   -> (C++) -> (B A^T)^T   -> (identity) -> A B^T "rt multiply t"
  if (flag==RIGHT_MULTIPLY || flag==RIGHT_MULTIPLY_TRANSPOSE)
    std::swap(A,B);
 
  // transa, transb values to pass to blas
  char
    transa[] = "n",
    transb[] = "n";
 
  // Integer values to pass to BLAS:
  //
  // M
  // In Fortran, the number of rows of op(A),
  // In the BLAS documentation, typically known as 'M'.
  //
  // In C/C++, we set:
  // M = n_cols(A) if (transa='n')
  //     n_rows(A) if (transa='t')
  PetscBLASInt M = static_cast<PetscBLASInt>( A->n() );
 
  // N
  // In Fortran, the number of cols of op(B), and also the number of cols of C.
  // In the BLAS documentation, typically known as 'N'.
  //
  // In C/C++, we set:
  // N = n_rows(B) if (transb='n')
  //     n_cols(B) if (transb='t')
  PetscBLASInt N = static_cast<PetscBLASInt>( B->m() );
 
  // K
  // In Fortran, the number of cols of op(A), and also
  // the number of rows of op(B). In the BLAS documentation,
  // typically known as 'K'.
  //
  // In C/C++, we set:
  // K = n_rows(A) if (transa='n')
  //     n_cols(A) if (transa='t')
  PetscBLASInt K = static_cast<PetscBLASInt>( A->m() );
 
  // LDA (leading dimension of A). In our cases,
  // LDA is always the number of columns of A.
  PetscBLASInt LDA = static_cast<PetscBLASInt>( A->n() );
 
  // LDB (leading dimension of B).  In our cases,
  // LDB is always the number of columns of B.
  PetscBLASInt LDB = static_cast<PetscBLASInt>( B->n() );
 
  if (flag == LEFT_MULTIPLY_TRANSPOSE)
    {
      transb[0] = 't';
      N = static_cast<PetscBLASInt>( B->n() );
    }
 
  else if (flag == RIGHT_MULTIPLY_TRANSPOSE)
    {
      transa[0] = 't';
      std::swap(M,K);
    }
 
  // LDC (leading dimension of C).  LDC is the
  // number of columns in the solution matrix.
  PetscBLASInt LDC = M;
 
  // Scalar values to pass to BLAS
  //
  // scalar multiplying the whole product AB
  T alpha = 1.;
 
  // scalar multiplying C, which is the original matrix.
  T beta  = 0.;
 
  // Storage for the result
  std::vector<T> result (result_size);
 
  // Finally ready to call the BLAS
  BLASgemm_(transa, transb, &M, &N, &K, pPS(&alpha),
            pPS(A->get_values().data()), &LDA,
            pPS(B->get_values().data()), &LDB, pPS(&beta),
            pPS(result.data()), &LDC);
 
  // Update the relevant dimension for this matrix.
  switch (flag)
    {
    case LEFT_MULTIPLY:            { this->_m = other.m(); break; }
    case RIGHT_MULTIPLY:           { this->_n = other.n(); break; }
    case LEFT_MULTIPLY_TRANSPOSE:  { this->_m = other.n(); break; }
    case RIGHT_MULTIPLY_TRANSPOSE: { this->_n = other.m(); break; }
    default:
      libmesh_error_msg("Unknown flag selected.");
    }
 
  // Swap my data vector with the result
  this->_val.swap(result);
}

References libMesh::DenseMatrix< T >::get_values(), libMesh::DenseMatrixBase< T >::m(), libMesh::DenseMatrixBase< T >::n(), and libMesh::pPS().

_right_multiply_transpose()

template<typename T >

template<typename T2 >

void libMesh::DenseMatrix< T >::_right_multiply_transpose ( const DenseMatrix< T2 > &  A )

private

Right multiplies by the transpose of the matrix A which may contain a different numerical type.

Definition at line 239 of file dense_matrix_impl.h.

{
  //Check to see if we are doing B*(B^T)
  if (this == &B)
    {
      //libmesh_here();
      DenseMatrix<T> A(*this);
 
      // Simple but inefficient way
      // return this->right_multiply_transpose(A);
 
      // More efficient, more code
      // If B is mxn, the result will be a square matrix of Size m x m.
      const unsigned int n_rows = B.m();
      const unsigned int n_cols = B.n();
 
      // resize() *this and also zero out all entries.
      this->resize(n_rows,n_rows);
 
      // Compute the lower-triangular part
      for (unsigned int i=0; i<n_rows; ++i)
        for (unsigned int j=0; j<=i; ++j)
          for (unsigned int k=0; k<n_cols; ++k) // inner products are over n_cols
            (*this)(i,j) += A(i,k)*A(j,k);
 
      // Copy lower-triangular part into upper-triangular part
      for (unsigned int i=0; i<n_rows; ++i)
        for (unsigned int j=i+1; j<n_rows; ++j)
          (*this)(i,j) = (*this)(j,i);
    }
 
  else
    {
      DenseMatrix<T> A(*this);
 
      this->resize (A.m(), B.m());
 
      libmesh_assert_equal_to (A.n(), B.n());
      libmesh_assert_equal_to (this->m(), A.m());
      libmesh_assert_equal_to (this->n(), B.m());
 
      const unsigned int m_s = A.m();
      const unsigned int p_s = A.n();
      const unsigned int n_s = this->n();
 
      // Do it this way because there is a
      // decent chance (at least for constraint matrices)
      // that B.transpose(k,j) = 0.
      for (unsigned int j=0; j<n_s; j++)
        for (unsigned int k=0; k<p_s; k++)
          if (B.transpose(k,j) != 0.)
            for (unsigned int i=0; i<m_s; i++)
              (*this)(i,j) += A(i,k)*B.transpose(k,j);
    }
}

References libMesh::DenseMatrixBase< T >::m(), and libMesh::DenseMatrixBase< T >::n().

_svd_helper()

template<typename T >

template LIBMESH_EXPORT void libMesh::DenseMatrix< T >::_svd_helper ( char  JOBU, char  JOBVT, std::vector< Real > &  sigma_val, std::vector< Number > &  U_val, std::vector< Number > &  VT_val  )

private

Helper function that actually performs the SVD.

[ Implementation in dense_matrix_blas_lapack.C ]

Definition at line 384 of file dense_matrix_blas_lapack.C.

{
 
  // M (input)
  //   The number of rows of the matrix A.  M >= 0.
  // In C/C++, pass the number of *cols* of A
  PetscBLASInt M = this->n();
 
  // N (input)
  //   The number of columns of the matrix A.  N >= 0.
  // In C/C++, pass the number of *rows* of A
  PetscBLASInt N = this->m();
 
  PetscBLASInt min_MN = (M < N) ? M : N;
  PetscBLASInt max_MN = (M > N) ? M : N;
 
  // A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
  //   On entry, the M-by-N matrix A.
  //   On exit,
  //   if JOBU  = 'O', A is overwritten with the first min(m,n)
  //                   columns of U (the left singular vectors,
  //                   stored columnwise);
  //   if JOBVT = 'O', A is overwritten with the first min(m,n)
  //                   rows of V**T (the right singular vectors,
  //                   stored rowwise);
  //   if JOBU != 'O' and JOBVT != 'O', the contents of A are destroyed.
 
  // LDA (input)
  //     The leading dimension of the array A.  LDA >= max(1,M).
  PetscBLASInt LDA = M;
 
  // S (output) DOUBLE PRECISION array, dimension (min(M,N))
  //   The singular values of A, sorted so that S(i) >= S(i+1).
  sigma_val.resize( min_MN );
 
  // LDU (input)
  //     The leading dimension of the array U.  LDU >= 1; if
  //     JOBU = 'S' or 'A', LDU >= M.
  PetscBLASInt LDU = M;
 
  // U (output) DOUBLE PRECISION array, dimension (LDU,UCOL)
  //   (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'.
  //   If JOBU = 'A', U contains the M-by-M orthogonal matrix U;
  //   if JOBU = 'S', U contains the first min(m,n) columns of U
  //   (the left singular vectors, stored columnwise);
  //   if JOBU = 'N' or 'O', U is not referenced.
  if (JOBU == 'S')
    U_val.resize( LDU*min_MN );
  else
    U_val.resize( LDU*M );
 
  // LDVT (input)
  //      The leading dimension of the array VT.  LDVT >= 1; if
  //      JOBVT = 'A', LDVT >= N; if JOBVT = 'S', LDVT >= min(M,N).
  PetscBLASInt LDVT = N;
  if (JOBVT == 'S')
    LDVT = min_MN;
 
  // VT (output) DOUBLE PRECISION array, dimension (LDVT,N)
  //    If JOBVT = 'A', VT contains the N-by-N orthogonal matrix
  //    V**T;
  //    if JOBVT = 'S', VT contains the first min(m,n) rows of
  //    V**T (the right singular vectors, stored rowwise);
  //    if JOBVT = 'N' or 'O', VT is not referenced.
  VT_val.resize( LDVT*N );
 
  // LWORK (input)
  //       The dimension of the array WORK.
  //       LWORK >= MAX(1,3*MIN(M,N)+MAX(M,N),5*MIN(M,N)).
  //       For good performance, LWORK should generally be larger.
  //
  //       If LWORK = -1, then a workspace query is assumed; the routine
  //       only calculates the optimal size of the WORK array, returns
  //       this value as the first entry of the WORK array, and no error
  //       message related to LWORK is issued by XERBLA.
  PetscBLASInt larger = (3*min_MN+max_MN > 5*min_MN) ? 3*min_MN+max_MN : 5*min_MN;
  PetscBLASInt LWORK  = (larger > 1) ? larger : 1;
 
 
  // WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
  //      On exit, if INFO = 0, WORK(1) returns the optimal LWORK;
  //      if INFO > 0, WORK(2:MIN(M,N)) contains the unconverged
  //      superdiagonal elements of an upper bidiagonal matrix B
  //      whose diagonal is in S (not necessarily sorted). B
  //      satisfies A = U * B * VT, so it has the same singular values
  //      as A, and singular vectors related by U and VT.
  std::vector<Number> WORK( LWORK );
 
  // INFO (output)
  //      = 0:  successful exit.
  //      < 0:  if INFO = -i, the i-th argument had an illegal value.
  //      > 0:  if DBDSQR did not converge, INFO specifies how many
  //            superdiagonals of an intermediate bidiagonal form B
  //            did not converge to zero. See the description of WORK
  //            above for details.
  PetscBLASInt INFO = 0;
 
  // Ready to call the actual factorization routine through PETSc's interface.
#ifdef LIBMESH_USE_REAL_NUMBERS
  // Note that the call to LAPACKgesvd_ may modify _val
  LAPACKgesvd_(&JOBU, &JOBVT, &M, &N, pPS(_val.data()), &LDA,
               pPS(sigma_val.data()), pPS(U_val.data()), &LDU,
               pPS(VT_val.data()), &LDVT, pPS(WORK.data()), &LWORK,
               &INFO);
#else
  // When we have LIBMESH_USE_COMPLEX_NUMBERS then we must pass an array of Complex
  // numbers to LAPACKgesvd_, but _val may contain Reals so we copy to Number below to
  // handle both the real-valued and complex-valued cases.
  std::vector<Number> val_copy(_val.size());
  for (auto i : make_range(_val.size()))
    val_copy[i] = _val[i];
 
  std::vector<Real> RWORK(5 * min_MN);
  LAPACKgesvd_(&JOBU, &JOBVT, &M, &N, val_copy.data(), &LDA, sigma_val.data(), U_val.data(),
               &LDU, VT_val.data(), &LDVT, WORK.data(), &LWORK, RWORK.data(), &INFO);
#endif
 
  // Check return value for errors
  libmesh_error_msg_if(INFO != 0, "INFO=" << INFO << ", Error during Lapack SVD calculation!");
}

References libMesh::make_range(), and libMesh::pPS().

_svd_lapack() [1/2]

template<typename T >

template LIBMESH_EXPORT void libMesh::DenseMatrix< T >::_svd_lapack ( DenseVector< Real > &  sigma )

private

Computes an SVD of the matrix using the Lapack routine "getsvd".

[ Implementation in dense_matrix_blas_lapack.C ]

Definition at line 277 of file dense_matrix_blas_lapack.C.

{
  // The calling sequence for dgetrf is:
  // DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, INFO )
 
  // JOBU (input)
  //      Specifies options for computing all or part of the matrix U:
  //      = 'A':  all M columns of U are returned in array U:
  //      = 'S':  the first min(m,n) columns of U (the left singular
  //              vectors) are returned in the array U;
  //      = 'O':  the first min(m,n) columns of U (the left singular
  //              vectors) are overwritten on the array A;
  //      = 'N':  no columns of U (no left singular vectors) are
  //              computed.
  char JOBU = 'N';
 
  // JOBVT (input)
  //       Specifies options for computing all or part of the matrix
  //       V**T:
  //       = 'A':  all N rows of V**T are returned in the array VT;
  //       = 'S':  the first min(m,n) rows of V**T (the right singular
  //               vectors) are returned in the array VT;
  //       = 'O':  the first min(m,n) rows of V**T (the right singular
  //               vectors) are overwritten on the array A;
  //       = 'N':  no rows of V**T (no right singular vectors) are
  //               computed.
  char JOBVT = 'N';
 
  std::vector<Real> sigma_val;
  std::vector<Number> U_val;
  std::vector<Number> VT_val;
 
  _svd_helper(JOBU, JOBVT, sigma_val, U_val, VT_val);
 
  // Copy the singular values into sigma, ignore U_val and VT_val
  sigma.resize(cast_int<unsigned int>(sigma_val.size()));
  for (auto i : make_range(sigma.size()))
    sigma(i) = sigma_val[i];
}

References libMesh::make_range(), libMesh::DenseVector< T >::resize(), and libMesh::DenseVector< T >::size().

_svd_lapack() [2/2]

template<typename T >

template LIBMESH_EXPORT void libMesh::DenseMatrix< T >::_svd_lapack ( DenseVector< Real > &  sigma, DenseMatrix< Number > &  U, DenseMatrix< Number > &  VT  )

private

Computes a "reduced" SVD of the matrix using the Lapack routine "getsvd".

[ Implementation in dense_matrix_blas_lapack.C ]

Definition at line 318 of file dense_matrix_blas_lapack.C.

{
  // The calling sequence for dgetrf is:
  // DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, INFO )
 
  // JOBU (input)
  //      Specifies options for computing all or part of the matrix U:
  //      = 'A':  all M columns of U are returned in array U:
  //      = 'S':  the first min(m,n) columns of U (the left singular
  //              vectors) are returned in the array U;
  //      = 'O':  the first min(m,n) columns of U (the left singular
  //              vectors) are overwritten on the array A;
  //      = 'N':  no columns of U (no left singular vectors) are
  //              computed.
  char JOBU = 'S';
 
  // JOBVT (input)
  //       Specifies options for computing all or part of the matrix
  //       V**T:
  //       = 'A':  all N rows of V**T are returned in the array VT;
  //       = 'S':  the first min(m,n) rows of V**T (the right singular
  //               vectors) are returned in the array VT;
  //       = 'O':  the first min(m,n) rows of V**T (the right singular
  //               vectors) are overwritten on the array A;
  //       = 'N':  no rows of V**T (no right singular vectors) are
  //               computed.
  char JOBVT = 'S';
 
  // Note: Lapack is going to compute the singular values of A^T.  If
  // A=U * S * V^T, then A^T = V * S * U^T, which means that the
  // values returned in the "U_val" array actually correspond to the
  // entries of the V matrix, and the values returned in the VT_val
  // array actually correspond to the entries of U^T.  Therefore, we
  // pass VT in the place of U and U in the place of VT below!
  std::vector<Real> sigma_val;
  int M = this->n();
  int N = this->m();
  int min_MN = (M < N) ? M : N;
 
  // Size user-provided storage appropriately. Inside svd_helper:
  // U_val is sized to (M x min_MN)
  // VT_val is sized to (min_MN x N)
  // So, we set up U to have the shape of "VT_val^T", and VT to
  // have the shape of "U_val^T".
  //
  // Finally, since the results are stored in column-major order by
  // Lapack, but we actually want the transpose of what Lapack
  // returns, this means (conveniently) that we don't even have to
  // copy anything after the call to _svd_helper, it should already be
  // in the correct order!
  U.resize(N, min_MN);
  VT.resize(min_MN, M);
 
  _svd_helper(JOBU, JOBVT, sigma_val, VT.get_values(), U.get_values());
 
  // Copy the singular values into sigma.
  sigma.resize(cast_int<unsigned int>(sigma_val.size()));
  for (auto i : make_range(sigma.size()))
    sigma(i) = sigma_val[i];
}

References libMesh::DenseMatrix< T >::get_values(), libMesh::make_range(), libMesh::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::resize(), and libMesh::DenseVector< T >::size().

_svd_solve_lapack()

template<typename T >

template LIBMESH_EXPORT void libMesh::DenseMatrix< T >::_svd_solve_lapack ( const DenseVector< T > &  rhs, DenseVector< T > &  x, Real  rcond  ) const

private

Called by svd_solve(rhs).

Definition at line 529 of file dense_matrix_blas_lapack.C.

{
  // Since BLAS is expecting column-major storage, we first need to
  // make a transposed copy of *this, then pass it to the gelss
  // routine instead of the original.  This extra copy is kind of a
  // bummer, it might be better if we could use the full SVD to
  // compute the least-squares solution instead...  Note that it isn't
  // completely terrible either, since A_trans gets overwritten by
  // Lapack, and we usually would end up making a copy of A outside
  // the function call anyway.
  DenseMatrix<T> A_trans;
  this->get_transpose(A_trans);
 
  // M
  // The number of rows of the input matrix. M >= 0.
  // This is actually the number of *columns* of A_trans.
  PetscBLASInt M = A_trans.n();
 
  // N
  // The number of columns of the matrix A. N >= 0.
  // This is actually the number of *rows* of A_trans.
  PetscBLASInt N = A_trans.m();
 
  // We'll use the min and max of (M,N) several times below.
  PetscBLASInt max_MN = std::max(M,N);
  PetscBLASInt min_MN = std::min(M,N);
 
  // NRHS
  // The number of right hand sides, i.e., the number of columns
  // of the matrices B and X. NRHS >= 0.
  // This could later be generalized to solve for multiple right-hand
  // sides...
  PetscBLASInt NRHS = 1;
 
  // A is double precision array, dimension (LDA,N)
  // On entry, the M-by-N matrix A.
  // On exit, the first min(m,n) rows of A are overwritten with
  // its right singular vectors, stored rowwise.
  //
  // The data vector that will be passed to Lapack.
  std::vector<T> & A_trans_vals = A_trans.get_values();
 
  // LDA
  // The leading dimension of the array A.  LDA >= max(1,M).
  PetscBLASInt LDA = M;
 
  // B is double precision array, dimension (LDB,NRHS)
  // On entry, the M-by-NRHS right hand side matrix B.
  // On exit, B is overwritten by the N-by-NRHS solution
  // matrix X.  If m >= n and RANK = n, the residual
  // sum-of-squares for the solution in the i-th column is given
  // by the sum of squares of elements n+1:m in that column.
  //
  // Since we don't want the user's rhs vector to be overwritten by
  // the solution, we copy the rhs values into the solution vector "x"
  // now.  x needs to be long enough to hold both the (Nx1) solution
  // vector or the (Mx1) rhs, so size it to the max of those.
  x.resize(max_MN);
  for (auto i : make_range(rhs.size()))
    x(i) = rhs(i);
 
  // Make the syntax below simpler by grabbing a reference to this array.
  std::vector<T> & B = x.get_values();
 
  // LDB
  // The leading dimension of the array B. LDB >= max(1,max(M,N)).
  PetscBLASInt LDB = x.size();
 
  // S is double precision array, dimension (min(M,N))
  // The singular values of A in decreasing order.
  // The condition number of A in the 2-norm = S(1)/S(min(m,n)).
  std::vector<T> S(min_MN);
 
  // RCOND
  // Used to determine the effective rank of A.  Singular values
  // S(i) <= RCOND*S(1) are treated as zero.  If RCOND < 0, machine
  // precision is used instead.
  PetscScalar RCOND = PS(rcond);
 
  // RANK
  // The effective rank of A, i.e., the number of singular values
  // which are greater than RCOND*S(1).
  PetscBLASInt RANK = 0;
 
  // LWORK
  // The dimension of the array WORK. LWORK >= 1, and also:
  // LWORK >= 3*min(M,N) + max( 2*min(M,N), max(M,N), NRHS )
  // For good performance, LWORK should generally be larger.
  //
  // If LWORK = -1, then a workspace query is assumed; the routine
  // only calculates the optimal size of the WORK array, returns
  // this value as the first entry of the WORK array, and no error
  // message related to LWORK is issued by XERBLA.
  //
  // The factor of 3/2 is arbitrary and is used to satisfy the "should
  // generally be larger" clause.
  PetscBLASInt LWORK = (3*min_MN + std::max(2*min_MN, std::max(max_MN, NRHS))) * 3/2;
 
  // WORK is double precision array, dimension (MAX(1,LWORK))
  // On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  std::vector<T> WORK(LWORK);
 
  // INFO
  // = 0:  successful exit
  // < 0:  if INFO = -i, the i-th argument had an illegal value.
  // > 0:  the algorithm for computing the SVD failed to converge;
  //       if INFO = i, i off-diagonal elements of an intermediate
  //       bidiagonal form did not converge to zero.
  PetscBLASInt INFO = 0;
 
  // LAPACKgelss_(const PetscBLASInt *, // M
  //              const PetscBLASInt *, // N
  //              const PetscBLASInt *, // NRHS
  //              PetscScalar *,        // A
  //              const PetscBLASInt *, // LDA
  //              PetscScalar *,        // B
  //              const PetscBLASInt *, // LDB
  //              PetscReal *,          // S(out) = singular values of A in increasing order
  //              const PetscReal *,    // RCOND = tolerance for singular values
  //              PetscBLASInt *,       // RANK(out) = number of "non-zero" singular values
  //              PetscScalar *,        // WORK
  //              const PetscBLASInt *, // LWORK
  //              PetscBLASInt *);      // INFO
  LAPACKgelss_(&M, &N, &NRHS, pPS(A_trans_vals.data()), &LDA,
               pPS(B.data()), &LDB, pPS(S.data()), &RCOND, &RANK,
               pPS(WORK.data()), &LWORK, &INFO);
 
  // Check for errors in the Lapack call
  libmesh_error_msg_if(INFO < 0, "Error, argument " << -INFO << " to LAPACKgelss_ had an illegal value.");
  libmesh_error_msg_if(INFO > 0, "The algorithm for computing the SVD failed to converge!");
 
  // Debugging: print singular values and information about condition number:
  // libMesh::err << "RCOND=" << RCOND << std::endl;
  // libMesh::err << "Singular values: " << std::endl;
  // for (auto i : index_range(S))
  //   libMesh::err << S[i] << std::endl;
  // libMesh::err << "The condition number of A is approximately: " << S[0]/S.back() << std::endl;
 
  // Lapack has already written the solution into B, but it will be
  // the wrong size for non-square problems, so we need to resize it
  // correctly.  The size of the solution vector should be the number
  // of columns of the original A matrix.  Unfortunately, resizing a
  // DenseVector currently also zeros it out (unlike a std::vector) so
  // we'll resize the underlying storage directly (the size is not
  // stored independently elsewhere).
  x.get_values().resize(this->n());
}

References libMesh::DenseMatrix< T >::get_values(), libMesh::DenseVector< T >::get_values(), libMesh::DenseMatrixBase< T >::m(), libMesh::make_range(), libMesh::DenseMatrixBase< T >::n(), libMesh::pPS(), libMesh::PS(), libMesh::DenseVector< T >::resize(), and libMesh::DenseVector< T >::size().

add() [1/2]

template<typename T >

template<typename T2 , typename T3 >

std::enable_if< ScalarTraits< T2 >::value, void >::type libMesh::DenseMatrix< T >::add ( const T2  factor, const DenseMatrix< T3 > &  mat  )

inline

Adds factor times mat to this matrix.

Returns

A reference to *this.

Definition at line 1018 of file dense_matrix.h.

{
  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);
}

References libMesh::DenseMatrixBase< T >::m(), libMesh::make_range(), and libMesh::DenseMatrixBase< T >::n().

Referenced by assemble_wave(), libMesh::FEMSystem::assembly(), libMesh::TransientRBEvaluation::rb_solve(), and libMesh::RBEvaluation::rb_solve().

add() [2/2]

template<typename T >

template<typename T2 , typename T3 >

std::enable_if< ScalarTraits< T2 >::value, void >::type libMesh::DenseMatrixBase< T >::add ( const T2  factor, const DenseMatrixBase< T3 > &  mat  )

inlineinherited

Adds factor to every element in the matrix.

This should only work if T += T2 * T3 is valid C++ and if T2 is scalar. Return type is void

Definition at line 195 of file dense_matrix_base.h.

{
  libmesh_assert_equal_to (this->m(), mat.m());
  libmesh_assert_equal_to (this->n(), mat.n());
 
  for (auto j : make_range(this->n()))
    for (auto i : make_range(this->m()))
      this->el(i,j) += factor*mat.el(i,j);
}

References libMesh::DenseMatrixBase< T >::el(), libMesh::DenseMatrixBase< T >::m(), libMesh::make_range(), and libMesh::DenseMatrixBase< T >::n().

cholesky_solve()

template<typename T >

template<typename T2 >

template LIBMESH_EXPORT void libMesh::DenseMatrix< T >::cholesky_solve	(	const DenseVector< T2 > &	b,
		DenseVector< T2 > &	x
	)

For symmetric positive definite (SPD) matrices.

A Cholesky factorization of A such that A = L L^T is about twice as fast as a standard LU factorization. Therefore you can use this method if you know a-priori that the matrix is SPD. If the matrix is not SPD, an error is generated. One nice property of Cholesky decompositions is that they do not require pivoting for stability.

Important note: once you call cholesky_solve(), you must not modify the entries of the matrix via calls to operator(i,j) and call cholesky_solve() again without first calling either zero() or resize(), otherwise the code will skip computing the decomposition of the matrix and go directly to the back substitution step. This is done on purpose for efficiency, so that the same decomposition can be used with multiple right-hand sides, but it does also make it possible to "shoot yourself in the foot", so be careful!

Note

This method may also be used when A is real-valued and x and b are complex-valued.

Definition at line 844 of file dense_matrix_impl.h.

{
  // Check for a previous decomposition
  switch(this->_decomposition_type)
    {
    case NONE:
      {
        this->_cholesky_decompose ();
        break;
      }
 
    case CHOLESKY:
      {
        // Already factored, just need to call back_substitute.
        break;
      }
 
    default:
      libmesh_error_msg("Error! This matrix already has a different decomposition...");
    }
 
  // Perform back substitution
  this->_cholesky_back_substitute (b, x);
}

References b.

Referenced by libMesh::FEGenericBase< OutputType >::coarsened_dof_values(), libMesh::FEGenericBase< OutputType >::compute_periodic_constraints(), libMesh::FEGenericBase< OutputType >::compute_proj_constraints(), LinearElasticity::compute_stresses(), libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction >::SubProjector::construct_projection(), libMesh::BoundaryProjectSolution::operator()(), and libMesh::HPCoarsenTest::select_refinement().

condense() [1/2]

template<typename T >

void libMesh::DenseMatrix< T >::condense ( const unsigned int  i, const unsigned int  j, const T  val, DenseVector< T > &  rhs  )

inline

Condense-out the (i,j) entry of the matrix, forcing it to take on the value val.

This is useful in numerical simulations for applying boundary conditions. Preserves the symmetry of the matrix.

Definition at line 395 of file dense_matrix.h.

  { DenseMatrixBase<T>::condense (i, j, val, rhs); }

References libMesh::DenseMatrixBase< T >::condense().

condense() [2/2]

template<typename T >

void libMesh::DenseMatrixBase< T >::condense ( const unsigned int  i, const unsigned int  j, const T  val, DenseVectorBase< T > &  rhs  )

protectedinherited

Condense-out the (i,j) entry of the matrix, forcing it to take on the value val.

This is useful in numerical simulations for applying boundary conditions. Preserves the symmetry of the matrix.

Definition at line 73 of file dense_matrix_base_impl.h.

  {
    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.el(i) -= this->el(i,jv)*val;
        this->el(i,jv) = 0.;
      }
 
    // zero the row
    for (auto j : make_range(this->n()))
      this->el(iv,j) = 0.;
 
    this->el(iv,jv) = 1.;
    rhs.el(iv) = val;
 
  }

References libMesh::DenseVectorBase< T >::el(), libMesh::make_range(), and libMesh::DenseVectorBase< T >::size().

Referenced by libMesh::DenseMatrix< T >::condense().

det()

template<typename T >

T libMesh::DenseMatrix< T >::det

(

)

Returns

The determinant of the matrix.

Note

Implemented by computing an LU decomposition and then taking the product of the diagonal terms. Therefore this is a non-const method which modifies the entries of the matrix.

Definition at line 783 of file dense_matrix_impl.h.

{
  switch(this->_decomposition_type)
    {
    case NONE:
      {
        // First LU decompose the matrix.
        // Note that the lu_decompose routine will check to see if the
        // matrix is square so we don't worry about it.
        if (this->use_blas_lapack)
          this->_lu_decompose_lapack();
        else
          this->_lu_decompose ();
      }
    case LU:
    case LU_BLAS_LAPACK:
      {
        // Already decomposed, don't do anything
        break;
      }
    default:
      libmesh_error_msg("Error! Can't compute the determinant under the current decomposition.");
    }
 
  // A variable to keep track of the running product of diagonal terms.
  T determinant = 1.;
 
  // Loop over diagonal terms, computing the product.  In practice,
  // be careful because this value could easily become too large to
  // fit in a double or float.  To be safe, one should keep track of
  // the power (of 10) of the determinant in a separate variable
  // and maintain an order 1 value for the determinant itself.
  unsigned int n_interchanges = 0;
  for (auto i : make_range(this->m()))
    {
      if (this->_decomposition_type==LU)
        if (_pivots[i] != static_cast<pivot_index_t>(i))
          n_interchanges++;
 
      // Lapack pivots are 1-based!
      if (this->_decomposition_type==LU_BLAS_LAPACK)
        if (_pivots[i] != static_cast<pivot_index_t>(i+1))
          n_interchanges++;
 
      determinant *= (*this)(i,i);
    }
 
  // Compute sign of determinant, depends on number of row interchanges!
  // The sign should be (-1)^{n}, where n is the number of interchanges.
  Real sign = n_interchanges % 2 == 0 ? 1. : -1.;
 
  return sign*determinant;
}

diagonal()

template<typename T >

DenseVector< T > libMesh::DenseMatrixBase< T >::diagonal ( ) const

inherited

Return the matrix diagonal.

Definition at line 37 of file dense_matrix_base_impl.h.

  {
    DenseVector<T> ret(_m);
    for (decltype(_m) i = 0; i < _m; ++i)
      ret(i) = el(i, i);
    return ret;
  }

Referenced by libMesh::DiagonalMatrix< T >::add_matrix().

el() [1/2]

template<typename T >

virtual T libMesh::DenseMatrix< T >::el ( const unsigned int  i, const unsigned int  j  ) const

inlinefinaloverridevirtual

Returns

The (i,j) element of the matrix. Since internal data representations may differ, you must redefine this function.

Implements libMesh::DenseMatrixBase< T >.

Definition at line 126 of file dense_matrix.h.

  { return (*this)(i,j); }

Referenced by libMesh::RBConstruction::add_scaled_matrix_and_vector(), libMesh::TransientRBConstruction::enrich_RB_space(), libMesh::RBEIMConstruction::train_eim_approximation_with_POD(), and libMesh::RBConstruction::train_reduced_basis_with_POD().

el() [2/2]

template<typename T >

virtual T & libMesh::DenseMatrix< T >::el ( const unsigned int  i, const unsigned int  j  )

inlinefinaloverridevirtual

Returns

The (i,j) element of the matrix as a writable reference. Since internal data representations may differ, you must redefine this function.

Implements libMesh::DenseMatrixBase< T >.

Definition at line 130 of file dense_matrix.h.

  { return (*this)(i,j); }

evd()

template<typename T >

void libMesh::DenseMatrix< T >::evd	(	DenseVector< T > &	lambda_real,
		DenseVector< T > &	lambda_imag
	)

Compute the eigenvalues (both real and imaginary parts) of a general matrix.

Warning: the contents of *this are overwritten by this function!

The implementation requires the LAPACKgeev_ function which is wrapped by PETSc.

Definition at line 739 of file dense_matrix_impl.h.

{
  // We use the LAPACK eigenvalue problem implementation
  _evd_lapack(lambda_real, lambda_imag);
}

evd_left()

template<typename T >

void libMesh::DenseMatrix< T >::evd_left	(	DenseVector< T > &	lambda_real,
		DenseVector< T > &	lambda_imag,
		DenseMatrix< T > &	VL
	)

Compute the eigenvalues (both real and imaginary parts) and left eigenvectors of a general matrix, \( A \).

Warning: the contents of *this are overwritten by this function!

The left eigenvector \( u_j \) of \( A \) satisfies: \( u_j^H A = lambda_j u_j^H \) where \( u_j^H \) denotes the conjugate-transpose of \( u_j \).

If the j-th and (j+1)-st eigenvalues form a complex conjugate pair, then the j-th and (j+1)-st columns of VL "share" their real-valued storage in the following way: u_j = VL(:,j) + i*VL(:,j+1) and u_{j+1} = VL(:,j) - i*VL(:,j+1).

The implementation requires the LAPACKgeev_ routine which is provided by PETSc.

Definition at line 749 of file dense_matrix_impl.h.

{
  // We use the LAPACK eigenvalue problem implementation
  _evd_lapack(lambda_real, lambda_imag, &VL, nullptr);
}

evd_left_and_right()

template<typename T >

void libMesh::DenseMatrix< T >::evd_left_and_right	(	DenseVector< T > &	lambda_real,
		DenseVector< T > &	lambda_imag,
		DenseMatrix< T > &	VL,
		DenseMatrix< T > &	VR
	)

Compute the eigenvalues (both real and imaginary parts) as well as the left and right eigenvectors of a general matrix.

Warning: the contents of *this are overwritten by this function!

See the documentation of the evd_left() and evd_right() functions for more information. The implementation requires the LAPACKgeev_ routine which is provided by PETSc.

Definition at line 771 of file dense_matrix_impl.h.

{
  // We use the LAPACK eigenvalue problem implementation
  _evd_lapack(lambda_real, lambda_imag, &VL, &VR);
}

Referenced by DenseMatrixTest::testEVD_helper().

evd_right()

template<typename T >

void libMesh::DenseMatrix< T >::evd_right	(	DenseVector< T > &	lambda_real,
		DenseVector< T > &	lambda_imag,
		DenseMatrix< T > &	VR
	)

Compute the eigenvalues (both real and imaginary parts) and right eigenvectors of a general matrix, \( A \).

Warning: the contents of *this are overwritten by this function!

The right eigenvector \( v_j \) of \( A \) satisfies: \( A v_j = lambda_j v_j \) where \( lambda_j \) is its corresponding eigenvalue.

Note

If the j-th and (j+1)-st eigenvalues form a complex conjugate pair, then the j-th and (j+1)-st columns of VR "share" their real-valued storage in the following way: v_j = VR(:,j) + i*VR(:,j+1) and v_{j+1} = VR(:,j) - i*VR(:,j+1).

The implementation requires the LAPACKgeev_ routine which is provided by PETSc.

Definition at line 760 of file dense_matrix_impl.h.

{
  // We use the LAPACK eigenvalue problem implementation
  _evd_lapack(lambda_real, lambda_imag, nullptr, &VR);
}

get_principal_submatrix() [1/2]

template<typename T >

void libMesh::DenseMatrix< T >::get_principal_submatrix	(	unsigned int	sub_m,
		DenseMatrix< T > &	dest
	)		const

Put the sub_m x sub_m principal submatrix into dest.

Definition at line 500 of file dense_matrix_impl.h.

{
  get_principal_submatrix(sub_m, sub_m, dest);
}

get_principal_submatrix() [2/2]

template<typename T >

void libMesh::DenseMatrix< T >::get_principal_submatrix	(	unsigned int	sub_m,
		unsigned int	sub_n,
		DenseMatrix< T > &	dest
	)		const

Put the sub_m x sub_n principal submatrix into dest.

Definition at line 470 of file dense_matrix_impl.h.

{
  libmesh_assert( (sub_m <= this->m()) && (sub_n <= this->n()) );
 
  dest.resize(sub_m, sub_n);
  for (unsigned int i=0; i<sub_m; i++)
    for (unsigned int j=0; j<sub_n; j++)
      dest(i,j) = (*this)(i,j);
}

References libMesh::libmesh_assert(), and libMesh::DenseMatrix< T >::resize().

Referenced by libMesh::RBEIMConstruction::compute_max_eim_error(), libMesh::RBEIMEvaluation::rb_eim_solve(), libMesh::RBEIMEvaluation::rb_eim_solves(), libMesh::TransientRBEvaluation::rb_solve(), libMesh::RBEvaluation::rb_solve(), and libMesh::TransientRBConstruction::update_RB_initial_condition_all_N().

get_transpose()

template<typename T >

void libMesh::DenseMatrix< T >::get_transpose

(

DenseMatrix< T > &

dest

)

const

Put the tranposed matrix into dest.

Definition at line 508 of file dense_matrix_impl.h.

{
  dest.resize(this->n(), this->m());
 
  for (auto i : make_range(dest.m()))
    for (auto j : make_range(dest.n()))
      dest(i,j) = (*this)(j,i);
}

References libMesh::DenseMatrixBase< T >::m(), libMesh::make_range(), libMesh::DenseMatrixBase< T >::n(), and libMesh::DenseMatrix< T >::resize().

Referenced by libMesh::RBConstruction::add_scaled_matrix_and_vector().

get_values() [1/2]

template<typename T >

std::vector< T > & libMesh::DenseMatrix< T >::get_values ( )

inline

Returns

A reference to the underlying data storage vector.

This should be used with caution (i.e. one should not change the size of the vector, etc.) but is useful for interoperating with low level BLAS routines which expect a simple array.

Definition at line 382 of file dense_matrix.h.

{ return _val; }

References libMesh::DenseMatrix< T >::_val.

Referenced by libMesh::DenseMatrix< T >::_matvec_blas(), libMesh::DenseMatrix< T >::_multiply_blas(), libMesh::DenseMatrix< T >::_svd_lapack(), libMesh::DenseMatrix< T >::_svd_solve_lapack(), libMesh::PetscMatrix< T >::add_block_matrix(), libMesh::PetscMatrix< T >::add_matrix(), libMesh::EpetraMatrix< T >::add_matrix(), and OverlappingCouplingGhostingTest::run_sparsity_pattern_test().

get_values() [2/2]

template<typename T >

const std::vector< T > & libMesh::DenseMatrix< T >::get_values ( ) const

inline

Returns

A constant reference to the underlying data storage vector.

Definition at line 387 of file dense_matrix.h.

{ return _val; }

References libMesh::DenseMatrix< T >::_val.

l1_norm()

template<typename T >

auto libMesh::DenseMatrix< T >::l1_norm ( ) const

inline

Returns

The l1-norm of the matrix, that is, the max column sum:

\( |M|_1 = max_{all columns j} \sum_{all rows i} |M_ij| \),

This is the natural matrix norm that is compatible to the l1-norm for vectors, i.e. \( |Mv|_1 \leq |M|_1 |v|_1 \).

Definition at line 1125 of file dense_matrix.h.

{
  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;
}

References libMesh::libmesh_assert().

Referenced by libMesh::FEMSystem::assembly().

left_multiply() [1/2]

template<typename T >

void libMesh::DenseMatrix< T >::left_multiply ( const DenseMatrixBase< T > &  M2 )

finaloverridevirtual

Performs the operation: (*this) <- M2 * (*this)

Implements libMesh::DenseMatrixBase< T >.

Definition at line 45 of file dense_matrix_impl.h.

{
  if (this->use_blas_lapack)
    this->_multiply_blas(M2, LEFT_MULTIPLY);
  else
    {
      // (*this) <- M2 * (*this)
      // Where:
      // (*this) = (m x n),
      // M2      = (m x p),
      // M3      = (p x n)
 
      // M3 is a copy of *this before it gets resize()d
      DenseMatrix<T> M3(*this);
 
      // Resize *this so that the result can fit
      this->resize (M2.m(), M3.n());
 
      // Call the multiply function in the base class
      this->multiply(*this, M2, M3);
    }
}

References libMesh::DenseMatrixBase< T >::m(), and libMesh::DenseMatrixBase< T >::n().

left_multiply() [2/2]

template<typename T >

template<typename T2 >

void libMesh::DenseMatrix< T >::left_multiply

(

const DenseMatrixBase< T2 > &

M2

)

Left multiplies by the matrix M2 of different type.

Definition at line 72 of file dense_matrix_impl.h.

{
  // (*this) <- M2 * (*this)
  // Where:
  // (*this) = (m x n),
  // M2      = (m x p),
  // M3      = (p x n)
 
  // M3 is a copy of *this before it gets resize()d
  DenseMatrix<T> M3(*this);
 
  // Resize *this so that the result can fit
  this->resize (M2.m(), M3.n());
 
  // Call the multiply function in the base class
  this->multiply(*this, M2, M3);
}

References libMesh::DenseMatrixBase< T >::m(), and libMesh::DenseMatrixBase< T >::n().

left_multiply_transpose() [1/2]

template<typename T >

void libMesh::DenseMatrix< T >::left_multiply_transpose

(

const DenseMatrix< T > &

A

)

Left multiplies by the transpose of the matrix A.

Definition at line 93 of file dense_matrix_impl.h.

{
  if (this->use_blas_lapack)
    this->_multiply_blas(A, LEFT_MULTIPLY_TRANSPOSE);
  else
    this->_left_multiply_transpose(A);
}

Referenced by libMesh::DofMap::heterogeneously_constrain_element_jacobian_and_residual().

left_multiply_transpose() [2/2]

template<typename T >

template<typename T2 >

void libMesh::DenseMatrix< T >::left_multiply_transpose

(

const DenseMatrix< T2 > &

A

)

Left multiplies by the transpose of the matrix A which contains a different numerical type.

Definition at line 104 of file dense_matrix_impl.h.

{
  this->_left_multiply_transpose(A);
}

linfty_norm()

template<typename T >

auto libMesh::DenseMatrix< T >::linfty_norm ( ) const

inline

Returns

The linfty-norm of the matrix, that is, the max row sum:

\( |M|_\infty = max_{all rows i} \sum_{all columns j} |M_ij| \),

This is the natural matrix norm that is compatible to the linfty-norm of vectors, i.e. \( |Mv|_\infty \leq |M|_\infty |v|_\infty \).

Definition at line 1153 of file dense_matrix.h.

{
  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;
}

References libMesh::libmesh_assert().

lu_solve()

template<typename T >

void libMesh::DenseMatrix< T >::lu_solve	(	const DenseVector< T > &	b,
		DenseVector< T > &	x
	)

Solve the system Ax=b given the input vector b.

Partial pivoting is performed by default in order to keep the algorithm stable to the effects of round-off error.

Important note: once you call lu_solve(), you must not modify the entries of the matrix via calls to operator(i,j) and call lu_solve() again without first calling either zero() or resize(), otherwise the code will skip computing the decomposition of the matrix and go directly to the back substitution step. This is done on purpose for efficiency, so that the same LU decomposition can be used with multiple right-hand sides, but it does also make it possible to "shoot yourself in the foot", so be careful!

Definition at line 521 of file dense_matrix_impl.h.

{
  // Check to be sure that the matrix is square before attempting
  // an LU-solve.  In general, one can compute the LU factorization of
  // a non-square matrix, but:
  //
  // Overdetermined systems (m>n) have a solution only if enough of
  // the equations are linearly-dependent.
  //
  // Underdetermined systems (m<n) typically have infinitely many
  // solutions.
  //
  // We don't want to deal with either of these ambiguous cases here...
  libmesh_assert_equal_to (this->m(), this->n());
 
  switch(this->_decomposition_type)
    {
    case NONE:
      {
        if (this->use_blas_lapack)
          this->_lu_decompose_lapack();
        else
          this->_lu_decompose ();
        break;
      }
 
    case LU_BLAS_LAPACK:
      {
        // Already factored, just need to call back_substitute.
        if (this->use_blas_lapack)
          break;
      }
      libmesh_fallthrough();
 
    case LU:
      {
        // Already factored, just need to call back_substitute.
        if (!(this->use_blas_lapack))
          break;
      }
      libmesh_fallthrough();
 
    default:
      libmesh_error_msg("Error! This matrix already has a different decomposition...");
    }
 
  if (this->use_blas_lapack)
    this->_lu_back_substitute_lapack (b, x);
  else
    this->_lu_back_substitute (b, x);
}

Referenced by compute_enriched_soln(), libMesh::RBEIMConstruction::compute_max_eim_error(), fe_assembly(), libMesh::PatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::SmoothnessEstimator::EstimateSmoothness::operator()(), libMesh::WeightedPatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::RBEIMEvaluation::rb_eim_solve(), libMesh::RBEIMEvaluation::rb_eim_solves(), libMesh::TransientRBEvaluation::rb_solve(), libMesh::RBEvaluation::rb_solve(), libMesh::TransientRBEvaluation::rb_solve_again(), libMesh::TransientRBConstruction::set_error_temporal_data(), and libMesh::TransientRBConstruction::update_RB_initial_condition_all_N().

m()

template<typename T >

unsigned int libMesh::DenseMatrixBase< T >::m ( ) const

inlineinherited

Returns

The row-dimension of the matrix.

Definition at line 104 of file dense_matrix_base.h.

{ return _m; }

References libMesh::DenseMatrixBase< T >::_m.

Referenced by libMesh::DenseMatrix< T >::_left_multiply_transpose(), libMesh::DenseMatrix< T >::_multiply_blas(), libMesh::DenseMatrix< T >::_right_multiply_transpose(), libMesh::DenseMatrix< T >::_svd_solve_lapack(), libMesh::DenseMatrix< T >::add(), libMesh::DenseMatrixBase< T >::add(), libMesh::PetscMatrix< T >::add_block_matrix(), libMesh::StaticCondensation::add_matrix(), libMesh::DiagonalMatrix< T >::add_matrix(), libMesh::EigenSparseMatrix< T >::add_matrix(), libMesh::LaspackMatrix< T >::add_matrix(), libMesh::PetscMatrix< T >::add_matrix(), libMesh::EpetraMatrix< T >::add_matrix(), libMesh::RBConstruction::add_scaled_matrix_and_vector(), libMesh::DofMap::build_constraint_matrix(), libMesh::DofMap::build_constraint_matrix_and_vector(), libMesh::DofMap::extract_local_vector(), libMesh::DenseMatrix< T >::get_transpose(), libMesh::DofMap::heterogeneously_constrain_element_jacobian_and_residual(), libMesh::DofMap::heterogeneously_constrain_element_residual(), libMesh::isfinite(), libMesh::isfinite(), libMesh::isinf(), libMesh::isinf(), libMesh::isnan(), libMesh::isnan(), libMesh::DenseMatrix< T >::left_multiply(), libMesh::DenseMatrix< T >::left_multiply(), libMesh::DofMap::max_constraint_error(), libMesh::DenseMatrixBase< T >::multiply(), libMesh::PatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::WeightedPatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::DenseMatrix< T >::operator=(), libMesh::RBEIMEvaluation::set_interpolation_matrix_entry(), DualShapeTest::testEdge2Lagrange(), DenseMatrixTest::testEVD_helper(), DenseMatrixTest::testSVD(), and MetaPhysicL::RawType< libMesh::DenseMatrix< T > >::value().

max()

template<typename T >

auto libMesh::DenseMatrix< T >::max ( ) const -> decltype(libmesh_real(T(0)))

inline

Returns

The maximum entry in the matrix, or the maximum real part in the case of complex numbers.

Definition at line 1104 of file dense_matrix.h.

{
  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;
}

References libMesh::libmesh_assert(), and libMesh::libmesh_real().

min()

template<typename T >

auto libMesh::DenseMatrix< T >::min ( ) const -> decltype(libmesh_real(T(0)))

inline

Returns

The minimum entry in the matrix, or the minimum real part in the case of complex numbers.

Definition at line 1083 of file dense_matrix.h.

{
  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;
}

References libMesh::libmesh_assert(), and libMesh::libmesh_real().

multiply()

template<typename T >

void libMesh::DenseMatrixBase< T >::multiply ( DenseMatrixBase< T > &  M1, const DenseMatrixBase< T > &  M2, const DenseMatrixBase< T > &  M3  )

staticprotectedinherited

Helper function - Performs the computation M1 = M2 * M3 where: M1 = (m x n) M2 = (m x p) M3 = (p x n)

Definition at line 46 of file dense_matrix_base_impl.h.

  {
    // Assertions to make sure we have been
    // passed matrices of the correct dimension.
    libmesh_assert_equal_to (M1.m(), M2.m());
    libmesh_assert_equal_to (M1.n(), M3.n());
    libmesh_assert_equal_to (M2.n(), M3.m());
 
    const unsigned int m_s = M2.m();
    const unsigned int p_s = M2.n();
    const unsigned int n_s = M1.n();
 
    // Do it this way because there is a
    // decent chance (at least for constraint matrices)
    // that M3(k,j) = 0. when right-multiplying.
    for (unsigned int k=0; k<p_s; k++)
      for (unsigned int j=0; j<n_s; j++)
        if (M3.el(k,j) != 0.)
          for (unsigned int i=0; i<m_s; i++)
            M1.el(i,j) += M2.el(i,k) * M3.el(k,j);
  }

References libMesh::DenseMatrixBase< T >::el(), libMesh::DenseMatrixBase< T >::m(), and libMesh::DenseMatrixBase< T >::n().

n()

template<typename T >

unsigned int libMesh::DenseMatrixBase< T >::n ( ) const

inlineinherited

Returns

The column-dimension of the matrix.

Definition at line 109 of file dense_matrix_base.h.

{ return _n; }

References libMesh::DenseMatrixBase< T >::_n.

Referenced by libMesh::DenseMatrix< T >::_left_multiply_transpose(), libMesh::DenseMatrix< T >::_multiply_blas(), libMesh::DenseMatrix< T >::_right_multiply_transpose(), libMesh::DenseMatrix< T >::_svd_solve_lapack(), libMesh::DenseMatrix< T >::add(), libMesh::DenseMatrixBase< T >::add(), libMesh::PetscMatrix< T >::add_block_matrix(), libMesh::StaticCondensation::add_matrix(), libMesh::DiagonalMatrix< T >::add_matrix(), libMesh::EigenSparseMatrix< T >::add_matrix(), libMesh::LaspackMatrix< T >::add_matrix(), libMesh::PetscMatrix< T >::add_matrix(), libMesh::EpetraMatrix< T >::add_matrix(), libMesh::RBConstruction::add_scaled_matrix_and_vector(), libMesh::DofMap::build_constraint_matrix(), libMesh::DofMap::build_constraint_matrix_and_vector(), libMesh::DofMap::constrain_element_residual(), libMesh::DofMap::extract_local_vector(), libMesh::DenseMatrix< T >::get_transpose(), libMesh::DofMap::heterogeneously_constrain_element_jacobian_and_residual(), libMesh::DofMap::heterogeneously_constrain_element_residual(), libMesh::isfinite(), libMesh::isfinite(), libMesh::isinf(), libMesh::isinf(), libMesh::isnan(), libMesh::isnan(), libMesh::DenseMatrix< T >::left_multiply(), libMesh::DenseMatrix< T >::left_multiply(), libMesh::DofMap::max_constraint_error(), libMesh::DenseMatrixBase< T >::multiply(), libMesh::PatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::WeightedPatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::DenseMatrix< T >::operator=(), libMesh::DenseMatrix< Real >::right_multiply(), libMesh::RBEIMEvaluation::set_interpolation_matrix_entry(), DualShapeTest::testEdge2Lagrange(), DenseMatrixTest::testSVD(), and MetaPhysicL::RawType< libMesh::DenseMatrix< T > >::value().

operator!=()

template<typename T >

bool libMesh::DenseMatrix< T >::operator!= ( const DenseMatrix< T > &  mat ) const

inline

Returns

true if mat is not exactly equal to this matrix, false otherwise.

Definition at line 1046 of file dense_matrix.h.

{
  for (auto i : index_range(_val))
    if (_val[i] != mat._val[i])
      return true;
 
  return false;
}

References libMesh::DenseMatrix< T >::_val, and libMesh::index_range().

operator()() [1/2]

template<typename T >

T & libMesh::DenseMatrix< T >::operator() ( const unsigned int  i, const unsigned int  j  )

inline

Returns

The (i,j) element of the matrix as a writable reference.

Definition at line 969 of file dense_matrix.h.

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

operator()() [2/2]

template<typename T >

T libMesh::DenseMatrix< T >::operator() ( const unsigned int  i, const unsigned int  j  ) const

inline

Returns

The (i,j) element of the matrix.

Definition at line 953 of file dense_matrix.h.

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

operator*=()

template<typename T >

DenseMatrix< T > & libMesh::DenseMatrix< T >::operator*= ( const T  factor )

inline

Multiplies every element in the matrix by factor.

Returns

A reference to *this.

Definition at line 1005 of file dense_matrix.h.

{
  this->scale(factor);
  return *this;
}

operator+=()

template<typename T >

DenseMatrix< T > & libMesh::DenseMatrix< T >::operator+= ( const DenseMatrix< T > &  mat )

inline

Adds mat to this matrix.

Returns

A reference to *this.

Definition at line 1059 of file dense_matrix.h.

{
  for (auto i : index_range(_val))
    _val[i] += mat._val[i];
 
  return *this;
}

References libMesh::DenseMatrix< T >::_val, and libMesh::index_range().

operator-=()

template<typename T >

DenseMatrix< T > & libMesh::DenseMatrix< T >::operator-= ( const DenseMatrix< T > &  mat )

inline

Subtracts mat from this matrix.

Returns

A reference to *this.

Definition at line 1071 of file dense_matrix.h.

{
  for (auto i : index_range(_val))
    _val[i] -= mat._val[i];
 
  return *this;
}

References libMesh::DenseMatrix< T >::_val, and libMesh::index_range().

operator=() [1/3]

template<typename T >

DenseMatrix & libMesh::DenseMatrix< T >::operator= ( const DenseMatrix< T > &  )

default

operator=() [2/3]

template<typename T >

template<typename T2 >

DenseMatrix< T > & libMesh::DenseMatrix< T >::operator= ( const DenseMatrix< T2 > &  other_matrix )

inline

Assignment-from-other-matrix-type operator.

Copies the dense matrix of type T2 into the present matrix. This is useful for copying real matrices into complex ones for further operations.

Returns

A reference to *this.

Definition at line 880 of file dense_matrix.h.

{
  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;
}

References libMesh::DenseMatrixBase< T >::m(), and libMesh::DenseMatrixBase< T >::n().

operator=() [3/3]

template<typename T >

DenseMatrix & libMesh::DenseMatrix< T >::operator= ( DenseMatrix< T > &&  )

default

operator==()

template<typename T >

bool libMesh::DenseMatrix< T >::operator== ( const DenseMatrix< T > &  mat ) const

inline

Returns

true if mat is exactly equal to this matrix, false otherwise.

Definition at line 1033 of file dense_matrix.h.

{
  for (auto i : index_range(_val))
    if (_val[i] != mat._val[i])
      return false;
 
  return true;
}

References libMesh::DenseMatrix< T >::_val, and libMesh::index_range().

outer_product()

template<typename T >

void libMesh::DenseMatrix< T >::outer_product	(	const DenseVector< T > &	a,
		const DenseVector< T > &	b
	)

Computes the outer (dyadic) product of two vectors and stores in (*this).

The outer product of two real-valued vectors \(\mathbf{a}\) and \(\mathbf{b}\) is

\[ (\mathbf{a}\mathbf{b}^T)_{i,j} = \mathbf{a}_i \mathbf{b}_j . \]

The outer product of two complex-valued vectors \(\mathbf{a}\) and \(\mathbf{b}\) is

\[ (\mathbf{a}\mathbf{b}^H)_{i,j} = \mathbf{a}_i \mathbf{b}^*_j , \]

where \(H\) denotes the conjugate transpose of the vector and \(*\) denotes the complex conjugate.

Parameters

[in]	a	Vector whose entries correspond to rows in the product matrix.
[in]	b	Vector whose entries correspond to columns in the product matrix.

Definition at line 485 of file dense_matrix_impl.h.

{
  const unsigned int m = a.size();
  const unsigned int n = b.size();
 
  this->resize(m, n);
  for (unsigned int i = 0; i < m; ++i)
    for (unsigned int j = 0; j < n; ++j)
      (*this)(i,j) = a(i) * libmesh_conj(b(j));
}

References b, libMesh::libmesh_conj(), and libMesh::DenseVector< T >::size().

Referenced by DenseMatrixTest::testOuterProduct().

print()

template<typename T >

void libMesh::DenseMatrixBase< T >::print ( std::ostream &  os = libMesh::out ) const

inherited

Pretty-print the matrix, by default to libMesh::out.

Definition at line 125 of file dense_matrix_base_impl.h.

  {
    for (auto i : make_range(this->m()))
      {
        for (auto j : make_range(this->n()))
          os << std::setw(8)
             << this->el(i,j) << " ";
 
        os << std::endl;
      }
 
    return;
  }

References libMesh::make_range().

print_scientific()

template<typename T >

void libMesh::DenseMatrixBase< T >::print_scientific ( std::ostream &  os, unsigned  precision = 8  ) const

inherited

Prints the matrix entries with more decimal places in scientific notation.

Definition at line 101 of file dense_matrix_base_impl.h.

  {
    // save the initial format flags
    std::ios_base::fmtflags os_flags = os.flags();
 
    // Print the matrix entries.
    for (auto i : make_range(this->m()))
      {
        for (auto j : make_range(this->n()))
          os << std::setw(15)
             << std::scientific
             << std::setprecision(precision)
             << this->el(i,j) << " ";
 
        os << std::endl;
      }
 
    // reset the original format flags
    os.flags(os_flags);
  }

References libMesh::make_range().

resize()

template<typename T >

void libMesh::DenseMatrix< T >::resize ( const unsigned int  new_m, const unsigned int  new_n  )

inline

Resizes the matrix to the specified size and calls zero().

Will never free memory, but may allocate more. Note: when the matrix is zero()'d, any decomposition (LU, Cholesky, etc.) is also cleared, forcing a new decomposition to be computed the next time e.g. lu_solve() is called.

Definition at line 895 of file dense_matrix.h.

{
  _val.resize(new_m*new_n);
 
  this->_m = new_m;
  this->_n = new_n;
 
  // zero and set decomposition_type to NONE
  this->zero();
}

References libMesh::DenseMatrix< T >::resize(), and libMesh::zero.

Referenced by libMesh::DenseMatrix< T >::_svd_lapack(), libMesh::HPCoarsenTest::add_projection(), alternative_fe_assembly(), LinearElasticity::assemble(), assemble(), assemble(), assemble_1D(), AssembleOptimization::assemble_A_and_F(), assemble_biharmonic(), assemble_cd(), assemble_cd(), assemble_divgrad(), assemble_elasticity(), assemble_ellipticdg(), assemble_func(), assemble_graddiv(), assemble_helmholtz(), libMesh::ClawSystem::assemble_jump_coupling_matrix(), assemble_laplace(), assemble_mass(), libMesh::ClawSystem::assemble_mass_matrix(), assemble_matrices(), assemble_matrix_and_rhs(), assemble_poisson(), assemble_poisson(), assemble_poisson(), assemble_SchroedingerEquation(), assemble_shell(), assemble_shell(), assemble_stokes(), assemble_temperature_jump(), assemble_wave(), assemble_wave(), libMesh::DofMap::build_constraint_matrix(), libMesh::DofMap::build_constraint_matrix_and_vector(), libMesh::TransientRBEvaluation::cache_online_residual_terms(), NonlinearNeoHookeCurrentConfig::calculate_tangent(), libMesh::RBEIMConstruction::clear(), libMesh::RBEIMEvaluation::clear(), libMesh::StaticCondensation::close(), libMesh::FEGenericBase< OutputType >::coarsened_dof_values(), compute_enriched_soln(), compute_jacobian(), libMesh::FEGenericBase< OutputType >::compute_periodic_constraints(), libMesh::FEGenericBase< OutputType >::compute_proj_constraints(), libMesh::HDGProblem::compute_stress(), compute_stresses(), libMesh::DenseMatrix< T >::DenseMatrix(), fe_assembly(), form_matrixA(), NonlinearNeoHookeCurrentConfig::get_linearized_stiffness(), libMesh::DenseMatrix< T >::get_principal_submatrix(), NonlinearNeoHookeCurrentConfig::get_residual(), libMesh::DenseMatrix< T >::get_transpose(), libMesh::FESubdivision::init_subdivision_matrix(), LargeDeformationElasticity::jacobian(), libMesh::HDGProblem::jacobian(), LaplaceYoung::jacobian(), libMesh::DGFEMContext::neighbor_side_fe_reinit(), libMesh::BoundaryProjectSolution::operator()(), periodic_bc_test_poisson(), libMesh::FEMContext::pre_fe_reinit(), libMesh::TransientRBEvaluation::rb_solve(), libMesh::RBEIMConstruction::reinit_eim_projection_matrix(), LinearElasticityWithContact::residual_and_jacobian(), Biharmonic::JR::residual_and_jacobian(), libMesh::DenseMatrix< T >::resize(), libMesh::RBEIMEvaluation::resize_data_structures(), libMesh::TransientRBEvaluation::resize_data_structures(), OverlappingCouplingGhostingTest::run_sparsity_pattern_test(), libMesh::HPCoarsenTest::select_refinement(), libMesh::StaticCondensation::StaticCondensation(), libMesh::RBEIMConstruction::train_eim_approximation_with_greedy(), and libMesh::RBEIMConstruction::train_eim_approximation_with_POD().

right_multiply() [1/2]

template<typename T >

void libMesh::DenseMatrix< T >::right_multiply ( const DenseMatrixBase< T > &  M3 )

finaloverridevirtual

Performs the operation: (*this) <- (*this) * M3.

Implements libMesh::DenseMatrixBase< T >.

Definition at line 171 of file dense_matrix_impl.h.

{
  if (this->use_blas_lapack)
    this->_multiply_blas(M3, RIGHT_MULTIPLY);
  else
    {
      // (*this) <- (*this) * M3
      // Where:
      // (*this) = (m x n),
      // M2      = (m x p),
      // M3      = (p x n)
 
      // M2 is a copy of *this before it gets resize()d
      DenseMatrix<T> M2(*this);
 
      // Resize *this so that the result can fit
      this->resize (M2.m(), M3.n());
 
      this->multiply(*this, M2, M3);
    }
}

Referenced by libMesh::DofMap::build_constraint_matrix(), libMesh::DofMap::build_constraint_matrix_and_vector(), NonlinearNeoHookeCurrentConfig::get_linearized_stiffness(), libMesh::DofMap::heterogeneously_constrain_element_jacobian_and_residual(), and libMesh::FESubdivision::init_shape_functions().

right_multiply() [2/2]

template<typename T >

template<typename T2 >

void libMesh::DenseMatrix< T >::right_multiply

(

const DenseMatrixBase< T2 > &

M2

)

Right multiplies by the matrix M2 of different type.

Definition at line 197 of file dense_matrix_impl.h.

{
  // (*this) <- M3 * (*this)
  // Where:
  // (*this) = (m x n),
  // M2      = (m x p),
  // M3      = (p x n)
 
  // M2 is a copy of *this before it gets resize()d
  DenseMatrix<T> M2(*this);
 
  // Resize *this so that the result can fit
  this->resize (M2.m(), M3.n());
 
  this->multiply(*this, M2, M3);
}

right_multiply_transpose() [1/2]

template<typename T >

void libMesh::DenseMatrix< T >::right_multiply_transpose

(

const DenseMatrix< T > &

A

)

Right multiplies by the transpose of the matrix A.

Definition at line 218 of file dense_matrix_impl.h.

{
  if (this->use_blas_lapack)
    this->_multiply_blas(B, RIGHT_MULTIPLY_TRANSPOSE);
  else
    this->_right_multiply_transpose(B);
}

Referenced by NonlinearNeoHookeCurrentConfig::get_linearized_stiffness().

right_multiply_transpose() [2/2]

template<typename T >

template<typename T2 >

void libMesh::DenseMatrix< T >::right_multiply_transpose

(

const DenseMatrix< T2 > &

A

)

Right multiplies by the transpose of the matrix A which contains a different numerical type.

Definition at line 230 of file dense_matrix_impl.h.

{
  this->_right_multiply_transpose(B);
}

scale()

template<typename T >

void libMesh::DenseMatrix< T >::scale ( const T  factor )

inline

Multiplies every element in the matrix by factor.

Definition at line 986 of file dense_matrix.h.

{
  for (auto & v : _val)
    v *= factor;
}

Referenced by compute_stresses(), and SolidSystem::element_time_derivative().

scale_column()

template<typename T >

void libMesh::DenseMatrix< T >::scale_column ( const unsigned int  col, const T  factor  )

inline

Multiplies every element in the column col matrix by factor.

Definition at line 995 of file dense_matrix.h.

{
  for (auto i : make_range(this->m()))
    (*this)(i, col) *= factor;
}

References libMesh::make_range().

sub_matrix()

template<typename T >

DenseMatrix< T > libMesh::DenseMatrix< T >::sub_matrix ( unsigned int  row_id, unsigned int  row_size, unsigned int  col_id, unsigned int  col_size  ) const

inline

Get submatrix with the smallest row and column indices and the submatrix size.

Definition at line 922 of file dense_matrix.h.

{
  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;
}

References libMesh::DenseMatrixBase< T >::_m, libMesh::DenseMatrixBase< T >::_n, and libMesh::DenseMatrix< T >::_val.

Referenced by DenseMatrixTest::testSubMatrix().

svd() [1/2]

template<typename T >

void libMesh::DenseMatrix< T >::svd

(

DenseVector< Real > &

sigma

)

Compute the singular value decomposition of the matrix.

On exit, sigma holds all of the singular values (in descending order).

The implementation uses PETSc's interface to BLAS/LAPACK. If this is not available, this function throws an error.

Definition at line 710 of file dense_matrix_impl.h.

{
  // We use the LAPACK svd implementation
  _svd_lapack(sigma);
}

Referenced by libMesh::TransientRBConstruction::enrich_RB_space(), DenseMatrixTest::testComplexSVD(), DenseMatrixTest::testSVD(), libMesh::RBEIMConstruction::train_eim_approximation_with_POD(), and libMesh::RBConstruction::train_reduced_basis_with_POD().

svd() [2/2]

template<typename T >

void libMesh::DenseMatrix< T >::svd	(	DenseVector< Real > &	sigma,
		DenseMatrix< Number > &	U,
		DenseMatrix< Number > &	VT
	)

Compute the "reduced" singular value decomposition of the matrix.

On exit, sigma holds all of the singular values (in descending order), U holds the left singular vectors, and VT holds the transpose of the right singular vectors. In the reduced SVD, U has min(m,n) columns and VT has min(m,n) rows. (In the "full" SVD, U and VT would be square.)

The implementation uses PETSc's interface to BLAS/LAPACK. If this is not available, this function throws an error.

Definition at line 718 of file dense_matrix_impl.h.

{
  // We use the LAPACK svd implementation
  _svd_lapack(sigma, U, VT);
}

svd_solve()

template<typename T >

void libMesh::DenseMatrix< T >::svd_solve	(	const DenseVector< T > &	rhs,
		DenseVector< T > &	x,
		Real	rcond = std::numeric_limits<Real>::epsilon()
	)		const

Solve the system of equations \( A x = rhs \) for \( x \) in the least-squares sense.

\( A \) may be non-square and/or rank-deficient. You can control which singular values are treated as zero by changing the "rcond" parameter. Singular values S(i) for which S(i) <= rcond*S(1) are treated as zero for purposes of the solve. Passing a negative number for rcond forces a "machine precision" value to be used instead.

This function is marked const, since due to various implementation details, we do not need to modify the contents of A in order to compute the SVD (a copy is made internally instead).

Requires PETSc >= 3.1 since this was the first version to provide the LAPACKgelss_ wrapper.

Definition at line 729 of file dense_matrix_impl.h.

{
  _svd_solve_lapack(rhs, x, rcond);
}

Referenced by fe_assembly().

swap()

template<typename T >

void libMesh::DenseMatrix< T >::swap ( DenseMatrix< T > &  other_matrix )

inline

STL-like swap method.

Definition at line 865 of file dense_matrix.h.

{
  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;
}

References libMesh::DenseMatrix< T >::_decomposition_type, libMesh::DenseMatrixBase< T >::_m, libMesh::DenseMatrixBase< T >::_n, and libMesh::DenseMatrix< T >::_val.

Referenced by libMesh::NewmarkSolver::_general_residual(), libMesh::Euler2Solver::_general_residual(), libMesh::EulerSolver::_general_residual(), and libMesh::FEMSystem::assembly().

transpose()

template<typename T >

T libMesh::DenseMatrix< T >::transpose ( const unsigned int  i, const unsigned int  j  ) const

inline

Returns

The (i,j) element of the transposed matrix.

Definition at line 1181 of file dense_matrix.h.

{
  // Implement in terms of operator()
  return (*this)(j,i);
}

vector_mult() [1/2]

template<typename T >

void libMesh::DenseMatrix< T >::vector_mult	(	DenseVector< T > &	dest,
		const DenseVector< T > &	arg
	)		const

Performs the matrix-vector multiplication, dest := (*this) * arg.

Definition at line 299 of file dense_matrix_impl.h.

{
  // Make sure the input sizes are compatible
  libmesh_assert_equal_to (this->n(), arg.size());
 
  // Resize and clear dest.
  // Note: DenseVector::resize() also zeros the vector.
  dest.resize(this->m());
 
  // Short-circuit if the matrix is empty
  if(this->m() == 0 || this->n() == 0)
    return;
 
  if (this->use_blas_lapack)
    this->_matvec_blas(1., 0., dest, arg);
  else
    {
      const unsigned int n_rows = this->m();
      const unsigned int n_cols = this->n();
 
      for (unsigned int i=0; i<n_rows; i++)
        for (unsigned int j=0; j<n_cols; j++)
          dest(i) += (*this)(i,j)*arg(j);
    }
}

References libMesh::DenseVector< T >::resize(), and libMesh::DenseVector< T >::size().

Referenced by libMesh::FEMap::compute_single_point_map(), NonlinearNeoHookeCurrentConfig::get_residual(), libMesh::TransientRBEvaluation::rb_solve(), and libMesh::TransientRBEvaluation::rb_solve_again().

vector_mult() [2/2]

template<typename T >

template<typename T2 >

void libMesh::DenseMatrix< T >::vector_mult	(	DenseVector< typename CompareTypes< T, T2 >::supertype > &	dest,
		const DenseVector< T2 > &	arg
	)		const

Performs the matrix-vector multiplication, dest := (*this) * arg on mixed types.

Definition at line 330 of file dense_matrix_impl.h.

{
  // Make sure the input sizes are compatible
  libmesh_assert_equal_to (this->n(), arg.size());
 
  // Resize and clear dest.
  // Note: DenseVector::resize() also zeros the vector.
  dest.resize(this->m());
 
  // Short-circuit if the matrix is empty
  if (this->m() == 0 || this->n() == 0)
    return;
 
  const unsigned int n_rows = this->m();
  const unsigned int n_cols = this->n();
 
  for (unsigned int i=0; i<n_rows; i++)
    for (unsigned int j=0; j<n_cols; j++)
      dest(i) += (*this)(i,j)*arg(j);
}

References libMesh::DenseVector< T >::size().

vector_mult_add() [1/2]

template<typename T >

void libMesh::DenseMatrix< T >::vector_mult_add	(	DenseVector< T > &	dest,
		const T	factor,
		const DenseVector< T > &	arg
	)		const

Performs the scaled matrix-vector multiplication, dest += factor * (*this) * arg.

Definition at line 425 of file dense_matrix_impl.h.

{
  // Short-circuit if the matrix is empty
  if (this->m() == 0)
    {
      dest.resize(0);
      return;
    }
 
  if (this->use_blas_lapack)
    this->_matvec_blas(factor, 1., dest, arg);
  else
    {
      DenseVector<T> temp(arg.size());
      this->vector_mult(temp, arg);
      dest.add(factor, temp);
    }
}

References libMesh::DenseVector< T >::resize().

Referenced by libMesh::DofMap::build_constraint_matrix_and_vector().

vector_mult_add() [2/2]

template<typename T >

template<typename T2 , typename T3 >

void libMesh::DenseMatrix< T >::vector_mult_add	(	DenseVector< typename CompareTypes< T, typename CompareTypes< T2, T3 >::supertype >::supertype > &	dest,
		const T2	factor,
		const DenseVector< T3 > &	arg
	)		const

Performs the scaled matrix-vector multiplication, dest += factor * (*this) * arg.

on mixed types

Definition at line 450 of file dense_matrix_impl.h.

{
  // Short-circuit if the matrix is empty
  if (this->m() == 0)
    {
      dest.resize(0);
      return;
    }
 
  DenseVector<typename CompareTypes<T,T3>::supertype>
    temp(arg.size());
  this->vector_mult(temp, arg);
  dest.add(factor, temp);
}

vector_mult_transpose() [1/2]

template<typename T >

void libMesh::DenseMatrix< T >::vector_mult_transpose	(	DenseVector< T > &	dest,
		const DenseVector< T > &	arg
	)		const

Performs the matrix-vector multiplication, dest := (*this)^T * arg.

Definition at line 355 of file dense_matrix_impl.h.

{
  // Make sure the input sizes are compatible
  libmesh_assert_equal_to (this->m(), arg.size());
 
  // Resize and clear dest.
  // Note: DenseVector::resize() also zeros the vector.
  dest.resize(this->n());
 
  // Short-circuit if the matrix is empty
  if (this->m() == 0)
    return;
 
  if (this->use_blas_lapack)
    {
      this->_matvec_blas(1., 0., dest, arg, /*trans=*/true);
    }
  else
    {
      const unsigned int n_rows = this->m();
      const unsigned int n_cols = this->n();
 
      // WORKS
      // for (unsigned int j=0; j<n_cols; j++)
      //   for (unsigned int i=0; i<n_rows; i++)
      //     dest(j) += (*this)(i,j)*arg(i);
 
      // ALSO WORKS, (i,j) just swapped
      for (unsigned int i=0; i<n_cols; i++)
        for (unsigned int j=0; j<n_rows; j++)
          dest(i) += (*this)(j,i)*arg(j);
    }
}

Referenced by libMesh::DofMap::constrain_element_residual(), libMesh::DofMap::heterogeneously_constrain_element_jacobian_and_residual(), and libMesh::DofMap::heterogeneously_constrain_element_residual().

vector_mult_transpose() [2/2]

template<typename T >

template<typename T2 >

void libMesh::DenseMatrix< T >::vector_mult_transpose	(	DenseVector< typename CompareTypes< T, T2 >::supertype > &	dest,
		const DenseVector< T2 > &	arg
	)		const

Performs the matrix-vector multiplication, dest := (*this)^T * arg.

on mixed types

Definition at line 394 of file dense_matrix_impl.h.

{
  // Make sure the input sizes are compatible
  libmesh_assert_equal_to (this->m(), arg.size());
 
  // Resize and clear dest.
  // Note: DenseVector::resize() also zeros the vector.
  dest.resize(this->n());
 
  // Short-circuit if the matrix is empty
  if (this->m() == 0)
    return;
 
  const unsigned int n_rows = this->m();
  const unsigned int n_cols = this->n();
 
  // WORKS
  // for (unsigned int j=0; j<n_cols; j++)
  //   for (unsigned int i=0; i<n_rows; i++)
  //     dest(j) += (*this)(i,j)*arg(i);
 
  // ALSO WORKS, (i,j) just swapped
  for (unsigned int i=0; i<n_cols; i++)
    for (unsigned int j=0; j<n_rows; j++)
      dest(i) += (*this)(j,i)*arg(j);
}

References libMesh::DenseVector< T >::resize(), and libMesh::DenseVector< T >::size().

zero()

template<typename T >

void libMesh::DenseMatrix< T >::zero ( )

inlinefinaloverridevirtual

Sets all elements of the matrix to 0 and resets any decomposition flag which may have been previously set.

This allows e.g. a new LU decomposition to be computed while reusing the same storage.

Implements libMesh::DenseMatrixBase< T >.

Definition at line 911 of file dense_matrix.h.

{
  _decomposition_type = NONE;
 
  std::fill (_val.begin(), _val.end(), static_cast<T>(0));
}

Referenced by libMesh::HPCoarsenTest::add_projection(), assemble_temperature_jump(), libMesh::FEMSystem::assembly(), libMesh::FEGenericBase< OutputType >::coarsened_dof_values(), libMesh::BoundaryProjectSolution::operator()(), libMesh::TransientRBEvaluation::rb_solve(), libMesh::RBEvaluation::rb_solve(), libMesh::HPCoarsenTest::select_refinement(), SparseMatrixTest< DerivedClass >::setValues(), DiagonalMatrixTest::testNumerics(), and libMesh::TransientRBConstruction::update_RB_system_matrices().

Member Data Documentation

_decomposition_type

template<typename T >

DecompositionType libMesh::DenseMatrix< T >::_decomposition_type

private

This flag keeps track of which type of decomposition has been performed on the matrix.

Definition at line 650 of file dense_matrix.h.

Referenced by libMesh::DenseMatrix< T >::swap().

_m

template<typename T >

unsigned int libMesh::DenseMatrixBase< T >::_m

protectedinherited

The row dimension.

Definition at line 176 of file dense_matrix_base.h.

Referenced by libMesh::DenseMatrixBase< T >::m(), libMesh::DenseMatrix< T >::sub_matrix(), and libMesh::DenseMatrix< T >::swap().

_n

template<typename T >

unsigned int libMesh::DenseMatrixBase< T >::_n

protectedinherited

The column dimension.

Definition at line 181 of file dense_matrix_base.h.

Referenced by libMesh::DenseMatrixBase< T >::n(), libMesh::DenseMatrix< T >::sub_matrix(), and libMesh::DenseMatrix< T >::swap().

_pivots

template<typename T >

std::vector<pivot_index_t> libMesh::DenseMatrix< T >::_pivots

private

Definition at line 740 of file dense_matrix.h.

_val

template<typename T >

std::vector<T> libMesh::DenseMatrix< T >::_val

private

The actual data values, stored as a 1D array.

Definition at line 600 of file dense_matrix.h.

Referenced by libMesh::DenseMatrix< T >::get_values(), libMesh::DenseMatrix< T >::get_values(), libMesh::DenseMatrix< T >::operator!=(), libMesh::DenseMatrix< T >::operator+=(), libMesh::DenseMatrix< T >::operator-=(), libMesh::DenseMatrix< T >::operator==(), libMesh::DenseMatrix< T >::sub_matrix(), and libMesh::DenseMatrix< T >::swap().

use_blas_lapack

template<typename T >

bool libMesh::DenseMatrix< T >::use_blas_lapack

Computes the inverse of the dense matrix (assuming it is invertible) by first computing the LU decomposition and then performing multiple back substitution steps.

Follows the algorithm from Numerical Recipes in C that is available on the web.

This routine is commented out since it is not really a memory- or computationally- efficient implementation. Also, you typically don't need the actual inverse for anything, and can use something like lu_solve() instead. Run-time selectable option to turn on/off BLAS support. This was primarily used for testing purposes, and could be removed...

Definition at line 585 of file dense_matrix.h.

---

The documentation for this class was generated from the following files:

• include/base/dof_map.h
• include/numerics/dense_matrix.h
• include/numerics/dense_matrix_impl.h
• src/numerics/dense_matrix.C
• src/numerics/dense_matrix_blas_lapack.C
• src/systems/system_projection.C