libMesh
Classes | Public Member Functions | Public Attributes | Protected Member Functions | Static Protected Member Functions | Protected Attributes | Private Types | Private Member Functions | Private Attributes | List of all members
libMesh::DenseMatrix< T > Class Template Reference

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

#include <dof_map.h>

Inheritance diagram for libMesh::DenseMatrix< T >:
[legend]

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. More...
 
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. More...
 
 DenseMatrix (DenseMatrix &&)=default
 The 5 special functions can be defaulted for this class, as it does not manage any memory itself. More...
 
 DenseMatrix (const DenseMatrix &)=default
 
DenseMatrixoperator= (const DenseMatrix &)=default
 
DenseMatrixoperator= (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. More...
 
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. More...
 
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) More...
 
template<typename T2 >
void left_multiply (const DenseMatrixBase< T2 > &M2)
 Left multiplies by the matrix M2 of different type. More...
 
virtual void right_multiply (const DenseMatrixBase< T > &M2) override final
 Performs the operation: (*this) <- (*this) * M3. More...
 
template<typename T2 >
void right_multiply (const DenseMatrixBase< T2 > &M2)
 Right multiplies by the matrix M2 of different type. More...
 
void vector_mult (DenseVector< T > &dest, const DenseVector< T > &arg) const
 Performs the matrix-vector multiplication, dest := (*this) * arg. More...
 
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. More...
 
void vector_mult_transpose (DenseVector< T > &dest, const DenseVector< T > &arg) const
 Performs the matrix-vector multiplication, dest := (*this)^T * arg. More...
 
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. More...
 
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. More...
 
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. More...
 
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. More...
 
void get_principal_submatrix (unsigned int sub_m, DenseMatrix< T > &dest) const
 Put the sub_m x sub_m principal submatrix into dest. More...
 
void outer_product (const DenseVector< T > &a, const DenseVector< T > &b)
 Computes the outer (dyadic) product of two vectors and stores in (*this). More...
 
template<typename T2 >
DenseMatrix< T > & operator= (const DenseMatrix< T2 > &other_matrix)
 Assignment-from-other-matrix-type operator. More...
 
void swap (DenseMatrix< T > &other_matrix)
 STL-like swap method. More...
 
void resize (const unsigned int new_m, const unsigned int new_n)
 Resizes the matrix to the specified size and calls zero(). More...
 
void scale (const T factor)
 Multiplies every element in the matrix by factor. More...
 
void scale_column (const unsigned int col, const T factor)
 Multiplies every element in the column col matrix by factor. More...
 
DenseMatrix< T > & operator*= (const T factor)
 Multiplies every element in the matrix by factor. More...
 
template<typename T2 , typename T3 >
boostcopy::enable_if_c< ScalarTraits< T2 >::value, void >::type add (const T2 factor, const DenseMatrix< T3 > &mat)
 Adds factor times mat to this matrix. More...
 
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. More...
 
DenseMatrix< T > & operator-= (const DenseMatrix< T > &mat)
 Subtracts mat from this matrix. More...
 
auto min () const -> decltype(libmesh_real(T(0)))
 
auto max () const -> decltype(libmesh_real(T(0)))
 
auto l1_norm () const -> decltype(std::abs(T(0)))
 
auto linfty_norm () const -> decltype(std::abs(T(0)))
 
void left_multiply_transpose (const DenseMatrix< T > &A)
 Left multiplies by the transpose of the matrix A. More...
 
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. More...
 
void right_multiply_transpose (const DenseMatrix< T > &A)
 Right multiplies by the transpose of the matrix A. More...
 
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. More...
 
transpose (const unsigned int i, const unsigned int j) const
 
void get_transpose (DenseMatrix< T > &dest) const
 Put the tranposed matrix into dest. More...
 
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. More...
 
void lu_solve (const DenseVector< T > &b, DenseVector< T > &x)
 Solve the system Ax=b given the input vector b. More...
 
template<typename T2 >
void cholesky_solve (const DenseVector< T2 > &b, DenseVector< T2 > &x)
 For symmetric positive definite (SPD) matrices. More...
 
void svd (DenseVector< Real > &sigma)
 Compute the singular value decomposition of the matrix. More...
 
void svd (DenseVector< Real > &sigma, DenseMatrix< Number > &U, DenseMatrix< Number > &VT)
 Compute the "reduced" singular value decomposition of the matrix. More...
 
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. More...
 
void evd (DenseVector< T > &lambda_real, DenseVector< T > &lambda_imag)
 Compute the eigenvalues (both real and imaginary parts) of a general matrix. More...
 
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 \). More...
 
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 \). More...
 
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. More...
 
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. More...
 
void print_scientific (std::ostream &os, unsigned precision=8) const
 Prints the matrix entries with more decimal places in scientific notation. More...
 
template<typename T2 , typename T3 >
boostcopy::enable_if_c< ScalarTraits< T2 >::value, void >::type add (const T2 factor, const DenseMatrixBase< T3 > &mat)
 Adds factor to every element in the matrix. More...
 
DenseVector< T > diagonal () const
 Return the matrix diagonal. More...
 

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. More...
 

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. More...
 

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) More...
 

Protected Attributes

unsigned int _m
 The row dimension. More...
 
unsigned int _n
 The column dimension. More...
 

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. More...
 
typedef int pivot_index_t
 

Private Member Functions

void _lu_decompose ()
 Form the LU decomposition of the matrix. More...
 
void _lu_back_substitute (const DenseVector< T > &b, DenseVector< T > &x) const
 Solves the system Ax=b through back substitution. More...
 
void _cholesky_decompose ()
 Decomposes a symmetric positive definite matrix into a product of two lower triangular matrices according to A = LL^T. More...
 
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. More...
 
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. More...
 
void _lu_decompose_lapack ()
 Computes an LU factorization of the matrix using the Lapack routine "getrf". More...
 
void _svd_lapack (DenseVector< Real > &sigma)
 Computes an SVD of the matrix using the Lapack routine "getsvd". More...
 
void _svd_lapack (DenseVector< Real > &sigma, DenseMatrix< Number > &U, DenseMatrix< Number > &VT)
 Computes a "reduced" SVD of the matrix using the Lapack routine "getsvd". More...
 
void _svd_solve_lapack (const DenseVector< T > &rhs, DenseVector< T > &x, Real rcond) const
 Called by svd_solve(rhs). More...
 
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. More...
 
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". More...
 
void _lu_back_substitute_lapack (const DenseVector< T > &b, DenseVector< T > &x)
 Companion function to _lu_decompose_lapack(). More...
 
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. More...
 
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. More...
 
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. More...
 

Private Attributes

std::vector< T > _val
 The actual data values, stored as a 1D array. More...
 
DecompositionType _decomposition_type
 This flag keeps track of which type of decomposition has been performed on the matrix. More...
 
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 75 of file dof_map.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.

◆ 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.

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.

837  :
838  DenseMatrixBase<T>(new_m,new_n),
840  _val(),
842 {
843  this->resize(new_m,new_n);
844 }
DecompositionType _decomposition_type
This flag keeps track of which type of decomposition has been performed on the matrix.
Definition: dense_matrix.h:650
bool use_blas_lapack
Computes the inverse of the dense matrix (assuming it is invertible) by first computing the LU decomp...
Definition: dense_matrix.h:585
void resize(const unsigned int new_m, const unsigned int new_n)
Resizes the matrix to the specified size and calls zero().
Definition: dense_matrix.h:895
std::vector< T > _val
The actual data values, stored as a 1D array.
Definition: dense_matrix.h:600

◆ 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.

850  :
851  DenseMatrixBase<T>(nrow, ncol),
853  _val(init_list.begin(), init_list.end()),
855 {
856  // Make sure the user passed us an amount of data which is
857  // consistent with the size of the matrix.
858  libmesh_assert_equal_to(nrow * ncol, init_list.size());
859 }
DecompositionType _decomposition_type
This flag keeps track of which type of decomposition has been performed on the matrix.
Definition: dense_matrix.h:650
bool use_blas_lapack
Computes the inverse of the dense matrix (assuming it is invertible) by first computing the LU decomp...
Definition: dense_matrix.h:585
std::vector< T > _val
The actual data values, stored as a 1D array.
Definition: dense_matrix.h:600

◆ 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 919 of file dense_matrix_impl.h.

921 {
922  // Shorthand notation for number of rows and columns.
923  const unsigned int
924  n_rows = this->m(),
925  n_cols = this->n();
926 
927  // Just to be really sure...
928  libmesh_assert_equal_to (n_rows, n_cols);
929 
930  // A convenient reference to *this
931  const DenseMatrix<T> & A = *this;
932 
933  // Now compute the solution to Ax =b using the factorization.
934  x.resize(n_rows);
935 
936  // Solve for Ly=b
937  for (unsigned int i=0; i<n_cols; ++i)
938  {
939  T2 temp = b(i);
940 
941  for (unsigned int k=0; k<i; ++k)
942  temp -= A(i,k)*x(k);
943 
944  x(i) = temp / A(i,i);
945  }
946 
947  // Solve for L^T x = y
948  for (unsigned int i=0; i<n_cols; ++i)
949  {
950  const unsigned int ib = (n_cols-1)-i;
951 
952  for (unsigned int k=(ib+1); k<n_cols; ++k)
953  x(ib) -= A(k,ib) * x(k);
954 
955  x(ib) /= A(ib,ib);
956  }
957 }
unsigned int m() const
unsigned int n() const

◆ _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 873 of file dense_matrix_impl.h.

874 {
875  // If we called this function, there better not be any
876  // previous decomposition of the matrix.
877  libmesh_assert_equal_to (this->_decomposition_type, NONE);
878 
879  // Shorthand notation for number of rows and columns.
880  const unsigned int
881  n_rows = this->m(),
882  n_cols = this->n();
883 
884  // Just to be really sure...
885  libmesh_assert_equal_to (n_rows, n_cols);
886 
887  // A convenient reference to *this
888  DenseMatrix<T> & A = *this;
889 
890  for (unsigned int i=0; i<n_rows; ++i)
891  {
892  for (unsigned int j=i; j<n_cols; ++j)
893  {
894  for (unsigned int k=0; k<i; ++k)
895  A(i,j) -= A(i,k) * A(j,k);
896 
897  if (i == j)
898  {
899 #ifndef LIBMESH_USE_COMPLEX_NUMBERS
900  libmesh_error_msg_if(A(i,j) <= 0.0,
901  "Error! Can only use Cholesky decomposition with symmetric positive definite matrices.");
902 #endif
903 
904  A(i,i) = std::sqrt(A(i,j));
905  }
906  else
907  A(j,i) = A(i,j) / A(i,i);
908  }
909  }
910 
911  // Set the flag for CHOLESKY decomposition
913 }
DecompositionType _decomposition_type
This flag keeps track of which type of decomposition has been performed on the matrix.
Definition: dense_matrix.h:650
unsigned int m() const
unsigned int n() const

◆ _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.

