libMesh
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]

Public Member Functions

 DenseMatrix (const unsigned int new_m=0, const unsigned int new_n=0)
 Constructor. 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
 Set every element in the matrix to 0. 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
 
virtual T & el (const unsigned int i, const unsigned int j) override
 
virtual void left_multiply (const DenseMatrixBase< T > &M2) override
 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
 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)
 Resize the matrix. 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...
 
Real min () const
 
Real max () const
 
Real l1_norm () const
 
Real linfty_norm () const
 
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...
 

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

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 73 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 668 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 670 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 588 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 576 of file dense_matrix.h.

Constructor & Destructor Documentation

◆ DenseMatrix() [1/3]

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

Constructor.

Creates a dense matrix of dimension m by n.

Definition at line 744 of file dense_matrix.h.

745  :
746  DenseMatrixBase<T>(new_m,new_n),
747 #if defined(LIBMESH_HAVE_PETSC) && defined(LIBMESH_USE_REAL_NUMBERS) && defined(LIBMESH_DEFAULT_DOUBLE_PRECISION)
748  use_blas_lapack(true),
749 #else
750  use_blas_lapack(false),
751 #endif
752  _val(),
754 {
755  this->resize(new_m,new_n);
756 }
DecompositionType _decomposition_type
This flag keeps track of which type of decomposition has been performed on the matrix.
Definition: dense_matrix.h:582
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:525
void resize(const unsigned int new_m, const unsigned int new_n)
Resize the matrix.
Definition: dense_matrix.h:792
std::vector< T > _val
The actual data values, stored as a 1D array.
Definition: dense_matrix.h:532

◆ DenseMatrix() [2/3]

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() [3/3]

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

1009 {
1010  // Shorthand notation for number of rows and columns.
1011  const unsigned int
1012  n_rows = this->m(),
1013  n_cols = this->n();
1014 
1015  // Just to be really sure...
1016  libmesh_assert_equal_to (n_rows, n_cols);
1017 
1018  // A convenient reference to *this
1019  const DenseMatrix<T> & A = *this;
1020 
1021  // Now compute the solution to Ax =b using the factorization.
1022  x.resize(n_rows);
1023 
1024  // Solve for Ly=b
1025  for (unsigned int i=0; i<n_cols; ++i)
1026  {
1027  T2 temp = b(i);
1028 
1029  for (unsigned int k=0; k<i; ++k)
1030  temp -= A(i,k)*x(k);
1031 
1032  x(i) = temp / A(i,i);
1033  }
1034 
1035  // Solve for L^T x = y
1036  for (unsigned int i=0; i<n_cols; ++i)
1037  {
1038  const unsigned int ib = (n_cols-1)-i;
1039 
1040  for (unsigned int k=(ib+1); k<n_cols; ++k)
1041  x(ib) -= A(k,ib) * x(k);
1042 
1043  x(ib) /= A(ib,ib);
1044  }
1045 }
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 961 of file dense_matrix_impl.h.

962 {
963  // If we called this function, there better not be any
964  // previous decomposition of the matrix.
965  libmesh_assert_equal_to (this->_decomposition_type, NONE);
966 
967  // Shorthand notation for number of rows and columns.
968  const unsigned int
969  n_rows = this->m(),
970  n_cols = this->n();
971 
972  // Just to be really sure...
973  libmesh_assert_equal_to (n_rows, n_cols);
974 
975  // A convenient reference to *this
976  DenseMatrix<T> & A = *this;
977 
978  for (unsigned int i=0; i<n_rows; ++i)
979  {
980  for (unsigned int j=i; j<n_cols; ++j)
981  {
982  for (unsigned int k=0; k<i; ++k)
983  A(i,j) -= A(i,k) * A(j,k);
984 
985  if (i == j)
986  {
987 #ifndef LIBMESH_USE_COMPLEX_NUMBERS
988  if (A(i,j) <= 0.0)
989  libmesh_error_msg("Error! Can only use Cholesky decomposition with symmetric positive definite matrices.");
990 #endif
991 
992  A(i,i) = std::sqrt(A(i,j));
993  }
994  else
995  A(j,i) = A(i,j) / A(i,i);
996  }
997  }
998 
999  // Set the flag for CHOLESKY decomposition
1000  this->_decomposition_type = CHOLESKY;
1001 }
DecompositionType _decomposition_type
This flag keeps track of which type of decomposition has been performed on the matrix.
Definition: dense_matrix.h:582
unsigned int m() const
unsigned int n() const

◆ _evd_lapack()

template<typename T>
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 ]

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

668 {
669  const unsigned int
670  n_cols = this->n();
671 
672  libmesh_assert_equal_to (this->m(), n_cols);
673  libmesh_assert_equal_to (this->m(), b.size());
674 
675  x.resize (n_cols);
676 
677  // A convenient reference to *this
678  const DenseMatrix<T> & A = *this;
679 
680  // Temporary vector storage. We use this instead of
681  // modifying the RHS.
682  DenseVector<T> z = b;
683 
684  // Lower-triangular "top to bottom" solve step, taking into account pivots
685  for (unsigned int i=0; i<n_cols; ++i)
686  {
687  // Swap
688  if (_pivots[i] != static_cast<pivot_index_t>(i))
689  std::swap( z(i), z(_pivots[i]) );
690 
691  x(i) = z(i);
692 
693  for (unsigned int j=0; j<i; ++j)
694  x(i) -= A(i,j)*x(j);
695 
696  x(i) /= A(i,i);
697  }
698 
699  // Upper-triangular "bottom to top" solve step
700  const unsigned int last_row = n_cols-1;
701 
702  for (int i=last_row; i>=0; --i)
703  {
704  for (int j=i+1; j<static_cast<int>(n_cols); ++j)
705  x(i) -= A(i,j)*x(j);
706  }
707 }
unsigned int m() const
void swap(Iterator &lhs, Iterator &rhs)
swap, used to implement op=
unsigned int n() const
std::vector< pivot_index_t > _pivots
Definition: dense_matrix.h:672