699 {
700  // This algorithm only works for square matrices, so verify this up front.
701  libmesh_error_msg_if(this->m() != this->n(), "Can only compute eigen-decompositions for square matrices.");
702 
703  // If the user requests left or right eigenvectors, we have to make
704  // sure and pass the transpose of this matrix to Lapack, otherwise
705  // it will compute the inverse transpose of what we are
706  // after... since we know the matrix is square, we can just swap
707  // entries in place. If the user does not request eigenvectors, we
708  // can skip this extra step, since the eigenvalues for A and A^T are
709  // the same.
710  if (VL || VR)
711  {
712  for (auto i : make_range(this->_m))
713  for (auto j : make_range(i))
714  std::swap((*this)(i,j), (*this)(j,i));
715  }
716 
717  // The calling sequence for dgeev is:
718  // DGEEV (character JOBVL,
719  // character JOBVR,
720  // integer N,
721  // double precision, dimension( lda, * ) A,
722  // integer LDA,
723  // double precision, dimension( * ) WR,
724  // double precision, dimension( * ) WI,
725  // double precision, dimension( ldvl, * ) VL,
726  // integer LDVL,
727  // double precision, dimension( ldvr, * ) VR,
728  // integer LDVR,
729  // double precision, dimension( * ) WORK,
730  // integer LWORK,
731  // integer INFO)
732 
733  // JOBVL (input)
734  // = 'N': left eigenvectors of A are not computed;
735  // = 'V': left eigenvectors of A are computed.
736  char JOBVL = VL ? 'V' : 'N';
737 
738  // JOBVR (input)
739  // = 'N': right eigenvectors of A are not computed;
740  // = 'V': right eigenvectors of A are computed.
741  char JOBVR = VR ? 'V' : 'N';
742 
743  // N (input)
744  // The number of rows/cols of the matrix A. N >= 0.
745  PetscBLASInt N = this->m();
746 
747  // A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
748  // On entry, the N-by-N matrix A.
749  // On exit, A has been overwritten.
750 
751  // LDA (input)
752  // The leading dimension of the array A. LDA >= max(1,N).
753  PetscBLASInt LDA = N;
754 
755  // WR (output) double precision array, dimension (N)
756  // WI (output) double precision array, dimension (N)
757  // WR and WI contain the real and imaginary parts,
758  // respectively, of the computed eigenvalues. Complex
759  // conjugate pairs of eigenvalues appear consecutively
760  // with the eigenvalue having the positive imaginary part
761  // first.
762  lambda_real.resize(N);
763  lambda_imag.resize(N);
764 
765  // VL (output) double precision array, dimension (LDVL,N)
766  // If JOBVL = 'V', the left eigenvectors u(j) are stored one
767  // after another in the columns of VL, in the same order
768  // as their eigenvalues.
769  // If JOBVL = 'N', VL is not referenced.
770  // If the j-th eigenvalue is real, then u(j) = VL(:,j),
771  // the j-th column of VL.
772  // If the j-th and (j+1)-st eigenvalues form a complex
773  // conjugate pair, then u(j) = VL(:,j) + i*VL(:,j+1) and
774  // u(j+1) = VL(:,j) - i*VL(:,j+1).
775  // Will be set below if needed.
776 
777  // LDVL (input)
778  // The leading dimension of the array VL. LDVL >= 1; if
779  // JOBVL = 'V', LDVL >= N.
780  PetscBLASInt LDVL = VL ? N : 1;
781 
782  // VR (output) DOUBLE PRECISION array, dimension (LDVR,N)
783  // If JOBVR = 'V', the right eigenvectors v(j) are stored one
784  // after another in the columns of VR, in the same order
785  // as their eigenvalues.
786  // If JOBVR = 'N', VR is not referenced.
787  // If the j-th eigenvalue is real, then v(j) = VR(:,j),
788  // the j-th column of VR.
789  // If the j-th and (j+1)-st eigenvalues form a complex
790  // conjugate pair, then v(j) = VR(:,j) + i*VR(:,j+1) and
791  // v(j+1) = VR(:,j) - i*VR(:,j+1).
792  // Will be set below if needed.
793 
794  // LDVR (input)
795  // The leading dimension of the array VR. LDVR >= 1; if
796  // JOBVR = 'V', LDVR >= N.
797  PetscBLASInt LDVR = VR ? N : 1;
798 
799  // WORK (workspace/output) double precision array, dimension (MAX(1,LWORK))
800  // On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
801  //
802  // LWORK (input)
803  // The dimension of the array WORK. LWORK >= max(1,3*N), and
804  // if JOBVL = 'V' or JOBVR = 'V', LWORK >= 4*N. For good
805  // performance, LWORK must generally be larger.
806  //
807  // If LWORK = -1, then a workspace query is assumed; the routine
808  // only calculates the optimal size of the WORK array, returns
809  // this value as the first entry of the WORK array, and no error
810  // message related to LWORK is issued by XERBLA.
811  PetscBLASInt LWORK = (VR || VL) ? 4*N : 3*N;
812  std::vector<T> WORK(LWORK);
813 
814  // INFO (output)
815  // = 0: successful exit
816  // < 0: if INFO = -i, the i-th argument had an illegal value.
817  // > 0: if INFO = i, the QR algorithm failed to compute all the
818  // eigenvalues, and no eigenvectors or condition numbers
819  // have been computed; elements 1:ILO-1 and i+1:N of WR
820  // and WI contain eigenvalues which have converged.
821  PetscBLASInt INFO = 0;
822 
823  // Get references to raw data
824  std::vector<T> & lambda_real_val = lambda_real.get_values();
825  std::vector<T> & lambda_imag_val = lambda_imag.get_values();
826 
827  // Set up eigenvector storage if necessary.
828  T * VR_ptr = nullptr;
829  if (VR)
830  {
831  VR->resize(N, N);
832  VR_ptr = VR->get_values().data();
833  }
834 
835  T * VL_ptr = nullptr;
836  if (VL)
837  {
838  VL->resize(N, N);
839  VL_ptr = VL->get_values().data();
840  }
841 
842  // Ready to call the Lapack routine through PETSc's interface
843  LAPACKgeev_(&JOBVL,
844  &JOBVR,
845  &N,
846  pPS(_val.data()),
847  &LDA,
848  pPS(lambda_real_val.data()),
849  pPS(lambda_imag_val.data()),
850  pPS(VL_ptr),
851  &LDVL,
852  pPS(VR_ptr),
853  &LDVR,
854  pPS(WORK.data()),
855  &LWORK,
856  &INFO);
857 
858  // Check return value for errors
859  libmesh_error_msg_if(INFO != 0, "INFO=" << INFO << ", Error during Lapack eigenvalue calculation!");
860 
861  // If the user requested either right or left eigenvectors, LAPACK
862  // has now computed the transpose of the desired matrix, i.e. V^T
863  // instead of V. We could leave this up to user code to handle, but
864  // rather than risking getting very unexpected results, we'll just
865  // transpose it in place before handing it back.
866  if (VR)
867  {
868  for (auto i : make_range(N))
869  for (auto j : make_range(i))
870  std::swap((*VR)(i,j), (*VR)(j,i));
871  }
872 
873  if (VL)
874  {
875  for (auto i : make_range(N))
876  for (auto j : make_range(i))
877  std::swap((*VL)(i,j), (*VL)(j,i));
878  }
879 }
unsigned int m() const
PetscScalar * pPS(T *ptr)
Definition: petsc_macro.h:174
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:140
std::vector< T > _val
The actual data values, stored as a 1D array.
Definition: dense_matrix.h:600
unsigned int n() const
unsigned int _m
The row dimension.

◆ _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.

113 {
114  //Check to see if we are doing (A^T)*A
115  if (this == &A)
116  {
117  //libmesh_here();
118  DenseMatrix<T> B(*this);
119 
120  // Simple but inefficient way
121  // return this->left_multiply_transpose(B);
122 
123  // More efficient, but more code way
124  // If A is mxn, the result will be a square matrix of Size n x n.
125  const unsigned int n_rows = A.m();
126  const unsigned int n_cols = A.n();
127 
128  // resize() *this and also zero out all entries.
129  this->resize(n_cols,n_cols);
130 
131  // Compute the lower-triangular part
132  for (unsigned int i=0; i<n_cols; ++i)
133  for (unsigned int j=0; j<=i; ++j)
134  for (unsigned int k=0; k<n_rows; ++k) // inner products are over n_rows
135  (*this)(i,j) += B(k,i)*B(k,j);
136 
137  // Copy lower-triangular part into upper-triangular part
138  for (unsigned int i=0; i<n_cols; ++i)
139  for (unsigned int j=i+1; j<n_cols; ++j)
140  (*this)(i,j) = (*this)(j,i);
141  }
142 
143  else
144  {
145  DenseMatrix<T> B(*this);
146 
147  this->resize (A.n(), B.n());
148 
149  libmesh_assert_equal_to (A.m(), B.m());
150  libmesh_assert_equal_to (this->m(), A.n());
151  libmesh_assert_equal_to (this->n(), B.n());
152 
153  const unsigned int m_s = A.n();
154  const unsigned int p_s = A.m();
155  const unsigned int n_s = this->n();
156 
157  // Do it this way because there is a
158  // decent chance (at least for constraint matrices)
159  // that A.transpose(i,k) = 0.
160  for (unsigned int i=0; i<m_s; i++)
161  for (unsigned int k=0; k<p_s; k++)
162  if (A.transpose(i,k) != 0.)
163  for (unsigned int j=0; j<n_s; j++)
164  (*this)(i,j) += A.transpose(i,k)*B(k,j);
165  }
166 }
unsigned int m() const
Definition: assembly.h:38
void resize(const unsigned int new_m, const unsigned int new_n)
Resizes the matrix to the specified size and calls zero().
Definition: dense_matrix.h:895
unsigned int n() const

◆ _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.

582 {
583  const unsigned int
584  n_cols = this->n();
585 
586  libmesh_assert_equal_to (this->m(), n_cols);
587  libmesh_assert_equal_to (this->m(), b.size());
588 
589  x.resize (n_cols);
590 
591  // A convenient reference to *this
592  const DenseMatrix<T> & A = *this;
593 
594  // Temporary vector storage. We use this instead of
595  // modifying the RHS.
596  DenseVector<T> z = b;
597 
598  // Lower-triangular "top to bottom" solve step, taking into account pivots
599  for (unsigned int i=0; i<n_cols; ++i)
600  {
601  // Swap
602  if (_pivots[i] != static_cast<pivot_index_t>(i))
603  std::swap( z(i), z(_pivots[i]) );
604 
605  x(i) = z(i);
606 
607  for (unsigned int j=0; j<i; ++j)
608  x(i) -= A(i,j)*x(j);
609 
610  x(i) /= A(i,i);
611  }
612 
613  // Upper-triangular "bottom to top" solve step
614  const unsigned int last_row = n_cols-1;
615 
616  for (int i=last_row; i>=0; --i)
617  {
618  for (int j=i+1; j<static_cast<int>(n_cols); ++j)
619  x(i) -= A(i,j)*x(j);
620  }
621 }
unsigned int m() const
unsigned int n() const
std::vector< pivot_index_t > _pivots
Definition: dense_matrix.h:740

◆ _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.

903 {
904  // The calling sequence for getrs is:
905  // dgetrs(TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
906 
907  // trans (input)
908  // 'n' for no transpose, 't' for transpose
909  char TRANS[] = "t";
910 
911  // N (input)
912  // The order of the matrix A. N >= 0.
913  PetscBLASInt N = this->m();
914 
915 
916  // NRHS (input)
917  // The number of right hand sides, i.e., the number of columns
918  // of the matrix B. NRHS >= 0.
919  PetscBLASInt NRHS = 1;
920 
921  // A (input) double precision array, dimension (LDA,N)
922  // The factors L and U from the factorization A = P*L*U
923  // as computed by dgetrf.
924 
925  // LDA (input)
926  // The leading dimension of the array A. LDA >= max(1,N).
927  PetscBLASInt LDA = N;
928 
929  // ipiv (input) int array, dimension (N)
930  // The pivot indices from DGETRF; for 1<=i<=N, row i of the
931  // matrix was interchanged with row IPIV(i).
932  // Here, we pass _pivots.data() which was computed in _lu_decompose_lapack
933 
934  // B (input/output) double precision array, dimension (LDB,NRHS)
935  // On entry, the right hand side matrix B.
936  // On exit, the solution matrix X.
937  // Here, we pass a copy of the rhs vector's data array in x, so that the
938  // passed right-hand side b is unmodified. I don't see a way around this
939  // copy if we want to maintain an unmodified rhs in LibMesh.
940  x = b;
941  std::vector<T> & x_vec = x.get_values();
942 
943  // We can avoid the copy if we don't care about overwriting the RHS: just
944  // pass b to the Lapack routine and then swap with x before exiting
945  // std::vector<T> & x_vec = b.get_values();
946 
947  // LDB (input)
948  // The leading dimension of the array B. LDB >= max(1,N).
949  PetscBLASInt LDB = N;
950 
951  // INFO (output)
952  // = 0: successful exit
953  // < 0: if INFO = -i, the i-th argument had an illegal value
954  PetscBLASInt INFO = 0;
955 
956  // Finally, ready to call the Lapack getrs function
957  LAPACKgetrs_(TRANS, &N, &NRHS, pPS(_val.data()), &LDA,
958  _pivots.data(), pPS(x_vec.data()), &LDB, &INFO);
959 
960  // Check return value for errors
961  libmesh_error_msg_if(INFO != 0, "INFO=" << INFO << ", Error during Lapack LU solve!");
962 
963  // Don't do this if you already made a copy of b above
964  // Swap b and x. The solution will then be in x, and whatever was originally
965  // in x, maybe garbage, maybe nothing, will be in b.
966  // FIXME: Rewrite the LU and Cholesky solves to just take one input, and overwrite
967  // the input. This *should* make user code simpler, as they don't have to create
968  // an extra vector just to pass it in to the solve function!
969  // b.swap(x);
970 }
unsigned int m() const
PetscScalar * pPS(T *ptr)
Definition: petsc_macro.h:174
std::vector< T > _val
The actual data values, stored as a 1D array.
Definition: dense_matrix.h:600
std::vector< pivot_index_t > _pivots
Definition: dense_matrix.h:740