◆ _lu_back_substitute_lapack()

template<typename T>
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 ]

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

718 {
719  // If this function was called, there better not be any
720  // previous decomposition of the matrix.
721  libmesh_assert_equal_to (this->_decomposition_type, NONE);
722 
723  // Get the matrix size and make sure it is square
724  const unsigned int
725  n_rows = this->m();
726 
727  // A convenient reference to *this
728  DenseMatrix<T> & A = *this;
729 
730  _pivots.resize(n_rows);
731 
732  for (unsigned int i=0; i<n_rows; ++i)
733  {
734  // Find the pivot row by searching down the i'th column
735  _pivots[i] = i;
736 
737  // std::abs(complex) must return a Real!
738  Real the_max = std::abs( A(i,i) );
739  for (unsigned int j=i+1; j<n_rows; ++j)
740  {
741  Real candidate_max = std::abs( A(j,i) );
742  if (the_max < candidate_max)
743  {
744  the_max = candidate_max;
745  _pivots[i] = j;
746  }
747  }
748 
749  // libMesh::out << "the_max=" << the_max << " found at row " << _pivots[i] << std::endl;
750 
751  // If the max was found in a different row, interchange rows.
752  // Here we interchange the *entire* row, in Gaussian elimination
753  // you would only interchange the subrows A(i,j) and A(p(i),j), for j>i
754  if (_pivots[i] != static_cast<pivot_index_t>(i))
755  {
756  for (unsigned int j=0; j<n_rows; ++j)
757  std::swap( A(i,j), A(_pivots[i], j) );
758  }
759 
760 
761  // If the max abs entry found is zero, the matrix is singular
762  if (A(i,i) == libMesh::zero)
763  libmesh_error_msg("Matrix A is singular!");
764 
765  // Scale upper triangle entries of row i by the diagonal entry
766  // Note: don't scale the diagonal entry itself!
767  const T diag_inv = 1. / A(i,i);
768  for (unsigned int j=i+1; j<n_rows; ++j)
769  A(i,j) *= diag_inv;
770 
771  // Update the remaining sub-matrix A[i+1:m][i+1:m]
772  // by subtracting off (the diagonal-scaled)
773  // upper-triangular part of row i, scaled by the
774  // i'th column entry of each row. In terms of
775  // row operations, this is:
776  // for each r > i
777  // SubRow(r) = SubRow(r) - A(r,i)*SubRow(i)
778  //
779  // If we were scaling the i'th column as well, like
780  // in Gaussian elimination, this would 'zero' the
781  // entry in the i'th column.
782  for (unsigned int row=i+1; row<n_rows; ++row)
783  for (unsigned int col=i+1; col<n_rows; ++col)
784  A(row,col) -= A(row,i) * A(i,col);
785 
786  } // end i loop
787 
788  // Set the flag for LU decomposition
789  this->_decomposition_type = LU;
790 }
double abs(double a)
DecompositionType _decomposition_type
This flag keeps track of which type of decomposition has been performed on the matrix.
Definition: dense_matrix.h:582
unsigned int m() const
const Number zero
.
Definition: libmesh.h:239
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void swap(Iterator &lhs, Iterator &rhs)
swap, used to implement op=
std::vector< pivot_index_t > _pivots
Definition: dense_matrix.h:672

◆ _lu_decompose_lapack()

template<typename T>
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 ]

◆ _matvec_blas()

template<typename T>
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 ]

◆ _multiply_blas()

template<typename T>
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 ]

◆ _svd_helper()

template<typename T>
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 ]

◆ _svd_lapack() [1/2]

template<typename T>
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 ]

◆ _svd_lapack() [2/2]

template<typename T>
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 ]

◆ _svd_solve_lapack()

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

Called by svd_solve(rhs).

◆ 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 
)
inherited

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 188 of file dense_matrix_base.h.

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

190 {
191  libmesh_assert_equal_to (this->m(), mat.m());
192  libmesh_assert_equal_to (this->n(), mat.n());
193 
194  for (unsigned int j=0; j<this->n(); j++)
195  for (unsigned int i=0; i<this->m(); i++)
196  this->el(i,j) += factor*mat.el(i,j);
197 }
unsigned int m() const
virtual T el(const unsigned int i, const unsigned int j) const =0
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 
)

Adds factor times mat to this matrix.

Returns
A reference to *this.

Definition at line 884 of file dense_matrix.h.

886 {
887  libmesh_assert_equal_to (this->m(), mat.m());
888  libmesh_assert_equal_to (this->n(), mat.n());
889 
890  for (unsigned int i=0; i<this->m(); i++)
891  for (unsigned int j=0; j<this->n(); j++)
892  (*this)(i,j) += factor * mat(i,j);
893 }
unsigned int m() const
unsigned int n() const

◆ cholesky_solve()

template<typename T >
template<typename T2 >
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.

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

Definition at line 929 of file dense_matrix_impl.h.

Referenced by libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction >::operator()().

931 {
932  // Check for a previous decomposition
933  switch(this->_decomposition_type)
934  {
935  case NONE:
936  {
937  this->_cholesky_decompose ();
938  break;
939  }
940 
941  case CHOLESKY:
942  {
943  // Already factored, just need to call back_substitute.
944  break;
945  }
946 
947  default:
948  libmesh_error_msg("Error! This matrix already has a different decomposition...");
949  }
950 
951  // Perform back substitution
952  this->_cholesky_back_substitute (b, x);
953 }
DecompositionType _decomposition_type
This flag keeps track of which type of decomposition has been performed on the matrix.
Definition: dense_matrix.h:582
void _cholesky_back_substitute(const DenseVector< T2 > &b, DenseVector< T2 > &x) const
Solves the equation Ax=b for the unknown value x and rhs b based on the Cholesky factorization of A...
void _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.

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

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

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

358  { 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 868 of file dense_matrix_impl.h.

869 {
870  switch(this->_decomposition_type)
871  {
872  case NONE:
873  {
874  // First LU decompose the matrix.
875  // Note that the lu_decompose routine will check to see if the
876  // matrix is square so we don't worry about it.
877  if (this->use_blas_lapack)
878  this->_lu_decompose_lapack();
879  else
880  this->_lu_decompose ();
881  }
882  case LU:
883  case LU_BLAS_LAPACK:
884  {
885  // Already decomposed, don't do anything
886  break;
887  }
888  default:
889  libmesh_error_msg("Error! Can't compute the determinant under the current decomposition.");
890  }
891 
892  // A variable to keep track of the running product of diagonal terms.
893  T determinant = 1.;
894 
895  // Loop over diagonal terms, computing the product. In practice,
896  // be careful because this value could easily become too large to
897  // fit in a double or float. To be safe, one should keep track of
898  // the power (of 10) of the determinant in a separate variable
899  // and maintain an order 1 value for the determinant itself.
900  unsigned int n_interchanges = 0;
901  for (unsigned int i=0; i<this->m(); i++)
902  {
903  if (this->_decomposition_type==LU)
904  if (_pivots[i] != static_cast<pivot_index_t>(i))
905  n_interchanges++;
906 
907  // Lapack pivots are 1-based!
909  if (_pivots[i] != static_cast<pivot_index_t>(i+1))
910  n_interchanges++;
911 
912  determinant *= (*this)(i,i);
913  }
914 
915  // Compute sign of determinant, depends on number of row interchanges!
916  // The sign should be (-1)^{n}, where n is the number of interchanges.
917  Real sign = n_interchanges % 2 == 0 ? 1. : -1.;
918 
919  return sign*determinant;
920 }
DecompositionType _decomposition_type
This flag keeps track of which type of decomposition has been performed on the matrix.
Definition: dense_matrix.h:582
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:525
void _lu_decompose_lapack()
Computes an LU factorization of the matrix using the Lapack routine "getrf".
void _lu_decompose()
Form the LU decomposition of the matrix.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
std::vector< pivot_index_t > _pivots
Definition: dense_matrix.h:672

◆ el() [1/2]

template<typename T>
virtual T libMesh::DenseMatrix< T >::el ( const unsigned int  i,
const unsigned int  j 
) const
overridevirtual
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 88 of file dense_matrix.h.

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

◆ el() [2/2]

template<typename T>
virtual T& libMesh::DenseMatrix< T >::el ( const unsigned int  i,
const unsigned int  j 
)
overridevirtual
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 92 of file dense_matrix.h.

94  { 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 824 of file dense_matrix_impl.h.

826 {
827  // We use the LAPACK eigenvalue problem implementation
828  _evd_lapack(lambda_real, lambda_imag);
829 }
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 834 of file dense_matrix_impl.h.

837 {
838  // We use the LAPACK eigenvalue problem implementation
839  _evd_lapack(lambda_real, lambda_imag, &VL, nullptr);
840 }
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 856 of file dense_matrix_impl.h.

860 {
861  // We use the LAPACK eigenvalue problem implementation
862  _evd_lapack(lambda_real, lambda_imag, &VL, &VR);
863 }
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 845 of file dense_matrix_impl.h.

848 {
849  // We use the LAPACK eigenvalue problem implementation
850  _evd_lapack(lambda_real, lambda_imag, nullptr, &VR);
851 }
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 556 of file dense_matrix_impl.h.

559 {
560  libmesh_assert( (sub_m <= this->m()) && (sub_n <= this->n()) );
561 
562  dest.resize(sub_m, sub_n);
563  for (unsigned int i=0; i<sub_m; i++)
564  for (unsigned int j=0; j<sub_n; j++)
565  dest(i,j) = (*this)(i,j);
566 }
unsigned int m() const
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 586 of file dense_matrix_impl.h.

587 {
588  get_principal_submatrix(sub_m, sub_m, dest);
589 }
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 594 of file dense_matrix_impl.h.

595 {
596  dest.resize(this->n(), this->m());
597 
598  for (unsigned int i=0; i<dest.m(); i++)
599  for (unsigned int j=0; j<dest.n(); j++)
600  dest(i,j) = (*this)(j,i);
601 }
unsigned int m() const
unsigned int n() const

◆ get_values() [1/2]

template<typename T>
std::vector<T>& libMesh::DenseMatrix< T >::get_values ( )
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 341 of file dense_matrix.h.

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

◆ get_values() [2/2]

template<typename T>
const std::vector<T>& libMesh::DenseMatrix< T >::get_values ( ) const
Returns
A constant reference to the underlying data storage vector.

Definition at line 346 of file dense_matrix.h.

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

◆ l1_norm()

template<typename T >
Real libMesh::DenseMatrix< T >::l1_norm ( ) const
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 991 of file dense_matrix.h.

992 {
993  libmesh_assert (this->_m);
994  libmesh_assert (this->_n);
995 
996  Real columnsum = 0.;
997  for (unsigned int i=0; i!=this->_m; i++)
998  {
999  columnsum += std::abs((*this)(i,0));
1000  }
1001  Real my_max = columnsum;
1002  for (unsigned int j=1; j!=this->_n; j++)
1003  {
1004  columnsum = 0.;
1005  for (unsigned int i=0; i!=this->_m; i++)
1006  {
1007  columnsum += std::abs((*this)(i,j));
1008  }
1009  my_max = (my_max > columnsum? my_max : columnsum);
1010  }
1011  return my_max;
1012 }
double abs(double a)
unsigned int _n
The column dimension.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
unsigned int _m
The row dimension.

◆ left_multiply() [1/2]

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

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

Implements libMesh::DenseMatrixBase< T >.

Definition at line 40 of file dense_matrix_impl.h.

41 {
42  if (this->use_blas_lapack)
43  this->_multiply_blas(M2, LEFT_MULTIPLY);
44  else
45  {
46  // (*this) <- M2 * (*this)
47  // Where:
48  // (*this) = (m x n),
49  // M2 = (m x p),
50  // M3 = (p x n)
51 
52  // M3 is a copy of *this before it gets resize()d
53  DenseMatrix<T> M3(*this);
54 
55  // Resize *this so that the result can fit
56  this->resize (M2.m(), M3.n());
57 
58  // Call the multiply function in the base class
59  this->multiply(*this, M2, M3);
60  }
61 }
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)...
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. ...
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:525
void resize(const unsigned int new_m, const unsigned int new_n)
Resize the matrix.
Definition: dense_matrix.h:792

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

68 {
69  // (*this) <- M2 * (*this)
70  // Where:
71  // (*this) = (m x n),
72  // M2 = (m x p),
73  // M3 = (p x n)
74 
75  // M3 is a copy of *this before it gets resize()d
76  DenseMatrix<T> M3(*this);
77 
78  // Resize *this so that the result can fit
79  this->resize (M2.m(), M3.n());
80 
81  // Call the multiply function in the base class
82  this->multiply(*this, M2, M3);
83 }
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)...
void resize(const unsigned int new_m, const unsigned int new_n)
Resize the matrix.
Definition: dense_matrix.h:792

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

89 {
90  if (this->use_blas_lapack)
92  else
93  {
94  //Check to see if we are doing (A^T)*A
95  if (this == &A)
96  {
97  //libmesh_here();
98  DenseMatrix<T> B(*this);
99 
100  // Simple but inefficient way
101  // return this->left_multiply_transpose(B);
102 
103  // More efficient, but more code way
104  // If A is mxn, the result will be a square matrix of Size n x n.
105  const unsigned int n_rows = A.m();
106  const unsigned int n_cols = A.n();
107 
108  // resize() *this and also zero out all entries.
109  this->resize(n_cols,n_cols);
110 
111  // Compute the lower-triangular part
112  for (unsigned int i=0; i<n_cols; ++i)
113  for (unsigned int j=0; j<=i; ++j)
114  for (unsigned int k=0; k<n_rows; ++k) // inner products are over n_rows
115  (*this)(i,j) += B(k,i)*B(k,j);
116 
117  // Copy lower-triangular part into upper-triangular part
118  for (unsigned int i=0; i<n_cols; ++i)
119  for (unsigned int j=i+1; j<n_cols; ++j)
120  (*this)(i,j) = (*this)(j,i);
121  }
122 
123  else
124  {
125  DenseMatrix<T> B(*this);
126 
127  this->resize (A.n(), B.n());
128 
129  libmesh_assert_equal_to (A.m(), B.m());
130  libmesh_assert_equal_to (this->m(), A.n());
131  libmesh_assert_equal_to (this->n(), B.n());
132 
133  const unsigned int m_s = A.n();
134  const unsigned int p_s = A.m();
135  const unsigned int n_s = this->n();
136 
137  // Do it this way because there is a
138  // decent chance (at least for constraint matrices)
139  // that A.transpose(i,k) = 0.
140  for (unsigned int i=0; i<m_s; i++)
141  for (unsigned int k=0; k<p_s; k++)
142  if (A.transpose(i,k) != 0.)
143  for (unsigned int j=0; j<n_s; j++)
144  (*this)(i,j) += A.transpose(i,k)*B(k,j);
145  }
146  }
147 
148 }
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. ...
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:525
Definition: assembly.h:38
void resize(const unsigned int new_m, const unsigned int new_n)
Resize the matrix.
Definition: dense_matrix.h:792
unsigned int n() const

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

155 {
156  //Check to see if we are doing (A^T)*A
157  if (this == &A)
158  {
159  //libmesh_here();
160  DenseMatrix<T> B(*this);
161 
162  // Simple but inefficient way
163  // return this->left_multiply_transpose(B);
164 
165  // More efficient, but more code way
166  // If A is mxn, the result will be a square matrix of Size n x n.
167  const unsigned int n_rows = A.m();
168  const unsigned int n_cols = A.n();
169 
170  // resize() *this and also zero out all entries.
171  this->resize(n_cols,n_cols);
172 
173  // Compute the lower-triangular part
174  for (unsigned int i=0; i<n_cols; ++i)
175  for (unsigned int j=0; j<=i; ++j)
176  for (unsigned int k=0; k<n_rows; ++k) // inner products are over n_rows
177  (*this)(i,j) += B(k,i)*B(k,j);
178 
179  // Copy lower-triangular part into upper-triangular part
180  for (unsigned int i=0; i<n_cols; ++i)
181  for (unsigned int j=i+1; j<n_cols; ++j)
182  (*this)(i,j) = (*this)(j,i);
183  }
184 
185  else
186  {
187  DenseMatrix<T> B(*this);
188 
189  this->resize (A.n(), B.n());
190 
191  libmesh_assert_equal_to (A.m(), B.m());
192  libmesh_assert_equal_to (this->m(), A.n());
193  libmesh_assert_equal_to (this->n(), B.n());
194 
195  const unsigned int m_s = A.n();
196  const unsigned int p_s = A.m();
197  const unsigned int n_s = this->n();
198 
199  // Do it this way because there is a
200  // decent chance (at least for constraint matrices)
201  // that A.transpose(i,k) = 0.
202  for (unsigned int i=0; i<m_s; i++)
203  for (unsigned int k=0; k<p_s; k++)
204  if (A.transpose(i,k) != 0.)
205  for (unsigned int j=0; j<n_s; j++)
206  (*this)(i,j) += A.transpose(i,k)*B(k,j);
207  }
208 }
unsigned int m() const
Definition: assembly.h:38
void resize(const unsigned int new_m, const unsigned int new_n)
Resize the matrix.
Definition: dense_matrix.h:792
unsigned int n() const

◆ linfty_norm()

template<typename T >
Real libMesh::DenseMatrix< T >::linfty_norm ( ) const
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 1018 of file dense_matrix.h.

1019 {
1020  libmesh_assert (this->_m);
1021  libmesh_assert (this->_n);
1022 
1023  Real rowsum = 0.;
1024  for (unsigned int j=0; j!=this->_n; j++)
1025  {
1026  rowsum += std::abs((*this)(0,j));
1027  }
1028  Real my_max = rowsum;
1029  for (unsigned int i=1; i!=this->_m; i++)
1030  {
1031  rowsum = 0.;
1032  for (unsigned int j=0; j!=this->_n; j++)
1033  {
1034  rowsum += std::abs((*this)(i,j));
1035  }
1036  my_max = (my_max > rowsum? my_max : rowsum);
1037  }
1038  return my_max;
1039 }
double abs(double a)
unsigned int _n
The column dimension.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
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.

Definition at line 607 of file dense_matrix_impl.h.

609 {
610  // Check to be sure that the matrix is square before attempting
611  // an LU-solve. In general, one can compute the LU factorization of
612  // a non-square matrix, but:
613  //
614  // Overdetermined systems (m>n) have a solution only if enough of
615  // the equations are linearly-dependent.
616  //
617  // Underdetermined systems (m<n) typically have infinitely many
618  // solutions.
619  //
620  // We don't want to deal with either of these ambiguous cases here...
621  libmesh_assert_equal_to (this->m(), this->n());
622 
623  switch(this->_decomposition_type)
624  {
625  case NONE:
626  {
627  if (this->use_blas_lapack)
628  this->_lu_decompose_lapack();
629  else
630  this->_lu_decompose ();
631  break;
632  }
633 
634  case LU_BLAS_LAPACK:
635  {
636  // Already factored, just need to call back_substitute.
637  if (this->use_blas_lapack)
638  break;
639  }
640  libmesh_fallthrough();
641 
642  case LU:
643  {
644  // Already factored, just need to call back_substitute.
645  if (!(this->use_blas_lapack))
646  break;
647  }
648  libmesh_fallthrough();
649 
650  default:
651  libmesh_error_msg("Error! This matrix already has a different decomposition...");
652  }
653 
654  if (this->use_blas_lapack)
655  this->_lu_back_substitute_lapack (b, x);
656  else
657  this->_lu_back_substitute (b, x);
658 }
void _lu_back_substitute(const DenseVector< T > &b, DenseVector< T > &x) const
Solves the system Ax=b through back substitution.
DecompositionType _decomposition_type
This flag keeps track of which type of decomposition has been performed on the matrix.
Definition: dense_matrix.h:582
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:525
void _lu_back_substitute_lapack(const DenseVector< T > &b, DenseVector< T > &x)
Companion function to _lu_decompose_lapack().
void _lu_decompose_lapack()
Computes an LU factorization of the matrix using the Lapack routine "getrf".
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
inherited

◆ max()

template<typename T >
Real libMesh::DenseMatrix< T >::max ( ) const
Returns
The maximum entry in the matrix, or the maximum real part in the case of complex numbers.

Definition at line 970 of file dense_matrix.h.

971 {
972  libmesh_assert (this->_m);
973  libmesh_assert (this->_n);
974  Real my_max = libmesh_real((*this)(0,0));
975 
976  for (unsigned int i=0; i!=this->_m; i++)
977  {
978  for (unsigned int j=0; j!=this->_n; j++)
979  {
980  Real current = libmesh_real((*this)(i,j));
981  my_max = (my_max > current? my_max : current);
982  }
983  }
984  return my_max;
985 }
T libmesh_real(T a)
unsigned int _n
The column dimension.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
unsigned int _m
The row dimension.

◆ min()

template<typename T >
Real libMesh::DenseMatrix< T >::min ( ) const
Returns
The minimum entry in the matrix, or the minimum real part in the case of complex numbers.

Definition at line 949 of file dense_matrix.h.

950 {
951  libmesh_assert (this->_m);
952  libmesh_assert (this->_n);
953  Real my_min = libmesh_real((*this)(0,0));
954 
955  for (unsigned int i=0; i!=this->_m; i++)
956  {
957  for (unsigned int j=0; j!=this->_n; j++)
958  {
959  Real current = libmesh_real((*this)(i,j));
960  my_min = (my_min < current? my_min : current);
961  }
962  }
963  return my_min;
964 }
T libmesh_real(T a)
unsigned int _n
The column dimension.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
unsigned int _m
The row dimension.

◆ multiply()

template<typename T>
static 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)

◆ n()

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

◆ operator!=()

template<typename T>
bool libMesh::DenseMatrix< T >::operator!= ( const DenseMatrix< T > &  mat) const
Returns
true if mat is not exactly equal to this matrix, false otherwise.

Definition at line 912 of file dense_matrix.h.

913 {
914  for (std::size_t i=0; i<_val.size(); i++)
915  if (_val[i] != mat._val[i])
916  return true;
917 
918  return false;
919 }
std::vector< T > _val
The actual data values, stored as a 1D array.
Definition: dense_matrix.h:532

◆ operator()() [1/2]

template<typename T >
T libMesh::DenseMatrix< T >::operator() ( const unsigned int  i,
const unsigned int  j 
) const
Returns
The (i,j) element of the matrix.

Definition at line 819 of file dense_matrix.h.

821 {
822  libmesh_assert_less (i*j, _val.size());
823  libmesh_assert_less (i, this->_m);
824  libmesh_assert_less (j, this->_n);
825 
826 
827  // return _val[(i) + (this->_m)*(j)]; // col-major
828  return _val[(i)*(this->_n) + (j)]; // row-major
829 }
unsigned int _n
The column dimension.
std::vector< T > _val
The actual data values, stored as a 1D array.
Definition: dense_matrix.h:532
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 
)
Returns
The (i,j) element of the matrix as a writable reference.

Definition at line 835 of file dense_matrix.h.

837 {
838  libmesh_assert_less (i*j, _val.size());
839  libmesh_assert_less (i, this->_m);
840  libmesh_assert_less (j, this->_n);
841 
842  //return _val[(i) + (this->_m)*(j)]; // col-major
843  return _val[(i)*(this->_n) + (j)]; // row-major
844 }
unsigned int _n
The column dimension.
std::vector< T > _val
The actual data values, stored as a 1D array.
Definition: dense_matrix.h:532
unsigned int _m
The row dimension.

◆ operator*=()

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

Multiplies every element in the matrix by factor.

Returns
A reference to *this.

Definition at line 871 of file dense_matrix.h.

872 {
873  this->scale(factor);
874  return *this;
875 }
void scale(const T factor)
Multiplies every element in the matrix by factor.
Definition: dense_matrix.h:852

◆ operator+=()

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

Adds mat to this matrix.

Returns
A reference to *this.

Definition at line 925 of file dense_matrix.h.

926 {
927  for (std::size_t i=0; i<_val.size(); i++)
928  _val[i] += mat._val[i];
929 
930  return *this;
931 }
std::vector< T > _val
The actual data values, stored as a 1D array.
Definition: dense_matrix.h:532

◆ operator-=()

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

Subtracts mat from this matrix.

Returns
A reference to *this.

Definition at line 937 of file dense_matrix.h.

938 {
939  for (std::size_t i=0; i<_val.size(); i++)
940  _val[i] -= mat._val[i];
941 
942  return *this;
943 }
std::vector< T > _val
The actual data values, stored as a 1D array.
Definition: dense_matrix.h:532

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

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

778 {
779  unsigned int mat_m = mat.m(), mat_n = mat.n();
780  this->resize(mat_m, mat_n);
781  for (unsigned int i=0; i<mat_m; i++)
782  for (unsigned int j=0; j<mat_n; j++)
783  (*this)(i,j) = mat(i,j);
784 
785  return *this;
786 }
void resize(const unsigned int new_m, const unsigned int new_n)
Resize the matrix.
Definition: dense_matrix.h:792

◆ operator==()

template<typename T>
bool libMesh::DenseMatrix< T >::operator== ( const DenseMatrix< T > &  mat) const
Returns
true if mat is exactly equal to this matrix, false otherwise.

Definition at line 899 of file dense_matrix.h.

900 {
901  for (std::size_t i=0; i<_val.size(); i++)
902  if (_val[i] != mat._val[i])
903  return false;
904 
905  return true;
906 }
std::vector< T > _val
The actual data values, stored as a 1D array.
Definition: dense_matrix.h:532

◆ 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 ${a}$ and ${b}$ is

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

The outer product of two complex-valued vectors ${a}$ and ${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 571 of file dense_matrix_impl.h.

573 {
574  const unsigned int m = a.size();
575  const unsigned int n = b.size();
576 
577  this->resize(m, n);
578  for (unsigned int i = 0; i < m; ++i)
579  for (unsigned int j = 0; j < n; ++j)
580  (*this)(i,j) = a(i) * libmesh_conj(b(j));
581 }
T libmesh_conj(T a)
unsigned int m() const
void resize(const unsigned int new_m, const unsigned int new_n)
Resize the matrix.
Definition: dense_matrix.h:792
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.

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

◆ resize()

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

Resize the matrix.

Will never free memory, but may allocate more. Sets all elements to 0.

Definition at line 792 of file dense_matrix.h.

Referenced by libMesh::DenseMatrix< Real >::_lu_decompose(), libMesh::DenseMatrix< Real >::get_principal_submatrix(), libMesh::DenseMatrix< Real >::get_transpose(), libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction >::operator()(), and libMesh::DenseMatrix< Real >::resize().

794 {
795  _val.resize(new_m*new_n);
796 
797  this->_m = new_m;
798  this->_n = new_n;
799 
800  // zero and set decomposition_type to NONE
801  this->zero();
802 }
unsigned int _n
The column dimension.
virtual void zero() override
Set every element in the matrix to 0.
Definition: dense_matrix.h:808
std::vector< T > _val
The actual data values, stored as a 1D array.
Definition: dense_matrix.h:532
unsigned int _m
The row dimension.

◆ right_multiply() [1/2]

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

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

Implements libMesh::DenseMatrixBase< T >.

Definition at line 213 of file dense_matrix_impl.h.

214 {
215  if (this->use_blas_lapack)
216  this->_multiply_blas(M3, RIGHT_MULTIPLY);
217  else
218  {
219  // (*this) <- M3 * (*this)
220  // Where:
221  // (*this) = (m x n),
222  // M2 = (m x p),
223  // M3 = (p x n)
224 
225  // M2 is a copy of *this before it gets resize()d
226  DenseMatrix<T> M2(*this);
227 
228  // Resize *this so that the result can fit
229  this->resize (M2.m(), M3.n());
230 
231  this->multiply(*this, M2, M3);
232  }
233 }
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)...
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. ...
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:525
void resize(const unsigned int new_m, const unsigned int new_n)
Resize the matrix.
Definition: dense_matrix.h:792

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

240 {
241  // (*this) <- M3 * (*this)
242  // Where:
243  // (*this) = (m x n),
244  // M2 = (m x p),
245  // M3 = (p x n)
246 
247  // M2 is a copy of *this before it gets resize()d
248  DenseMatrix<T> M2(*this);
249 
250  // Resize *this so that the result can fit
251  this->resize (M2.m(), M3.n());
252 
253  this->multiply(*this, M2, M3);
254 }
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)...
void resize(const unsigned int new_m, const unsigned int new_n)
Resize the matrix.
Definition: dense_matrix.h:792

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

261 {
262  if (this->use_blas_lapack)
264  else
265  {
266  //Check to see if we are doing B*(B^T)
267  if (this == &B)
268  {
269  //libmesh_here();
270  DenseMatrix<T> A(*this);
271 
272  // Simple but inefficient way
273  // return this->right_multiply_transpose(A);
274 
275  // More efficient, more code
276  // If B is mxn, the result will be a square matrix of Size m x m.
277  const unsigned int n_rows = B.m();
278  const unsigned int n_cols = B.n();
279 
280  // resize() *this and also zero out all entries.
281  this->resize(n_rows,n_rows);
282 
283  // Compute the lower-triangular part
284  for (unsigned int i=0; i<n_rows; ++i)
285  for (unsigned int j=0; j<=i; ++j)
286  for (unsigned int k=0; k<n_cols; ++k) // inner products are over n_cols
287  (*this)(i,j) += A(i,k)*A(j,k);
288 
289  // Copy lower-triangular part into upper-triangular part
290  for (unsigned int i=0; i<n_rows; ++i)
291  for (unsigned int j=i+1; j<n_rows; ++j)
292  (*this)(i,j) = (*this)(j,i);
293  }
294 
295  else
296  {
297  DenseMatrix<T> A(*this);
298 
299  this->resize (A.m(), B.m());
300 
301  libmesh_assert_equal_to (A.n(), B.n());
302  libmesh_assert_equal_to (this->m(), A.m());
303  libmesh_assert_equal_to (this->n(), B.m());
304 
305  const unsigned int m_s = A.m();
306  const unsigned int p_s = A.n();
307  const unsigned int n_s = this->n();
308 
309  // Do it this way because there is a
310  // decent chance (at least for constraint matrices)
311  // that B.transpose(k,j) = 0.
312  for (unsigned int j=0; j<n_s; j++)
313  for (unsigned int k=0; k<p_s; k++)
314  if (B.transpose(k,j) != 0.)
315  for (unsigned int i=0; i<m_s; i++)
316  (*this)(i,j) += A(i,k)*B.transpose(k,j);
317  }
318  }
319 }
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. ...
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:525
Definition: assembly.h:38
void resize(const unsigned int new_m, const unsigned int new_n)
Resize the matrix.
Definition: dense_matrix.h:792
unsigned int n() const

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

326 {
327  //Check to see if we are doing B*(B^T)
328  if (this == &B)
329  {
330  //libmesh_here();
331  DenseMatrix<T> A(*this);
332 
333  // Simple but inefficient way
334  // return this->right_multiply_transpose(A);
335 
336  // More efficient, more code
337  // If B is mxn, the result will be a square matrix of Size m x m.
338  const unsigned int n_rows = B.m();
339  const unsigned int n_cols = B.n();
340 
341  // resize() *this and also zero out all entries.
342  this->resize(n_rows,n_rows);
343 
344  // Compute the lower-triangular part
345  for (unsigned int i=0; i<n_rows; ++i)
346  for (unsigned int j=0; j<=i; ++j)
347  for (unsigned int k=0; k<n_cols; ++k) // inner products are over n_cols
348  (*this)(i,j) += A(i,k)*A(j,k);
349 
350  // Copy lower-triangular part into upper-triangular part
351  for (unsigned int i=0; i<n_rows; ++i)
352  for (unsigned int j=i+1; j<n_rows; ++j)
353  (*this)(i,j) = (*this)(j,i);
354  }
355 
356  else
357  {
358  DenseMatrix<T> A(*this);
359 
360  this->resize (A.m(), B.m());
361 
362  libmesh_assert_equal_to (A.n(), B.n());
363  libmesh_assert_equal_to (this->m(), A.m());
364  libmesh_assert_equal_to (this->n(), B.m());
365 
366  const unsigned int m_s = A.m();
367  const unsigned int p_s = A.n();
368  const unsigned int n_s = this->n();
369 
370  // Do it this way because there is a
371  // decent chance (at least for constraint matrices)
372  // that B.transpose(k,j) = 0.
373  for (unsigned int j=0; j<n_s; j++)
374  for (unsigned int k=0; k<p_s; k++)
375  if (B.transpose(k,j) != 0.)
376  for (unsigned int i=0; i<m_s; i++)
377  (*this)(i,j) += A(i,k)*B.transpose(k,j);
378  }
379 }
unsigned int m() const
Definition: assembly.h:38
void resize(const unsigned int new_m, const unsigned int new_n)
Resize the matrix.
Definition: dense_matrix.h:792
unsigned int n() const

◆ scale()

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

Multiplies every element in the matrix by factor.

Definition at line 852 of file dense_matrix.h.

853 {
854  for (std::size_t i=0; i<_val.size(); i++)
855  _val[i] *= factor;
856 }
std::vector< T > _val
The actual data values, stored as a 1D array.
Definition: dense_matrix.h:532

◆ scale_column()

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

Multiplies every element in the column col matrix by factor.

Definition at line 861 of file dense_matrix.h.

862 {
863  for (unsigned int i=0; i<this->m(); i++)
864  (*this)(i, col) *= factor;
865 }
unsigned int m() const

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

796 {
797  // We use the LAPACK svd implementation
798  _svd_lapack(sigma);
799 }
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 803 of file dense_matrix_impl.h.

806 {
807  // We use the LAPACK svd implementation
808  _svd_lapack(sigma, U, VT);
809 }
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 814 of file dense_matrix_impl.h.

817 {
818  _svd_solve_lapack(rhs, x, rcond);
819 }
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)

STL-like swap method.

Definition at line 762 of file dense_matrix.h.

763 {
764  std::swap(this->_m, other_matrix._m);
765  std::swap(this->_n, other_matrix._n);
766  _val.swap(other_matrix._val);
768  _decomposition_type = other_matrix._decomposition_type;
769  other_matrix._decomposition_type = _temp;
770 }
DecompositionType _decomposition_type
This flag keeps track of which type of decomposition has been performed on the matrix.
Definition: dense_matrix.h:582
unsigned int _n
The column dimension.
void swap(Iterator &lhs, Iterator &rhs)
swap, used to implement op=
std::vector< T > _val
The actual data values, stored as a 1D array.
Definition: dense_matrix.h:532
DecompositionType
The decomposition schemes above change the entries of the matrix A.
Definition: dense_matrix.h:576
unsigned int _m
The row dimension.

◆ transpose()

template<typename T >
T libMesh::DenseMatrix< T >::transpose ( const unsigned int  i,
const unsigned int  j 
) const
Returns
The (i,j) element of the transposed matrix.

Definition at line 1045 of file dense_matrix.h.

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

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

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

387 {
388  // Make sure the input sizes are compatible
389  libmesh_assert_equal_to (this->n(), arg.size());
390 
391  // Resize and clear dest.
392  // Note: DenseVector::resize() also zeros the vector.
393  dest.resize(this->m());
394 
395  // Short-circuit if the matrix is empty
396  if(this->m() == 0 || this->n() == 0)
397  return;
398 
399  if (this->use_blas_lapack)
400  this->_matvec_blas(1., 0., dest, arg);
401  else
402  {
403  const unsigned int n_rows = this->m();
404  const unsigned int n_cols = this->n();
405 
406  for (unsigned int i=0; i<n_rows; i++)
407  for (unsigned int j=0; j<n_cols; j++)
408  dest(i) += (*this)(i,j)*arg(j);
409  }
410 }
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:525
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 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 416 of file dense_matrix_impl.h.

418 {
419  // Make sure the input sizes are compatible
420  libmesh_assert_equal_to (this->n(), arg.size());
421 
422  // Resize and clear dest.
423  // Note: DenseVector::resize() also zeros the vector.
424  dest.resize(this->m());
425 
426  // Short-circuit if the matrix is empty
427  if (this->m() == 0 || this->n() == 0)
428  return;
429 
430  const unsigned int n_rows = this->m();
431  const unsigned int n_cols = this->n();
432 
433  for (unsigned int i=0; i<n_rows; i++)
434  for (unsigned int j=0; j<n_cols; j++)
435  dest(i) += (*this)(i,j)*arg(j);
436 }
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 511 of file dense_matrix_impl.h.

514 {
515  // Short-circuit if the matrix is empty
516  if (this->m() == 0)
517  {
518  dest.resize(0);
519  return;
520  }
521 
522  if (this->use_blas_lapack)
523  this->_matvec_blas(factor, 1., dest, arg);
524  else
525  {
526  DenseVector<T> temp(arg.size());
527  this->vector_mult(temp, arg);
528  dest.add(factor, temp);
529  }
530 }
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:525
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.
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 536 of file dense_matrix_impl.h.

539 {
540  // Short-circuit if the matrix is empty
541  if (this->m() == 0)
542  {
543  dest.resize(0);
544  return;
545  }
546 
547  DenseVector<typename CompareTypes<T,T3>::supertype>
548  temp(arg.size());
549  this->vector_mult(temp, arg);
550  dest.add(factor, temp);
551 }
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 441 of file dense_matrix_impl.h.

443 {
444  // Make sure the input sizes are compatible
445  libmesh_assert_equal_to (this->m(), arg.size());
446 
447  // Resize and clear dest.
448  // Note: DenseVector::resize() also zeros the vector.
449  dest.resize(this->n());
450 
451  // Short-circuit if the matrix is empty
452  if (this->m() == 0)
453  return;
454 
455  if (this->use_blas_lapack)
456  {
457  this->_matvec_blas(1., 0., dest, arg, /*trans=*/true);
458  }
459  else
460  {
461  const unsigned int n_rows = this->m();
462  const unsigned int n_cols = this->n();
463 
464  // WORKS
465  // for (unsigned int j=0; j<n_cols; j++)
466  // for (unsigned int i=0; i<n_rows; i++)
467  // dest(j) += (*this)(i,j)*arg(i);
468 
469  // ALSO WORKS, (i,j) just swapped
470  for (unsigned int i=0; i<n_cols; i++)
471  for (unsigned int j=0; j<n_rows; j++)
472  dest(i) += (*this)(j,i)*arg(j);
473  }
474 }
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:525
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 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 480 of file dense_matrix_impl.h.

482 {
483  // Make sure the input sizes are compatible
484  libmesh_assert_equal_to (this->m(), arg.size());
485 
486  // Resize and clear dest.
487  // Note: DenseVector::resize() also zeros the vector.
488  dest.resize(this->n());
489 
490  // Short-circuit if the matrix is empty
491  if (this->m() == 0)
492  return;
493 
494  const unsigned int n_rows = this->m();
495  const unsigned int n_cols = this->n();
496 
497  // WORKS
498  // for (unsigned int j=0; j<n_cols; j++)
499  // for (unsigned int i=0; i<n_rows; i++)
500  // dest(j) += (*this)(i,j)*arg(i);
501 
502  // ALSO WORKS, (i,j) just swapped
503  for (unsigned int i=0; i<n_cols; i++)
504  for (unsigned int j=0; j<n_rows; j++)
505  dest(i) += (*this)(j,i)*arg(j);
506 }
unsigned int m() const
unsigned int n() const

◆ zero()

template<typename T >
void libMesh::DenseMatrix< T >::zero ( )
overridevirtual

Set every element in the matrix to 0.

You must redefine what you mean by zeroing the matrix since it depends on how your values are stored.

Implements libMesh::DenseMatrixBase< T >.

Definition at line 808 of file dense_matrix.h.

Referenced by libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction >::operator()().

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

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

◆ _m

template<typename T>
unsigned int libMesh::DenseMatrixBase< T >::_m
protectedinherited

The row dimension.

Definition at line 169 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 174 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 672 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 532 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 525 of file dense_matrix.h.


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