◆ _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.

632 {
633  // If this function was called, there better not be any
634  // previous decomposition of the matrix.
635  libmesh_assert_equal_to (this->_decomposition_type, NONE);
636 
637  // Get the matrix size and make sure it is square
638  const unsigned int
639  n_rows = this->m();
640 
641  // A convenient reference to *this
642  DenseMatrix<T> & A = *this;
643 
644  _pivots.resize(n_rows);
645 
646  for (unsigned int i=0; i<n_rows; ++i)
647  {
648  // Find the pivot row by searching down the i'th column
649  _pivots[i] = i;
650 
651  // std::abs(complex) must return a Real!
652  auto the_max = std::abs( A(i,i) );
653  for (unsigned int j=i+1; j<n_rows; ++j)
654  {
655  auto candidate_max = std::abs( A(j,i) );
656  if (the_max < candidate_max)
657  {
658  the_max = candidate_max;
659  _pivots[i] = j;
660  }
661  }
662 
663  // libMesh::out << "the_max=" << the_max << " found at row " << _pivots[i] << std::endl;
664 
665  // If the max was found in a different row, interchange rows.
666  // Here we interchange the *entire* row, in Gaussian elimination
667  // you would only interchange the subrows A(i,j) and A(p(i),j), for j>i
668  if (_pivots[i] != static_cast<pivot_index_t>(i))
669  {
670  for (unsigned int j=0; j<n_rows; ++j)
671  std::swap( A(i,j), A(_pivots[i], j) );
672  }
673 
674  // If the max abs entry found is zero, the matrix is singular
675  libmesh_error_msg_if(A(i,i) == libMesh::zero, "Matrix A is singular!");
676 
677  // Scale upper triangle entries of row i by the diagonal entry
678  // Note: don't scale the diagonal entry itself!
679  const T diag_inv = 1. / A(i,i);
680  for (unsigned int j=i+1; j<n_rows; ++j)
681  A(i,j) *= diag_inv;
682 
683  // Update the remaining sub-matrix A[i+1:m][i+1:m]
684  // by subtracting off (the diagonal-scaled)
685  // upper-triangular part of row i, scaled by the
686  // i'th column entry of each row. In terms of
687  // row operations, this is:
688  // for each r > i
689  // SubRow(r) = SubRow(r) - A(r,i)*SubRow(i)
690  //
691  // If we were scaling the i'th column as well, like
692  // in Gaussian elimination, this would 'zero' the
693  // entry in the i'th column.
694  for (unsigned int row=i+1; row<n_rows; ++row)
695  for (unsigned int col=i+1; col<n_rows; ++col)
696  A(row,col) -= A(row,i) * A(i,col);
697 
698  } // end i loop
699 
700  // Set the flag for LU decomposition
701  this->_decomposition_type = LU;
702 }
DecompositionType _decomposition_type
This flag keeps track of which type of decomposition has been performed on the matrix.
Definition: dense_matrix.h:650
unsigned int m() const
const Number zero
.
Definition: libmesh.h:304
std::vector< pivot_index_t > _pivots
Definition: dense_matrix.h:740

◆ _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.

211 {
212  // If this function was called, there better not be any
213  // previous decomposition of the matrix.
214  libmesh_assert_equal_to (this->_decomposition_type, NONE);
215 
216  // The calling sequence for dgetrf is:
217  // dgetrf(M, N, A, lda, ipiv, info)
218 
219  // M (input)
220  // The number of rows of the matrix A. M >= 0.
221  // In C/C++, pass the number of *cols* of A
222  PetscBLASInt M = this->n();
223 
224  // N (input)
225  // The number of columns of the matrix A. N >= 0.
226  // In C/C++, pass the number of *rows* of A
227  PetscBLASInt N = this->m();
228 
229  // A (input/output) double precision array, dimension (LDA,N)
230  // On entry, the M-by-N matrix to be factored.
231  // On exit, the factors L and U from the factorization
232  // A = P*L*U; the unit diagonal elements of L are not stored.
233 
234  // LDA (input)
235  // The leading dimension of the array A. LDA >= max(1,M).
236  PetscBLASInt LDA = M;
237 
238  // ipiv (output) integer array, dimension (min(m,n))
239  // The pivot indices; for 1 <= i <= min(m,n), row i of the
240  // matrix was interchanged with row IPIV(i).
241  // Here, we pass _pivots.data(), a private class member used to store pivots
242  this->_pivots.resize( std::min(M,N) );
243 
244  // info (output)
245  // = 0: successful exit
246  // < 0: if INFO = -i, the i-th argument had an illegal value
247  // > 0: if INFO = i, U(i,i) is exactly zero. The factorization
248  // has been completed, but the factor U is exactly
249  // singular, and division by zero will occur if it is used
250  // to solve a system of equations.
251  PetscBLASInt INFO = 0;
252 
253  // Ready to call the actual factorization routine through PETSc's interface
254  LAPACKgetrf_(&M, &N, pPS(this->_val.data()), &LDA, _pivots.data(),
255  &INFO);
256 
257  // Check return value for errors
258  libmesh_error_msg_if(INFO != 0, "INFO=" << INFO << ", Error during Lapack LU factorization!");
259 
260  // Set the flag for LU decomposition
262 }
DecompositionType _decomposition_type
This flag keeps track of which type of decomposition has been performed on the matrix.
Definition: dense_matrix.h:650
unsigned int m() const
PetscScalar * pPS(T *ptr)
Definition: petsc_macro.h:174
std::vector< T > _val
The actual data values, stored as a 1D array.
Definition: dense_matrix.h:600
unsigned int n() const
std::vector< pivot_index_t > _pivots
Definition: dense_matrix.h:740

◆ _matvec_blas()

template<typename T>
template LIBMESH_EXPORT void libMesh::DenseMatrix< T >::_matvec_blas ( alpha,
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.

995 {
996  // Ensure that dest and arg sizes are compatible
997  if (!trans)
998  {
999  // dest ~ A * arg
1000  // (mx1) (mxn) * (nx1)
1001  libmesh_error_msg_if((dest.size() != this->m()) || (arg.size() != this->n()),
1002  "Improper input argument sizes!");
1003  }
1004 
1005  else // trans == true
1006  {
1007  // Ensure that dest and arg are proper size
1008  // dest ~ A^T * arg
1009  // (nx1) (nxm) * (mx1)
1010  libmesh_error_msg_if((dest.size() != this->n()) || (arg.size() != this->m()),
1011  "Improper input argument sizes!");
1012  }
1013 
1014  // Calling sequence for dgemv:
1015  //
1016  // dgemv(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
1017 
1018  // TRANS (input)
1019  // 't' for transpose, 'n' for non-transpose multiply
1020  // We store everything in row-major order, so pass the transpose flag for
1021  // non-transposed matvecs and the 'n' flag for transposed matvecs
1022  char TRANS[] = "t";
1023  if (trans)
1024  TRANS[0] = 'n';
1025 
1026  // M (input)
1027  // On entry, M specifies the number of rows of the matrix A.
1028  // In C/C++, pass the number of *cols* of A
1029  PetscBLASInt M = this->n();
1030 
1031  // N (input)
1032  // On entry, N specifies the number of columns of the matrix A.
1033  // In C/C++, pass the number of *rows* of A
1034  PetscBLASInt N = this->m();
1035 
1036  // ALPHA (input)
1037  // The scalar constant passed to this function
1038 
1039  // A (input) double precision array of DIMENSION ( LDA, n ).
1040  // Before entry, the leading m by n part of the array A must
1041  // contain the matrix of coefficients.
1042  // The matrix, *this. Note that _matvec_blas is called from
1043  // a const function, vector_mult(), and so we have made this function const
1044  // as well. Since BLAS knows nothing about const, we have to cast it away
1045  // now.
1046  DenseMatrix<T> & a_ref = const_cast<DenseMatrix<T> &> ( *this );
1047  std::vector<T> & a = a_ref.get_values();
1048 
1049  // LDA (input)
1050  // On entry, LDA specifies the first dimension of A as declared
1051  // in the calling (sub) program. LDA must be at least
1052  // max( 1, m ).
1053  PetscBLASInt LDA = M;
1054 
1055  // X (input) double precision array of DIMENSION at least
1056  // ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
1057  // and at least
1058  // ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
1059  // Before entry, the incremented array X must contain the
1060  // vector x.
1061  // Here, we must cast away the const-ness of "arg" since BLAS knows
1062  // nothing about const
1063  DenseVector<T> & x_ref = const_cast<DenseVector<T> &> ( arg );
1064  std::vector<T> & x = x_ref.get_values();
1065 
1066  // INCX (input)
1067  // On entry, INCX specifies the increment for the elements of
1068  // X. INCX must not be zero.
1069  PetscBLASInt INCX = 1;
1070 
1071  // BETA (input)
1072  // On entry, BETA specifies the scalar beta. When BETA is
1073  // supplied as zero then Y need not be set on input.
1074  // The second scalar constant passed to this function
1075 
1076  // Y (input) double precision array of DIMENSION at least
1077  // ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
1078  // and at least
1079  // ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
1080  // Before entry with BETA non-zero, the incremented array Y
1081  // must contain the vector y. On exit, Y is overwritten by the
1082  // updated vector y.
1083  // The input vector "dest"
1084  std::vector<T> & y = dest.get_values();
1085 
1086  // INCY (input)
1087  // On entry, INCY specifies the increment for the elements of
1088  // Y. INCY must not be zero.
1089  PetscBLASInt INCY = 1;
1090 
1091  // Finally, ready to call the BLAS function
1092  BLASgemv_(TRANS, &M, &N, pPS(&alpha), pPS(a.data()), &LDA,
1093  pPS(x.data()), &INCX, pPS(&beta), pPS(y.data()), &INCY);
1094 }
unsigned int m() const
PetscScalar * pPS(T *ptr)
Definition: petsc_macro.h:174
unsigned int n() const

◆ _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.

39 {
40  int result_size = 0;
41 
42  // For each case, determine the size of the final result make sure
43  // that the inner dimensions match
44  switch (flag)
45  {
46  case LEFT_MULTIPLY:
47  {
48  result_size = other.m() * this->n();
49  if (other.n() == this->m())
50  break;
51  }
52  libmesh_fallthrough();
53  case RIGHT_MULTIPLY:
54  {
55  result_size = other.n() * this->m();
56  if (other.m() == this->n())
57  break;
58  }
59  libmesh_fallthrough();
61  {
62  result_size = other.n() * this->n();
63  if (other.m() == this->m())
64  break;
65  }
66  libmesh_fallthrough();
68  {
69  result_size = other.m() * this->m();
70  if (other.n() == this->n())
71  break;
72  }
73  libmesh_fallthrough();
74  default:
75  libmesh_error_msg("Unknown flag selected or matrices are incompatible for multiplication.");
76  }
77 
78  // For this to work, the passed arg. must actually be a DenseMatrix<T>
79  const DenseMatrix<T> * const_that = cast_ptr<const DenseMatrix<T> *>(&other);
80 
81  // Also, although 'that' is logically const in this BLAS routine,
82  // the PETSc BLAS interface does not specify that any of the inputs are
83  // const. To use it, I must cast away const-ness.
84  DenseMatrix<T> * that = const_cast<DenseMatrix<T> *> (const_that);
85 
86  // Initialize A, B pointers for LEFT_MULTIPLY* cases
87  DenseMatrix<T> * A = this;
88  DenseMatrix<T> * B = that;
89 
90  // For RIGHT_MULTIPLY* cases, swap the meaning of A and B.
91  // Here is a full table of combinations we can pass to BLASgemm, and what the answer is when finished:
92  // pass A B -> (Fortran) -> A^T B^T -> (C++) -> (A^T B^T)^T -> (identity) -> B A "lt multiply"
93  // pass B A -> (Fortran) -> B^T A^T -> (C++) -> (B^T A^T)^T -> (identity) -> A B "rt multiply"
94  // pass A B^T -> (Fortran) -> A^T B -> (C++) -> (A^T B)^T -> (identity) -> B^T A "lt multiply t"
95  // pass B^T A -> (Fortran) -> B A^T -> (C++) -> (B A^T)^T -> (identity) -> A B^T "rt multiply t"
96  if (flag==RIGHT_MULTIPLY || flag==RIGHT_MULTIPLY_TRANSPOSE)
97  std::swap(A,B);
98 
99  // transa, transb values to pass to blas
100  char
101  transa[] = "n",
102  transb[] = "n";
103 
104  // Integer values to pass to BLAS:
105  //
106  // M
107  // In Fortran, the number of rows of op(A),
108  // In the BLAS documentation, typically known as 'M'.
109  //
110  // In C/C++, we set:
111  // M = n_cols(A) if (transa='n')
112  // n_rows(A) if (transa='t')
113  PetscBLASInt M = static_cast<PetscBLASInt>( A->n() );
114 
115  // N
116  // In Fortran, the number of cols of op(B), and also the number of cols of C.
117  // In the BLAS documentation, typically known as 'N'.
118  //
119  // In C/C++, we set:
120  // N = n_rows(B) if (transb='n')
121  // n_cols(B) if (transb='t')
122  PetscBLASInt N = static_cast<PetscBLASInt>( B->m() );
123 
124  // K
125  // In Fortran, the number of cols of op(A), and also
126  // the number of rows of op(B). In the BLAS documentation,
127  // typically known as 'K'.
128  //
129  // In C/C++, we set:
130  // K = n_rows(A) if (transa='n')
131  // n_cols(A) if (transa='t')
132  PetscBLASInt K = static_cast<PetscBLASInt>( A->m() );
133 
134  // LDA (leading dimension of A). In our cases,
135  // LDA is always the number of columns of A.
136  PetscBLASInt LDA = static_cast<PetscBLASInt>( A->n() );
137 
138  // LDB (leading dimension of B). In our cases,
139  // LDB is always the number of columns of B.
140  PetscBLASInt LDB = static_cast<PetscBLASInt>( B->n() );
141 
142  if (flag == LEFT_MULTIPLY_TRANSPOSE)
143  {
144  transb[0] = 't';
145  N = static_cast<PetscBLASInt>( B->n() );
146  }
147 
148  else if (flag == RIGHT_MULTIPLY_TRANSPOSE)
149  {
150  transa[0] = 't';
151  std::swap(M,K);
152  }
153 
154  // LDC (leading dimension of C). LDC is the
155  // number of columns in the solution matrix.
156  PetscBLASInt LDC = M;
157 
158  // Scalar values to pass to BLAS
159  //
160  // scalar multiplying the whole product AB
161  T alpha = 1.;
162 
163  // scalar multiplying C, which is the original matrix.
164  T beta = 0.;
165 
166  // Storage for the result
167  std::vector<T> result (result_size);
168 
169  // Finally ready to call the BLAS
170  BLASgemm_(transa, transb, &M, &N, &K, pPS(&alpha),
171  pPS(A->get_values().data()), &LDA,
172  pPS(B->get_values().data()), &LDB, pPS(&beta),
173  pPS(result.data()), &LDC);
174 
175  // Update the relevant dimension for this matrix.
176  switch (flag)
177  {
178  case LEFT_MULTIPLY: { this->_m = other.m(); break; }
179  case RIGHT_MULTIPLY: { this->_n = other.n(); break; }
180  case LEFT_MULTIPLY_TRANSPOSE: { this->_m = other.n(); break; }
181  case RIGHT_MULTIPLY_TRANSPOSE: { this->_n = other.m(); break; }
182  default:
183  libmesh_error_msg("Unknown flag selected.");
184  }
185 
186  // Swap my data vector with the result
187  this->_val.swap(result);
188 }
unsigned int m() const
PetscScalar * pPS(T *ptr)
Definition: petsc_macro.h:174
Definition: assembly.h:38
unsigned int _n
The column dimension.
std::vector< T > _val
The actual data values, stored as a 1D array.
Definition: dense_matrix.h:600
unsigned int n() const
unsigned int _m
The row dimension.

◆ _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.

240 {
241  //Check to see if we are doing B*(B^T)
242  if (this == &B)
243  {
244  //libmesh_here();
245  DenseMatrix<T> A(*this);
246 
247  // Simple but inefficient way
248  // return this->right_multiply_transpose(A);
249 
250  // More efficient, more code
251  // If B is mxn, the result will be a square matrix of Size m x m.
252  const unsigned int n_rows = B.m();
253  const unsigned int n_cols = B.n();
254 
255  // resize() *this and also zero out all entries.
256  this->resize(n_rows,n_rows);
257 
258  // Compute the lower-triangular part
259  for (unsigned int i=0; i<n_rows; ++i)
260  for (unsigned int j=0; j<=i; ++j)
261  for (unsigned int k=0; k<n_cols; ++k) // inner products are over n_cols
262  (*this)(i,j) += A(i,k)*A(j,k);
263 
264  // Copy lower-triangular part into upper-triangular part
265  for (unsigned int i=0; i<n_rows; ++i)
266  for (unsigned int j=i+1; j<n_rows; ++j)
267  (*this)(i,j) = (*this)(j,i);
268  }
269 
270  else
271  {
272  DenseMatrix<T> A(*this);
273 
274  this->resize (A.m(), B.m());
275 
276  libmesh_assert_equal_to (A.n(), B.n());
277  libmesh_assert_equal_to (this->m(), A.m());
278  libmesh_assert_equal_to (this->n(), B.m());
279 
280  const unsigned int m_s = A.m();
281  const unsigned int p_s = A.n();
282  const unsigned int n_s = this->n();
283 
284  // Do it this way because there is a
285  // decent chance (at least for constraint matrices)
286  // that B.transpose(k,j) = 0.
287  for (unsigned int j=0; j<n_s; j++)
288  for (unsigned int k=0; k<p_s; k++)
289  if (B.transpose(k,j) != 0.)
290  for (unsigned int i=0; i<m_s; i++)
291  (*this)(i,j) += A(i,k)*B.transpose(k,j);
292  }
293 }
unsigned int m() const
Definition: assembly.h:38
void resize(const unsigned int new_m, const unsigned int new_n)
Resizes the matrix to the specified size and calls zero().
Definition: dense_matrix.h:895
unsigned int n() const

◆ _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.

389 {
390 
391  // M (input)
392  // The number of rows of the matrix A. M >= 0.
393  // In C/C++, pass the number of *cols* of A
394  PetscBLASInt M = this->n();
395 
396  // N (input)
397  // The number of columns of the matrix A. N >= 0.
398  // In C/C++, pass the number of *rows* of A
399  PetscBLASInt N = this->m();
400 
401  PetscBLASInt min_MN = (M < N) ? M : N;
402  PetscBLASInt max_MN = (M > N) ? M : N;
403 
404  // A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
405  // On entry, the M-by-N matrix A.
406  // On exit,
407  // if JOBU = 'O', A is overwritten with the first min(m,n)
408  // columns of U (the left singular vectors,
409  // stored columnwise);
410  // if JOBVT = 'O', A is overwritten with the first min(m,n)
411  // rows of V**T (the right singular vectors,
412  // stored rowwise);
413  // if JOBU != 'O' and JOBVT != 'O', the contents of A are destroyed.
414 
415  // LDA (input)
416  // The leading dimension of the array A. LDA >= max(1,M).
417  PetscBLASInt LDA = M;
418 
419  // S (output) DOUBLE PRECISION array, dimension (min(M,N))
420  // The singular values of A, sorted so that S(i) >= S(i+1).
421  sigma_val.resize( min_MN );
422 
423  // LDU (input)
424  // The leading dimension of the array U. LDU >= 1; if
425  // JOBU = 'S' or 'A', LDU >= M.
426  PetscBLASInt LDU = M;
427 
428  // U (output) DOUBLE PRECISION array, dimension (LDU,UCOL)
429  // (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'.
430  // If JOBU = 'A', U contains the M-by-M orthogonal matrix U;
431  // if JOBU = 'S', U contains the first min(m,n) columns of U
432  // (the left singular vectors, stored columnwise);
433  // if JOBU = 'N' or 'O', U is not referenced.
434  if (JOBU == 'S')
435  U_val.resize( LDU*min_MN );
436  else
437  U_val.resize( LDU*M );
438 
439  // LDVT (input)
440  // The leading dimension of the array VT. LDVT >= 1; if
441  // JOBVT = 'A', LDVT >= N; if JOBVT = 'S', LDVT >= min(M,N).
442  PetscBLASInt LDVT = N;
443  if (JOBVT == 'S')
444  LDVT = min_MN;
445 
446  // VT (output) DOUBLE PRECISION array, dimension (LDVT,N)
447  // If JOBVT = 'A', VT contains the N-by-N orthogonal matrix
448  // V**T;
449  // if JOBVT = 'S', VT contains the first min(m,n) rows of
450  // V**T (the right singular vectors, stored rowwise);
451  // if JOBVT = 'N' or 'O', VT is not referenced.
452  VT_val.resize( LDVT*N );
453 
454  // LWORK (input)
455  // The dimension of the array WORK.
456  // LWORK >= MAX(1,3*MIN(M,N)+MAX(M,N),5*MIN(M,N)).
457  // For good performance, LWORK should generally be larger.
458  //
459  // If LWORK = -1, then a workspace query is assumed; the routine
460  // only calculates the optimal size of the WORK array, returns
461  // this value as the first entry of the WORK array, and no error
462  // message related to LWORK is issued by XERBLA.
463  PetscBLASInt larger = (3*min_MN+max_MN > 5*min_MN) ? 3*min_MN+max_MN : 5*min_MN;
464  PetscBLASInt LWORK = (larger > 1) ? larger : 1;
465 
466 
467  // WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
468  // On exit, if INFO = 0, WORK(1) returns the optimal LWORK;
469  // if INFO > 0, WORK(2:MIN(M,N)) contains the unconverged
470  // superdiagonal elements of an upper bidiagonal matrix B
471  // whose diagonal is in S (not necessarily sorted). B
472  // satisfies A = U * B * VT, so it has the same singular values
473  // as A, and singular vectors related by U and VT.
474  std::vector<Number> WORK( LWORK );
475 
476  // INFO (output)
477  // = 0: successful exit.
478  // < 0: if INFO = -i, the i-th argument had an illegal value.
479  // > 0: if DBDSQR did not converge, INFO specifies how many
480  // superdiagonals of an intermediate bidiagonal form B
481  // did not converge to zero. See the description of WORK
482  // above for details.
483  PetscBLASInt INFO = 0;
484 
485  // Ready to call the actual factorization routine through PETSc's interface.
486 #ifdef LIBMESH_USE_REAL_NUMBERS
487  // Note that the call to LAPACKgesvd_ may modify _val
488  LAPACKgesvd_(&JOBU, &JOBVT, &M, &N, pPS(_val.data()), &LDA,
489  pPS(sigma_val.data()), pPS(U_val.data()), &LDU,
490  pPS(VT_val.data()), &LDVT, pPS(WORK.data()), &LWORK,
491  &INFO);
492 #else
493  // When we have LIBMESH_USE_COMPLEX_NUMBERS then we must pass an array of Complex
494  // numbers to LAPACKgesvd_, but _val may contain Reals so we copy to Number below to
495  // handle both the real-valued and complex-valued cases.
496  std::vector<Number> val_copy(_val.size());
497  for (auto i : make_range(_val.size()))
498  val_copy[i] = _val[i];
499 
500  std::vector<Real> RWORK(5 * min_MN);
501  LAPACKgesvd_(&JOBU, &JOBVT, &M, &N, val_copy.data(), &LDA, sigma_val.data(), U_val.data(),
502  &LDU, VT_val.data(), &LDVT, WORK.data(), &LWORK, RWORK.data(), &INFO);
503 #endif
504 
505  // Check return value for errors
506  libmesh_error_msg_if(INFO != 0, "INFO=" << INFO << ", Error during Lapack SVD calculation!");
507 }
unsigned int m() const
PetscScalar * pPS(T *ptr)
Definition: petsc_macro.h:174
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:140
std::vector< T > _val
The actual data values, stored as a 1D array.
Definition: dense_matrix.h:600
unsigned int n() const

◆ _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.

278 {
279  // The calling sequence for dgetrf is:
280  // DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, INFO )
281 
282  // JOBU (input)
283  // Specifies options for computing all or part of the matrix U:
284  // = 'A': all M columns of U are returned in array U:
285  // = 'S': the first min(m,n) columns of U (the left singular
286  // vectors) are returned in the array U;
287  // = 'O': the first min(m,n) columns of U (the left singular
288  // vectors) are overwritten on the array A;
289  // = 'N': no columns of U (no left singular vectors) are
290  // computed.
291  char JOBU = 'N';
292 
293  // JOBVT (input)
294  // Specifies options for computing all or part of the matrix
295  // V**T:
296  // = 'A': all N rows of V**T are returned in the array VT;
297  // = 'S': the first min(m,n) rows of V**T (the right singular
298  // vectors) are returned in the array VT;
299  // = 'O': the first min(m,n) rows of V**T (the right singular
300  // vectors) are overwritten on the array A;
301  // = 'N': no rows of V**T (no right singular vectors) are
302  // computed.
303  char JOBVT = 'N';
304 
305  std::vector<Real> sigma_val;
306  std::vector<Number> U_val;
307  std::vector<Number> VT_val;
308 
309  _svd_helper(JOBU, JOBVT, sigma_val, U_val, VT_val);
310 
311  // Copy the singular values into sigma, ignore U_val and VT_val
312  sigma.resize(cast_int<unsigned int>(sigma_val.size()));
313  for (auto i : make_range(sigma.size()))
314  sigma(i) = sigma_val[i];
315 }
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:140
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.

◆ _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.

321 {
322  // The calling sequence for dgetrf is:
323  // DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, INFO )
324 
325  // JOBU (input)
326  // Specifies options for computing all or part of the matrix U:
327  // = 'A': all M columns of U are returned in array U:
328  // = 'S': the first min(m,n) columns of U (the left singular
329  // vectors) are returned in the array U;
330  // = 'O': the first min(m,n) columns of U (the left singular
331  // vectors) are overwritten on the array A;
332  // = 'N': no columns of U (no left singular vectors) are
333  // computed.
334  char JOBU = 'S';
335 
336  // JOBVT (input)
337  // Specifies options for computing all or part of the matrix
338  // V**T:
339  // = 'A': all N rows of V**T are returned in the array VT;
340  // = 'S': the first min(m,n) rows of V**T (the right singular
341  // vectors) are returned in the array VT;
342  // = 'O': the first min(m,n) rows of V**T (the right singular
343  // vectors) are overwritten on the array A;
344  // = 'N': no rows of V**T (no right singular vectors) are
345  // computed.
346  char JOBVT = 'S';
347 
348  // Note: Lapack is going to compute the singular values of A^T. If
349  // A=U * S * V^T, then A^T = V * S * U^T, which means that the
350  // values returned in the "U_val" array actually correspond to the
351  // entries of the V matrix, and the values returned in the VT_val
352  // array actually correspond to the entries of U^T. Therefore, we
353  // pass VT in the place of U and U in the place of VT below!
354  std::vector<Real> sigma_val;
355  int M = this->n();
356  int N = this->m();
357  int min_MN = (M < N) ? M : N;
358 
359  // Size user-provided storage appropriately. Inside svd_helper:
360  // U_val is sized to (M x min_MN)
361  // VT_val is sized to (min_MN x N)
362  // So, we set up U to have the shape of "VT_val^T", and VT to
363  // have the shape of "U_val^T".
364  //
365  // Finally, since the results are stored in column-major order by
366  // Lapack, but we actually want the transpose of what Lapack
367  // returns, this means (conveniently) that we don't even have to
368  // copy anything after the call to _svd_helper, it should already be
369  // in the correct order!
370  U.resize(N, min_MN);
371  VT.resize(min_MN, M);
372 
373  _svd_helper(JOBU, JOBVT, sigma_val, VT.get_values(), U.get_values());
374 
375  // Copy the singular values into sigma.
376  sigma.resize(cast_int<unsigned int>(sigma_val.size()));
377  for (auto i : make_range(sigma.size()))
378  sigma(i) = sigma_val[i];
379 }
unsigned int m() const
std::vector< T > & get_values()
Definition: dense_matrix.h:382
void resize(const unsigned int new_m, const unsigned int new_n)
Resizes the matrix to the specified size and calls zero().
Definition: dense_matrix.h:895
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:140
unsigned int n() const
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.

◆ _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.

532 {
533  // Since BLAS is expecting column-major storage, we first need to
534  // make a transposed copy of *this, then pass it to the gelss
535  // routine instead of the original. This extra copy is kind of a
536  // bummer, it might be better if we could use the full SVD to
537  // compute the least-squares solution instead... Note that it isn't
538  // completely terrible either, since A_trans gets overwritten by
539  // Lapack, and we usually would end up making a copy of A outside
540  // the function call anyway.
541  DenseMatrix<T> A_trans;
542  this->get_transpose(A_trans);
543 
544  // M
545  // The number of rows of the input matrix. M >= 0.
546  // This is actually the number of *columns* of A_trans.
547  PetscBLASInt M = A_trans.n();
548 
549  // N
550  // The number of columns of the matrix A. N >= 0.
551  // This is actually the number of *rows* of A_trans.
552  PetscBLASInt N = A_trans.m();
553 
554  // We'll use the min and max of (M,N) several times below.
555  PetscBLASInt max_MN = std::max(M,N);
556  PetscBLASInt min_MN = std::min(M,N);
557 
558  // NRHS
559  // The number of right hand sides, i.e., the number of columns
560  // of the matrices B and X. NRHS >= 0.
561  // This could later be generalized to solve for multiple right-hand
562  // sides...
563  PetscBLASInt NRHS = 1;
564 
565  // A is double precision array, dimension (LDA,N)
566  // On entry, the M-by-N matrix A.
567  // On exit, the first min(m,n) rows of A are overwritten with
568  // its right singular vectors, stored rowwise.
569  //
570  // The data vector that will be passed to Lapack.
571  std::vector<T> & A_trans_vals = A_trans.get_values();
572 
573  // LDA
574  // The leading dimension of the array A. LDA >= max(1,M).
575  PetscBLASInt LDA = M;
576 
577  // B is double precision array, dimension (LDB,NRHS)
578  // On entry, the M-by-NRHS right hand side matrix B.
579  // On exit, B is overwritten by the N-by-NRHS solution
580  // matrix X. If m >= n and RANK = n, the residual
581  // sum-of-squares for the solution in the i-th column is given
582  // by the sum of squares of elements n+1:m in that column.
583  //
584  // Since we don't want the user's rhs vector to be overwritten by
585  // the solution, we copy the rhs values into the solution vector "x"
586  // now. x needs to be long enough to hold both the (Nx1) solution
587  // vector or the (Mx1) rhs, so size it to the max of those.
588  x.resize(max_MN);
589  for (auto i : make_range(rhs.size()))
590  x(i) = rhs(i);
591 
592  // Make the syntax below simpler by grabbing a reference to this array.
593  std::vector<T> & B = x.get_values();
594 
595  // LDB
596  // The leading dimension of the array B. LDB >= max(1,max(M,N)).
597  PetscBLASInt LDB = x.size();
598 
599  // S is double precision array, dimension (min(M,N))
600  // The singular values of A in decreasing order.
601  // The condition number of A in the 2-norm = S(1)/S(min(m,n)).
602  std::vector<T> S(min_MN);
603 
604  // RCOND
605  // Used to determine the effective rank of A. Singular values
606  // S(i) <= RCOND*S(1) are treated as zero. If RCOND < 0, machine
607  // precision is used instead.
608  PetscScalar RCOND = PS(rcond);
609 
610  // RANK
611  // The effective rank of A, i.e., the number of singular values
612  // which are greater than RCOND*S(1).
613  PetscBLASInt RANK = 0;
614 
615  // LWORK
616  // The dimension of the array WORK. LWORK >= 1, and also:
617  // LWORK >= 3*min(M,N) + max( 2*min(M,N), max(M,N), NRHS )
618  // For good performance, LWORK should generally be larger.
619  //
620  // If LWORK = -1, then a workspace query is assumed; the routine
621  // only calculates the optimal size of the WORK array, returns
622  // this value as the first entry of the WORK array, and no error
623  // message related to LWORK is issued by XERBLA.
624  //
625  // The factor of 3/2 is arbitrary and is used to satisfy the "should
626  // generally be larger" clause.
627  PetscBLASInt LWORK = (3*min_MN + std::max(2*min_MN, std::max(max_MN, NRHS))) * 3/2;
628 
629  // WORK is double precision array, dimension (MAX(1,LWORK))
630  // On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
631  std::vector<T> WORK(LWORK);
632 
633  // INFO
634  // = 0: successful exit
635  // < 0: if INFO = -i, the i-th argument had an illegal value.
636  // > 0: the algorithm for computing the SVD failed to converge;
637  // if INFO = i, i off-diagonal elements of an intermediate
638  // bidiagonal form did not converge to zero.
639  PetscBLASInt INFO = 0;
640 
641  // LAPACKgelss_(const PetscBLASInt *, // M
642  // const PetscBLASInt *, // N
643  // const PetscBLASInt *, // NRHS
644  // PetscScalar *, // A
645  // const PetscBLASInt *, // LDA
646  // PetscScalar *, // B
647  // const PetscBLASInt *, // LDB
648  // PetscReal *, // S(out) = singular values of A in increasing order
649  // const PetscReal *, // RCOND = tolerance for singular values
650  // PetscBLASInt *, // RANK(out) = number of "non-zero" singular values
651  // PetscScalar *, // WORK
652  // const PetscBLASInt *, // LWORK
653  // PetscBLASInt *); // INFO
654  LAPACKgelss_(&M, &N, &NRHS, pPS(A_trans_vals.data()), &LDA,
655  pPS(B.data()), &LDB, pPS(S.data()), &RCOND, &RANK,
656  pPS(WORK.data()), &LWORK, &INFO);
657 
658  // Check for errors in the Lapack call
659  libmesh_error_msg_if(INFO < 0, "Error, argument " << -INFO << " to LAPACKgelss_ had an illegal value.");
660  libmesh_error_msg_if(INFO > 0, "The algorithm for computing the SVD failed to converge!");
661 
662  // Debugging: print singular values and information about condition number:
663  // libMesh::err << "RCOND=" << RCOND << std::endl;
664  // libMesh::err << "Singular values: " << std::endl;
665  // for (auto i : index_range(S))
666  // libMesh::err << S[i] << std::endl;
667  // libMesh::err << "The condition number of A is approximately: " << S[0]/S.back() << std::endl;
668 
669  // Lapack has already written the solution into B, but it will be
670  // the wrong size for non-square problems, so we need to resize it
671  // correctly. The size of the solution vector should be the number
672  // of columns of the original A matrix. Unfortunately, resizing a
673  // DenseVector currently also zeros it out (unlike a std::vector) so
674  // we'll resize the underlying storage directly (the size is not
675  // stored independently elsewhere).
676  x.get_values().resize(this->n());
677 }
PetscScalar * pPS(T *ptr)
Definition: petsc_macro.h:174
Definition: assembly.h:38
PetscScalar PS(T val)
Definition: petsc_macro.h:168
void get_transpose(DenseMatrix< T > &dest) const
Put the tranposed matrix into dest.
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:140
unsigned int n() const

◆ add() [1/2]

template<typename T >
template<typename T2 , typename T3 >
boostcopy::enable_if_c< 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.

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

197 {
198  libmesh_assert_equal_to (this->m(), mat.m());
199  libmesh_assert_equal_to (this->n(), mat.n());
200 
201  for (auto j : make_range(this->n()))
202  for (auto i : make_range(this->m()))
203  this->el(i,j) += factor*mat.el(i,j);
204 }
unsigned int m() const
virtual T el(const unsigned int i, const unsigned int j) const =0
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:140
unsigned int n() const

◆ add() [2/2]

template<typename T >
template<typename T2 , typename T3 >
boostcopy::enable_if_c< 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.

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

1020 {
1021  libmesh_assert_equal_to (this->m(), mat.m());
1022  libmesh_assert_equal_to (this->n(), mat.n());
1023 
1024  for (auto i : make_range(this->m()))
1025  for (auto j : make_range(this->n()))
1026  (*this)(i,j) += factor * mat(i,j);
1027 }
unsigned int m() const
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:140
unsigned int n() const

◆ 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 841 of file dense_matrix_impl.h.

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

843 {
844  // Check for a previous decomposition
845  switch(this->_decomposition_type)
846  {
847  case NONE:
848  {
849  this->_cholesky_decompose ();
850  break;
851  }
852 
853  case CHOLESKY:
854  {
855  // Already factored, just need to call back_substitute.
856  break;
857  }
858 
859  default:
860  libmesh_error_msg("Error! This matrix already has a different decomposition...");
861  }
862 
863  // Perform back substitution
864  this->_cholesky_back_substitute (b, x);
865 }
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...
DecompositionType _decomposition_type
This flag keeps track of which type of decomposition has been performed on the matrix.
Definition: dense_matrix.h:650
void _cholesky_decompose()
Decomposes a symmetric positive definite matrix into a product of two lower triangular matrices accor...

◆ condense() [1/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.

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

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

77  {
78  libmesh_assert_equal_to (this->_m, rhs.size());
79  libmesh_assert_equal_to (iv, jv);
80 
81 
82  // move the known value into the RHS
83  // and zero the column
84  for (auto i : make_range(this->m()))
85  {
86  rhs.el(i) -= this->el(i,jv)*val;
87  this->el(i,jv) = 0.;
88  }
89 
90  // zero the row
91  for (auto j : make_range(this->n()))
92  this->el(iv,j) = 0.;
93 
94  this->el(iv,jv) = 1.;
95  rhs.el(iv) = val;
96 
97  }
unsigned int m() const
virtual T el(const unsigned int i, const unsigned int j) const =0
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:140
unsigned int n() const
unsigned int _m
The row dimension.

◆ condense() [2/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.

399  { DenseMatrixBase<T>::condense (i, j, val, rhs); }
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.

◆ 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 780 of file dense_matrix_impl.h.

781 {
782  switch(this->_decomposition_type)
783  {
784  case NONE:
785  {
786  // First LU decompose the matrix.
787  // Note that the lu_decompose routine will check to see if the
788  // matrix is square so we don't worry about it.
789  if (this->use_blas_lapack)
790  this->_lu_decompose_lapack();
791  else
792  this->_lu_decompose ();
793  }
794  case LU:
795  case LU_BLAS_LAPACK:
796  {
797  // Already decomposed, don't do anything
798  break;
799  }
800  default:
801  libmesh_error_msg("Error! Can't compute the determinant under the current decomposition.");
802  }
803 
804  // A variable to keep track of the running product of diagonal terms.
805  T determinant = 1.;
806 
807  // Loop over diagonal terms, computing the product. In practice,
808  // be careful because this value could easily become too large to
809  // fit in a double or float. To be safe, one should keep track of
810  // the power (of 10) of the determinant in a separate variable
811  // and maintain an order 1 value for the determinant itself.
812  unsigned int n_interchanges = 0;
813  for (auto i : make_range(this->m()))
814  {
815  if (this->_decomposition_type==LU)
816  if (_pivots[i] != static_cast<pivot_index_t>(i))
817  n_interchanges++;
818 
819  // Lapack pivots are 1-based!
821  if (_pivots[i] != static_cast<pivot_index_t>(i+1))
822  n_interchanges++;
823 
824  determinant *= (*this)(i,i);
825  }
826 
827  // Compute sign of determinant, depends on number of row interchanges!
828  // The sign should be (-1)^{n}, where n is the number of interchanges.
829  Real sign = n_interchanges % 2 == 0 ? 1. : -1.;
830 
831  return sign*determinant;
832 }
void _lu_decompose_lapack()
Computes an LU factorization of the matrix using the Lapack routine "getrf".
DecompositionType _decomposition_type
This flag keeps track of which type of decomposition has been performed on the matrix.
Definition: dense_matrix.h:650
unsigned int m() const
bool use_blas_lapack
Computes the inverse of the dense matrix (assuming it is invertible) by first computing the LU decomp...
Definition: dense_matrix.h:585
void _lu_decompose()
Form the LU decomposition of the matrix.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:140
std::vector< pivot_index_t > _pivots
Definition: dense_matrix.h:740

◆ 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.

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

38  {
39  DenseVector<T> ret(_m);
40  for (decltype(_m) i = 0; i < _m; ++i)
41  ret(i) = el(i, i);
42  return ret;
43  }
virtual T el(const unsigned int i, const unsigned int j) const =0
unsigned int _m
The row dimension.

◆ 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.

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().

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

◆ 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.

132  { 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 736 of file dense_matrix_impl.h.

738 {
739  // We use the LAPACK eigenvalue problem implementation
740  _evd_lapack(lambda_real, lambda_imag);
741 }
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".

◆ 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 746 of file dense_matrix_impl.h.

749 {
750  // We use the LAPACK eigenvalue problem implementation
751  _evd_lapack(lambda_real, lambda_imag, &VL, nullptr);
752 }
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".

◆ 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 768 of file dense_matrix_impl.h.

Referenced by DenseMatrixTest::testEVD_helper().

772 {
773  // We use the LAPACK eigenvalue problem implementation
774  _evd_lapack(lambda_real, lambda_imag, &VL, &VR);
775 }
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".

◆ 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 757 of file dense_matrix_impl.h.

760 {
761  // We use the LAPACK eigenvalue problem implementation
762  _evd_lapack(lambda_real, lambda_imag, nullptr, &VR);
763 }
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".

◆ get_principal_submatrix() [1/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.

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

473 {
474  libmesh_assert( (sub_m <= this->m()) && (sub_n <= this->n()) );
475 
476  dest.resize(sub_m, sub_n);
477  for (unsigned int i=0; i<sub_m; i++)
478  for (unsigned int j=0; j<sub_n; j++)
479  dest(i,j) = (*this)(i,j);
480 }
unsigned int m() const
libmesh_assert(ctx)
unsigned int n() const

◆ get_principal_submatrix() [2/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.

501 {
502  get_principal_submatrix(sub_m, sub_m, dest);
503 }
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.

◆ 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.

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

509 {
510  dest.resize(this->n(), this->m());
511 
512  for (auto i : make_range(dest.m()))
513  for (auto j : make_range(dest.n()))
514  dest(i,j) = (*this)(j,i);
515 }
unsigned int m() const
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:140
unsigned int n() const

◆ 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.

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

382 { return _val; }
std::vector< T > _val
The actual data values, stored as a 1D array.
Definition: dense_matrix.h:600

◆ 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.

387 { return _val; }
std::vector< T > _val
The actual data values, stored as a 1D array.
Definition: dense_matrix.h:600

◆ l1_norm()

template<typename T>
auto libMesh::DenseMatrix< T >::l1_norm ( ) const -> decltype(std::abs(T(0)))
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.

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

1126 {
1127  libmesh_assert (this->_m);
1128  libmesh_assert (this->_n);
1129 
1130  auto columnsum = std::abs(T(0));
1131  for (unsigned int i=0; i!=this->_m; i++)
1132  {
1133  columnsum += std::abs((*this)(i,0));
1134  }
1135  auto my_max = columnsum;
1136  for (unsigned int j=1; j!=this->_n; j++)
1137  {
1138  columnsum = 0.;
1139  for (unsigned int i=0; i!=this->_m; i++)
1140  {
1141  columnsum += std::abs((*this)(i,j));
1142  }
1143  my_max = (my_max > columnsum? my_max : columnsum);
1144  }
1145  return my_max;
1146 }
libmesh_assert(ctx)
unsigned int _n
The column dimension.
unsigned int _m
The row dimension.

◆ 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.

46 {
47  if (this->use_blas_lapack)
48  this->_multiply_blas(M2, LEFT_MULTIPLY);
49  else
50  {
51  // (*this) <- M2 * (*this)
52  // Where:
53  // (*this) = (m x n),
54  // M2 = (m x p),
55  // M3 = (p x n)
56 
57  // M3 is a copy of *this before it gets resize()d
58  DenseMatrix<T> M3(*this);
59 
60  // Resize *this so that the result can fit
61  this->resize (M2.m(), M3.n());
62 
63  // Call the multiply function in the base class
64  this->multiply(*this, M2, M3);
65  }
66 }
bool use_blas_lapack
Computes the inverse of the dense matrix (assuming it is invertible) by first computing the LU decomp...
Definition: dense_matrix.h:585
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 resize(const unsigned int new_m, const unsigned int new_n)
Resizes the matrix to the specified size and calls zero().
Definition: dense_matrix.h:895
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)...

◆ 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.

73 {
74  // (*this) <- M2 * (*this)
75  // Where:
76  // (*this) = (m x n),
77  // M2 = (m x p),
78  // M3 = (p x n)
79 
80  // M3 is a copy of *this before it gets resize()d
81  DenseMatrix<T> M3(*this);
82 
83  // Resize *this so that the result can fit
84  this->resize (M2.m(), M3.n());
85 
86  // Call the multiply function in the base class
87  this->multiply(*this, M2, M3);
88 }
void resize(const unsigned int new_m, const unsigned int new_n)
Resizes the matrix to the specified size and calls zero().
Definition: dense_matrix.h:895
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)...

◆ 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.

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

94 {
95  if (this->use_blas_lapack)
97  else
98  this->_left_multiply_transpose(A);
99 }
bool use_blas_lapack
Computes the inverse of the dense matrix (assuming it is invertible) by first computing the LU decomp...
Definition: dense_matrix.h:585
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 _left_multiply_transpose(const DenseMatrix< T2 > &A)
Left multiplies by the transpose of the matrix A which may contain a different numerical type...

◆ 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.

105 {
106  this->_left_multiply_transpose(A);
107 }
void _left_multiply_transpose(const DenseMatrix< T2 > &A)
Left multiplies by the transpose of the matrix A which may contain a different numerical type...

◆ linfty_norm()

template<typename T>
auto libMesh::DenseMatrix< T >::linfty_norm ( ) const -> decltype(std::abs(T(0)))
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 1152 of file dense_matrix.h.

1153 {
1154  libmesh_assert (this->_m);
1155  libmesh_assert (this->_n);
1156 
1157  auto rowsum = std::abs(T(0));
1158  for (unsigned int j=0; j!=this->_n; j++)
1159  {
1160  rowsum += std::abs((*this)(0,j));
1161  }
1162  auto my_max = rowsum;
1163  for (unsigned int i=1; i!=this->_m; i++)
1164  {
1165  rowsum = 0.;
1166  for (unsigned int j=0; j!=this->_n; j++)
1167  {
1168  rowsum += std::abs((*this)(i,j));
1169  }
1170  my_max = (my_max > rowsum? my_max : rowsum);
1171  }
1172  return my_max;
1173 }
libmesh_assert(ctx)
unsigned int _n
The column dimension.
unsigned int _m
The row dimension.

◆ 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.

Referenced by compute_enriched_soln(), libMesh::RBEIMConstruction::compute_max_eim_error(), fe_assembly(), libMesh::WeightedPatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::PatchRecoveryErrorEstimator::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().

523 {
524  // Check to be sure that the matrix is square before attempting
525  // an LU-solve. In general, one can compute the LU factorization of
526  // a non-square matrix, but:
527  //
528  // Overdetermined systems (m>n) have a solution only if enough of
529  // the equations are linearly-dependent.
530  //
531  // Underdetermined systems (m<n) typically have infinitely many
532  // solutions.
533  //
534  // We don't want to deal with either of these ambiguous cases here...
535  libmesh_assert_equal_to (this->m(), this->n());
536 
537  switch(this->_decomposition_type)
538  {
539  case NONE:
540  {
541  if (this->use_blas_lapack)
542  this->_lu_decompose_lapack();
543  else
544  this->_lu_decompose ();
545  break;
546  }
547 
548  case LU_BLAS_LAPACK:
549  {
550  // Already factored, just need to call back_substitute.
551  if (this->use_blas_lapack)
552  break;
553  }
554  libmesh_fallthrough();
555 
556  case LU:
557  {
558  // Already factored, just need to call back_substitute.
559  if (!(this->use_blas_lapack))
560  break;
561  }
562  libmesh_fallthrough();
563 
564  default:
565  libmesh_error_msg("Error! This matrix already has a different decomposition...");
566  }
567 
568  if (this->use_blas_lapack)
569  this->_lu_back_substitute_lapack (b, x);
570  else
571  this->_lu_back_substitute (b, x);
572 }
void _lu_back_substitute(const DenseVector< T > &b, DenseVector< T > &x) const
Solves the system Ax=b through back substitution.
void _lu_decompose_lapack()
Computes an LU factorization of the matrix using the Lapack routine "getrf".
DecompositionType _decomposition_type
This flag keeps track of which type of decomposition has been performed on the matrix.
Definition: dense_matrix.h:650
unsigned int m() const
bool use_blas_lapack
Computes the inverse of the dense matrix (assuming it is invertible) by first computing the LU decomp...
Definition: dense_matrix.h:585
void _lu_back_substitute_lapack(const DenseVector< T > &b, DenseVector< T > &x)
Companion function to _lu_decompose_lapack().
void _lu_decompose()
Form the LU decomposition of the matrix.
unsigned int n() const

◆ 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.

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

Referenced by libMesh::DenseMatrix< Real >::_left_multiply_transpose(), libMesh::DenseMatrix< Real >::_multiply_blas(), libMesh::DenseMatrix< Real >::_right_multiply_transpose(), libMesh::DenseMatrix< Real >::_svd_solve_lapack(), libMesh::DenseMatrixBase< T >::add(), libMesh::PetscMatrix< T >::add_block_matrix(), libMesh::SparseMatrix< ValOut >::add_block_matrix(), libMesh::StaticCondensation::add_matrix(), libMesh::EigenSparseMatrix< T >::add_matrix(), libMesh::LaspackMatrix< T >::add_matrix(), libMesh::DiagonalMatrix< T >::add_matrix(), libMesh::EpetraMatrix< T >::add_matrix(), libMesh::PetscMatrix< 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< Real >::get_transpose(), libMesh::DofMap::heterogeneously_constrain_element_jacobian_and_residual(), libMesh::DofMap::heterogeneously_constrain_element_residual(), libMesh::DenseMatrix< Real >::left_multiply(), libMesh::DofMap::max_constraint_error(), libMesh::DenseMatrixBase< T >::multiply(), libMesh::WeightedPatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::PatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::DenseMatrix< Real >::right_multiply(), libMesh::RBEIMEvaluation::set_interpolation_matrix_entry(), DualShapeTest::testEdge2Lagrange(), DenseMatrixTest::testEVD_helper(), DenseMatrixTest::testSVD(), and MetaPhysicL::RawType< libMesh::DenseMatrix< T > >::value().

104 { return _m; }
unsigned int _m
The row dimension.

◆ 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.

1105 {
1106  libmesh_assert (this->_m);
1107  libmesh_assert (this->_n);
1108  auto my_max = libmesh_real((*this)(0,0));
1109 
1110  for (unsigned int i=0; i!=this->_m; i++)
1111  {
1112  for (unsigned int j=0; j!=this->_n; j++)
1113  {
1114  auto current = libmesh_real((*this)(i,j));
1115  my_max = (my_max > current? my_max : current);
1116  }
1117  }
1118  return my_max;
1119 }
T libmesh_real(T a)
libmesh_assert(ctx)
unsigned int _n
The column dimension.
unsigned int _m
The row dimension.

◆ 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.

1084 {
1085  libmesh_assert (this->_m);
1086  libmesh_assert (this->_n);
1087  auto my_min = libmesh_real((*this)(0,0));
1088 
1089  for (unsigned int i=0; i!=this->_m; i++)
1090  {
1091  for (unsigned int j=0; j!=this->_n; j++)
1092  {
1093  auto current = libmesh_real((*this)(i,j));
1094  my_min = (my_min < current? my_min : current);
1095  }
1096  }
1097  return my_min;
1098 }
T libmesh_real(T a)
libmesh_assert(ctx)
unsigned int _n
The column dimension.
unsigned int _m
The row dimension.

◆ 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.

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

49  {
50  // Assertions to make sure we have been
51  // passed matrices of the correct dimension.
52  libmesh_assert_equal_to (M1.m(), M2.m());
53  libmesh_assert_equal_to (M1.n(), M3.n());
54  libmesh_assert_equal_to (M2.n(), M3.m());
55 
56  const unsigned int m_s = M2.m();
57  const unsigned int p_s = M2.n();
58  const unsigned int n_s = M1.n();
59 
60  // Do it this way because there is a
61  // decent chance (at least for constraint matrices)
62  // that M3(k,j) = 0. when right-multiplying.
63  for (unsigned int k=0; k<p_s; k++)
64  for (unsigned int j=0; j<n_s; j++)
65  if (M3.el(k,j) != 0.)
66  for (unsigned int i=0; i<m_s; i++)
67  M1.el(i,j) += M2.el(i,k) * M3.el(k,j);
68  }

◆ 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.

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

Referenced by libMesh::DenseMatrix< Real >::_left_multiply_transpose(), libMesh::DenseMatrix< Real >::_multiply_blas(), libMesh::DenseMatrix< Real >::_right_multiply_transpose(), libMesh::DenseMatrix< Real >::_svd_solve_lapack(), libMesh::DenseMatrixBase< T >::add(), libMesh::PetscMatrix< T >::add_block_matrix(), libMesh::SparseMatrix< ValOut >::add_block_matrix(), libMesh::StaticCondensation::add_matrix(), libMesh::EigenSparseMatrix< T >::add_matrix(), libMesh::LaspackMatrix< T >::add_matrix(), libMesh::DiagonalMatrix< T >::add_matrix(), libMesh::EpetraMatrix< T >::add_matrix(), libMesh::PetscMatrix< 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< Real >::get_transpose(), libMesh::DofMap::heterogeneously_constrain_element_jacobian_and_residual(), libMesh::DofMap::heterogeneously_constrain_element_residual(), libMesh::DenseMatrix< Real >::left_multiply(), libMesh::DofMap::max_constraint_error(), libMesh::DenseMatrixBase< T >::multiply(), libMesh::WeightedPatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::PatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::DenseMatrix< Real >::right_multiply(), libMesh::RBEIMEvaluation::set_interpolation_matrix_entry(), DualShapeTest::testEdge2Lagrange(), DenseMatrixTest::testSVD(), and MetaPhysicL::RawType< libMesh::DenseMatrix< T > >::value().

109 { return _n; }
unsigned int _n
The column dimension.

◆ 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.

1047 {
1048  for (auto i : index_range(_val))
1049  if (_val[i] != mat._val[i])
1050  return true;
1051 
1052  return false;
1053 }
std::vector< T > _val
The actual data values, stored as a 1D array.
Definition: dense_matrix.h:600
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition: int_range.h:117

◆ operator()() [1/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.

955 {
956  libmesh_assert_less (i*j, _val.size());
957  libmesh_assert_less (i, this->_m);
958  libmesh_assert_less (j, this->_n);
959 
960 
961  // return _val[(i) + (this->_m)*(j)]; // col-major
962  return _val[(i)*(this->_n) + (j)]; // row-major
963 }
unsigned int _n
The column dimension.
std::vector< T > _val
The actual data values, stored as a 1D array.
Definition: dense_matrix.h:600
unsigned int _m
The row dimension.

◆ operator()() [2/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.

971 {
972  libmesh_assert_less (i*j, _val.size());
973  libmesh_assert_less (i, this->_m);
974  libmesh_assert_less (j, this->_n);
975 
976  //return _val[(i) + (this->_m)*(j)]; // col-major
977  return _val[(i)*(this->_n) + (j)]; // row-major
978 }
unsigned int _n
The column dimension.
std::vector< T > _val
The actual data values, stored as a 1D array.
Definition: dense_matrix.h:600
unsigned int _m
The row dimension.

◆ 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.

1006 {
1007  this->scale(factor);
1008  return *this;
1009 }
void scale(const T factor)
Multiplies every element in the matrix by factor.
Definition: dense_matrix.h:986

◆ 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.

1060 {
1061  for (auto i : index_range(_val))
1062  _val[i] += mat._val[i];
1063 
1064  return *this;
1065 }
std::vector< T > _val
The actual data values, stored as a 1D array.
Definition: dense_matrix.h:600
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition: int_range.h:117

◆ 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.

1072 {
1073  for (auto i : index_range(_val))
1074  _val[i] -= mat._val[i];
1075 
1076  return *this;
1077 }
std::vector< T > _val
The actual data values, stored as a 1D array.
Definition: dense_matrix.h:600
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition: int_range.h:117

◆ operator=() [1/3]

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

◆ operator=() [2/3]

template<typename T>
DenseMatrix& libMesh::DenseMatrix< T >::operator= ( DenseMatrix< T > &&  )
default

◆ operator=() [3/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.

881 {
882  unsigned int mat_m = mat.m(), mat_n = mat.n();
883  this->resize(mat_m, mat_n);
884  for (unsigned int i=0; i<mat_m; i++)
885  for (unsigned int j=0; j<mat_n; j++)
886  (*this)(i,j) = mat(i,j);
887 
888  return *this;
889 }
void resize(const unsigned int new_m, const unsigned int new_n)
Resizes the matrix to the specified size and calls zero().
Definition: dense_matrix.h:895

◆ 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.

1034 {
1035  for (auto i : index_range(_val))
1036  if (_val[i] != mat._val[i])
1037  return false;
1038 
1039  return true;
1040 }
std::vector< T > _val
The actual data values, stored as a 1D array.
Definition: dense_matrix.h:600
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition: int_range.h:117

◆ 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]aVector whose entries correspond to rows in the product matrix.
[in]bVector whose entries correspond to columns in the product matrix.

Definition at line 485 of file dense_matrix_impl.h.

Referenced by DenseMatrixTest::testOuterProduct().

487 {
488  const unsigned int m = a.size();
489  const unsigned int n = b.size();
490 
491  this->resize(m, n);
492  for (unsigned int i = 0; i < m; ++i)
493  for (unsigned int j = 0; j < n; ++j)
494  (*this)(i,j) = a(i) * libmesh_conj(b(j));
495 }
T libmesh_conj(T a)
unsigned int m() const
void resize(const unsigned int new_m, const unsigned int new_n)
Resizes the matrix to the specified size and calls zero().
Definition: dense_matrix.h:895
unsigned int n() const

◆ 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.

References libMesh::make_range().

126  {
127  for (auto i : make_range(this->m()))
128  {
129  for (auto j : make_range(this->n()))
130  os << std::setw(8)
131  << this->el(i,j) << " ";
132 
133  os << std::endl;
134  }
135 
136  return;
137  }
unsigned int m() const
virtual T el(const unsigned int i, const unsigned int j) const =0
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:140
unsigned int n() const

◆ 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.

References libMesh::make_range().

102  {
103  // save the initial format flags
104  std::ios_base::fmtflags os_flags = os.flags();
105 
106  // Print the matrix entries.
107  for (auto i : make_range(this->m()))
108  {
109  for (auto j : make_range(this->n()))
110  os << std::setw(15)
111  << std::scientific
112  << std::setprecision(precision)
113  << this->el(i,j) << " ";
114 
115  os << std::endl;
116  }
117 
118  // reset the original format flags
119  os.flags(os_flags);
120  }
unsigned int m() const
virtual T el(const unsigned int i, const unsigned int j) const =0
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:140
unsigned int n() const

◆ 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.

Referenced by libMesh::DenseMatrix< Real >::_evd_lapack(), libMesh::DenseMatrix< Real >::_lu_decompose(), libMesh::DenseMatrix< Real >::_svd_lapack(), libMesh::HPCoarsenTest::add_projection(), alternative_fe_assembly(), assemble(), LinearElasticity::assemble(), assemble_1D(), AssembleOptimization::assemble_A_and_F(), assemble_biharmonic(), 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_SchroedingerEquation(), assemble_shell(), assemble_stokes(), 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< FEOutputType< T >::type >::coarsened_dof_values(), compute_enriched_soln(), compute_jacobian(), libMesh::FEGenericBase< FEOutputType< T >::type >::compute_periodic_constraints(), libMesh::FEGenericBase< FEOutputType< T >::type >::compute_proj_constraints(), compute_stresses(), fe_assembly(), form_matrixA(), NonlinearNeoHookeCurrentConfig::get_linearized_stiffness(), libMesh::DenseMatrix< Real >::get_principal_submatrix(), NonlinearNeoHookeCurrentConfig::get_residual(), libMesh::DenseMatrix< Real >::get_transpose(), libMesh::FESubdivision::init_subdivision_matrix(), libMesh::HDGProblem::jacobian(), LaplaceYoung::jacobian(), LargeDeformationElasticity::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(), Biharmonic::JR::residual_and_jacobian(), LinearElasticityWithContact::residual_and_jacobian(), libMesh::DenseMatrix< Real >::resize(), libMesh::TransientRBEvaluation::resize_data_structures(), libMesh::RBEvaluation::resize_data_structures(), libMesh::RBEIMEvaluation::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().

897 {
898  _val.resize(new_m*new_n);
899 
900  this->_m = new_m;
901  this->_n = new_n;
902 
903  // zero and set decomposition_type to NONE
904  this->zero();
905 }
virtual void zero() override final
Sets all elements of the matrix to 0 and resets any decomposition flag which may have been previously...
Definition: dense_matrix.h:911
unsigned int _n
The column dimension.
std::vector< T > _val
The actual data values, stored as a 1D array.
Definition: dense_matrix.h:600
unsigned int _m
The row dimension.

◆ 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.

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().

172 {
173  if (this->use_blas_lapack)
174  this->_multiply_blas(M3, RIGHT_MULTIPLY);
175  else
176  {
177  // (*this) <- (*this) * M3
178  // Where:
179  // (*this) = (m x n),
180  // M2 = (m x p),
181  // M3 = (p x n)
182 
183  // M2 is a copy of *this before it gets resize()d
184  DenseMatrix<T> M2(*this);
185 
186  // Resize *this so that the result can fit
187  this->resize (M2.m(), M3.n());
188 
189  this->multiply(*this, M2, M3);
190  }
191 }
bool use_blas_lapack
Computes the inverse of the dense matrix (assuming it is invertible) by first computing the LU decomp...
Definition: dense_matrix.h:585
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 resize(const unsigned int new_m, const unsigned int new_n)
Resizes the matrix to the specified size and calls zero().
Definition: dense_matrix.h:895
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)...

◆ 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.

198 {
199  // (*this) <- M3 * (*this)
200  // Where:
201  // (*this) = (m x n),
202  // M2 = (m x p),
203  // M3 = (p x n)
204 
205  // M2 is a copy of *this before it gets resize()d
206  DenseMatrix<T> M2(*this);
207 
208  // Resize *this so that the result can fit
209  this->resize (M2.m(), M3.n());
210 
211  this->multiply(*this, M2, M3);
212 }
void resize(const unsigned int new_m, const unsigned int new_n)
Resizes the matrix to the specified size and calls zero().
Definition: dense_matrix.h:895
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)...

◆ 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.

Referenced by NonlinearNeoHookeCurrentConfig::get_linearized_stiffness().

219 {
220  if (this->use_blas_lapack)
222  else
224 }
bool use_blas_lapack
Computes the inverse of the dense matrix (assuming it is invertible) by first computing the LU decomp...
Definition: dense_matrix.h:585
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. ...
Definition: assembly.h:38
void _right_multiply_transpose(const DenseMatrix< T2 > &A)
Right multiplies by the transpose of the matrix A which may contain a different numerical type...

◆ 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.

231 {
233 }
Definition: assembly.h:38
void _right_multiply_transpose(const DenseMatrix< T2 > &A)
Right multiplies by the transpose of the matrix A which may contain a different numerical type...

◆ 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.

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

987 {
988  for (auto & v : _val)
989  v *= factor;
990 }
std::vector< T > _val
The actual data values, stored as a 1D array.
Definition: dense_matrix.h:600

◆ 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.

996 {
997  for (auto i : make_range(this->m()))
998  (*this)(i, col) *= factor;
999 }
unsigned int m() const
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:140

◆ 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.

Referenced by DenseMatrixTest::testSubMatrix().

924 {
925  libmesh_assert_less (row_id + row_size - 1, this->_m);
926  libmesh_assert_less (col_id + col_size - 1, this->_n);
927 
928  DenseMatrix<T> sub;
929  sub._m = row_size;
930  sub._n = col_size;
931  sub._val.resize(row_size * col_size);
932 
933  unsigned int end_col = this->_n - col_size - col_id;
934  unsigned int p = row_id * this->_n;
935  unsigned int q = 0;
936  for (unsigned int i=0; i<row_size; i++)
937  {
938  // skip the beginning columns
939  p += col_id;
940  for (unsigned int j=0; j<col_size; j++)
941  sub._val[q++] = _val[p++];
942  // skip the rest columns
943  p += end_col;
944  }
945 
946  return sub;
947 }
unsigned int _n
The column dimension.
std::vector< T > _val
The actual data values, stored as a 1D array.
Definition: dense_matrix.h:600
unsigned int _m
The row dimension.

◆ 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 707 of file dense_matrix_impl.h.

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().

708 {
709  // We use the LAPACK svd implementation
710  _svd_lapack(sigma);
711 }
void _svd_lapack(DenseVector< Real > &sigma)
Computes an SVD of the matrix using the Lapack routine "getsvd".

◆ 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 715 of file dense_matrix_impl.h.

718 {
719  // We use the LAPACK svd implementation
720  _svd_lapack(sigma, U, VT);
721 }
void _svd_lapack(DenseVector< Real > &sigma)
Computes an SVD of the matrix using the Lapack routine "getsvd".

◆ 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 726 of file dense_matrix_impl.h.

Referenced by fe_assembly().

729 {
730  _svd_solve_lapack(rhs, x, rcond);
731 }
void _svd_solve_lapack(const DenseVector< T > &rhs, DenseVector< T > &x, Real rcond) const
Called by svd_solve(rhs).

◆ 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.

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

866 {
867  std::swap(this->_m, other_matrix._m);
868  std::swap(this->_n, other_matrix._n);
869  _val.swap(other_matrix._val);
871  _decomposition_type = other_matrix._decomposition_type;
872  other_matrix._decomposition_type = _temp;
873 }
DecompositionType _decomposition_type
This flag keeps track of which type of decomposition has been performed on the matrix.
Definition: dense_matrix.h:650
unsigned int _n
The column dimension.
std::vector< T > _val
The actual data values, stored as a 1D array.
Definition: dense_matrix.h:600
DecompositionType
The decomposition schemes above change the entries of the matrix A.
Definition: dense_matrix.h:644
unsigned int _m
The row dimension.

◆ 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 1179 of file dense_matrix.h.

Referenced by libMesh::DenseMatrix< Real >::_left_multiply_transpose().

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

◆ 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.

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

301 {
302  // Make sure the input sizes are compatible
303  libmesh_assert_equal_to (this->n(), arg.size());
304 
305  // Resize and clear dest.
306  // Note: DenseVector::resize() also zeros the vector.
307  dest.resize(this->m());
308 
309  // Short-circuit if the matrix is empty
310  if(this->m() == 0 || this->n() == 0)
311  return;
312 
313  if (this->use_blas_lapack)
314  this->_matvec_blas(1., 0., dest, arg);
315  else
316  {
317  const unsigned int n_rows = this->m();
318  const unsigned int n_cols = this->n();
319 
320  for (unsigned int i=0; i<n_rows; i++)
321  for (unsigned int j=0; j<n_cols; j++)
322  dest(i) += (*this)(i,j)*arg(j);
323  }
324 }
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.
unsigned int m() const
bool use_blas_lapack
Computes the inverse of the dense matrix (assuming it is invertible) by first computing the LU decomp...
Definition: dense_matrix.h:585
unsigned int n() const

◆ 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.

332 {
333  // Make sure the input sizes are compatible
334  libmesh_assert_equal_to (this->n(), arg.size());
335 
336  // Resize and clear dest.
337  // Note: DenseVector::resize() also zeros the vector.
338  dest.resize(this->m());
339 
340  // Short-circuit if the matrix is empty
341  if (this->m() == 0 || this->n() == 0)
342  return;
343 
344  const unsigned int n_rows = this->m();
345  const unsigned int n_cols = this->n();
346 
347  for (unsigned int i=0; i<n_rows; i++)
348  for (unsigned int j=0; j<n_cols; j++)
349  dest(i) += (*this)(i,j)*arg(j);
350 }
unsigned int m() const
unsigned int n() const

◆ 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.

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

428 {
429  // Short-circuit if the matrix is empty
430  if (this->m() == 0)
431  {
432  dest.resize(0);
433  return;
434  }
435 
436  if (this->use_blas_lapack)
437  this->_matvec_blas(factor, 1., dest, arg);
438  else
439  {
440  DenseVector<T> temp(arg.size());
441  this->vector_mult(temp, arg);
442  dest.add(factor, temp);
443  }
444 }
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.
unsigned int m() const
bool use_blas_lapack
Computes the inverse of the dense matrix (assuming it is invertible) by first computing the LU decomp...
Definition: dense_matrix.h:585
void vector_mult(DenseVector< T > &dest, const DenseVector< T > &arg) const
Performs the matrix-vector multiplication, dest := (*this) * arg.

◆ 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.

453 {
454  // Short-circuit if the matrix is empty
455  if (this->m() == 0)
456  {
457  dest.resize(0);
458  return;
459  }
460 
461  DenseVector<typename CompareTypes<T,T3>::supertype>
462  temp(arg.size());
463  this->vector_mult(temp, arg);
464  dest.add(factor, temp);
465 }
unsigned int m() const
void vector_mult(DenseVector< T > &dest, const DenseVector< T > &arg) const
Performs the matrix-vector multiplication, dest := (*this) * arg.

◆ 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.

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

357 {
358  // Make sure the input sizes are compatible
359  libmesh_assert_equal_to (this->m(), arg.size());
360 
361  // Resize and clear dest.
362  // Note: DenseVector::resize() also zeros the vector.
363  dest.resize(this->n());
364 
365  // Short-circuit if the matrix is empty
366  if (this->m() == 0)
367  return;
368 
369  if (this->use_blas_lapack)
370  {
371  this->_matvec_blas(1., 0., dest, arg, /*trans=*/true);
372  }
373  else
374  {
375  const unsigned int n_rows = this->m();
376  const unsigned int n_cols = this->n();
377 
378  // WORKS
379  // for (unsigned int j=0; j<n_cols; j++)
380  // for (unsigned int i=0; i<n_rows; i++)
381  // dest(j) += (*this)(i,j)*arg(i);
382 
383  // ALSO WORKS, (i,j) just swapped
384  for (unsigned int i=0; i<n_cols; i++)
385  for (unsigned int j=0; j<n_rows; j++)
386  dest(i) += (*this)(j,i)*arg(j);
387  }
388 }
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.
unsigned int m() const
bool use_blas_lapack
Computes the inverse of the dense matrix (assuming it is invertible) by first computing the LU decomp...
Definition: dense_matrix.h:585
unsigned int n() const

◆ 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.

396 {
397  // Make sure the input sizes are compatible
398  libmesh_assert_equal_to (this->m(), arg.size());
399 
400  // Resize and clear dest.
401  // Note: DenseVector::resize() also zeros the vector.
402  dest.resize(this->n());
403 
404  // Short-circuit if the matrix is empty
405  if (this->m() == 0)
406  return;
407 
408  const unsigned int n_rows = this->m();
409  const unsigned int n_cols = this->n();
410 
411  // WORKS
412  // for (unsigned int j=0; j<n_cols; j++)
413  // for (unsigned int i=0; i<n_rows; i++)
414  // dest(j) += (*this)(i,j)*arg(i);
415 
416  // ALSO WORKS, (i,j) just swapped
417  for (unsigned int i=0; i<n_cols; i++)
418  for (unsigned int j=0; j<n_rows; j++)
419  dest(i) += (*this)(j,i)*arg(j);
420 }
unsigned int m() const
unsigned int n() const

◆ 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.

Referenced by libMesh::HPCoarsenTest::add_projection(), libMesh::FEMSystem::assembly(), libMesh::FEGenericBase< FEOutputType< T >::type >::coarsened_dof_values(), libMesh::BoundaryProjectSolution::operator()(), libMesh::TransientRBEvaluation::rb_solve(), libMesh::RBEvaluation::rb_solve(), libMesh::HPCoarsenTest::select_refinement(), and DiagonalMatrixTest::testNumerics().

912 {
914 
915  std::fill (_val.begin(), _val.end(), static_cast<T>(0));
916 }
DecompositionType _decomposition_type
This flag keeps track of which type of decomposition has been performed on the matrix.
Definition: dense_matrix.h:650
std::vector< T > _val
The actual data values, stored as a 1D array.
Definition: dense_matrix.h:600

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.

◆ _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().

◆ _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().

◆ _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< Real >::get_values().

◆ 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: