libMesh::DenseMatrix< T > Class Template Reference

A matrix object used for finite element assembly and numerics. More...

#include <dense_matrix.h>

Inheritance diagram for libMesh::DenseMatrix< T >:

Public Member Functions

 DenseMatrix (const unsigned int new_m=0, const unsigned int new_n=0)
 
virtual ~DenseMatrix ()
 
virtual void zero () libmesh_override
 
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 libmesh_override
 
virtual T & el (const unsigned int i, const unsigned int j) libmesh_override
 
virtual void left_multiply (const DenseMatrixBase< T > &M2) libmesh_override
 
template<typename T2 >
void left_multiply (const DenseMatrixBase< T2 > &M2)
 
virtual void right_multiply (const DenseMatrixBase< T > &M2) libmesh_override
 
template<typename T2 >
void right_multiply (const DenseMatrixBase< T2 > &M2)
 
void vector_mult (DenseVector< T > &dest, const DenseVector< T > &arg) const
 
template<typename T2 >
void vector_mult (DenseVector< typename CompareTypes< T, T2 >::supertype > &dest, const DenseVector< T2 > &arg) const
 
void vector_mult_transpose (DenseVector< T > &dest, const DenseVector< T > &arg) const
 
template<typename T2 >
void vector_mult_transpose (DenseVector< typename CompareTypes< T, T2 >::supertype > &dest, const DenseVector< T2 > &arg) const
 
void vector_mult_add (DenseVector< T > &dest, const T factor, const DenseVector< T > &arg) const
 
template<typename T2 , typename T3 >
void vector_mult_add (DenseVector< typename CompareTypes< T, typename CompareTypes< T2, T3 >::supertype >::supertype > &dest, const T2 factor, const DenseVector< T3 > &arg) const
 
void get_principal_submatrix (unsigned int sub_m, unsigned int sub_n, DenseMatrix< T > &dest) const
 
void get_principal_submatrix (unsigned int sub_m, DenseMatrix< T > &dest) const
 
DenseMatrix< T > & operator= (const DenseMatrix< T > &other_matrix)
 
template<typename T2 >
DenseMatrix< T > & operator= (const DenseMatrix< T2 > &other_matrix)
 
void swap (DenseMatrix< T > &other_matrix)
 
void resize (const unsigned int new_m, const unsigned int new_n)
 
void scale (const T factor)
 
void scale_column (const unsigned int col, const T factor)
 
DenseMatrix< T > & operator*= (const T factor)
 
template<typename T2 , typename T3 >
boostcopy::enable_if_c< ScalarTraits< T2 >::value, void >::type add (const T2 factor, const DenseMatrix< T3 > &mat)
 
bool operator== (const DenseMatrix< T > &mat) const
 
bool operator!= (const DenseMatrix< T > &mat) const
 
DenseMatrix< T > & operator+= (const DenseMatrix< T > &mat)
 
DenseMatrix< T > & operator-= (const DenseMatrix< T > &mat)
 
Real min () const
 
Real max () const
 
Real l1_norm () const
 
Real linfty_norm () const
 
void left_multiply_transpose (const DenseMatrix< T > &A)
 
template<typename T2 >
void left_multiply_transpose (const DenseMatrix< T2 > &A)
 
void right_multiply_transpose (const DenseMatrix< T > &A)
 
template<typename T2 >
void right_multiply_transpose (const DenseMatrix< T2 > &A)
 
transpose (const unsigned int i, const unsigned int j) const
 
void get_transpose (DenseMatrix< T > &dest) const
 
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)
 
void lu_solve (const DenseVector< T > &b, DenseVector< T > &x)
 
template<typename T2 >
void cholesky_solve (const DenseVector< T2 > &b, DenseVector< T2 > &x)
 
void svd (DenseVector< Real > &sigma)
 
void svd (DenseVector< Real > &sigma, DenseMatrix< Number > &U, DenseMatrix< Number > &VT)
 
void svd_solve (const DenseVector< T > &rhs, DenseVector< T > &x, Real rcond=std::numeric_limits< Real >::epsilon()) const
 
void evd (DenseVector< T > &lambda_real, DenseVector< T > &lambda_imag)
 
void evd_left (DenseVector< T > &lambda_real, DenseVector< T > &lambda_imag, DenseMatrix< T > &VL)
 
void evd_right (DenseVector< T > &lambda_real, DenseVector< T > &lambda_imag, DenseMatrix< T > &VR)
 
void evd_left_and_right (DenseVector< T > &lambda_real, DenseVector< T > &lambda_imag, DenseMatrix< T > &VL, DenseMatrix< T > &VR)
 
det ()
 
unsigned int m () const
 
unsigned int n () const
 
void print (std::ostream &os=libMesh::out) const
 
void print_scientific (std::ostream &os, unsigned precision=8) const
 
template<typename T2 , typename T3 >
boostcopy::enable_if_c< ScalarTraits< T2 >::value, void >::type add (const T2 factor, const DenseMatrixBase< T3 > &mat)
 

Public Attributes

bool use_blas_lapack
 

Protected Member Functions

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

Static Protected Member Functions

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

Protected Attributes

unsigned int _m
 
unsigned int _n
 

Private Types

enum  DecompositionType { LU =0, CHOLESKY =1, LU_BLAS_LAPACK, NONE }
 
enum  _BLAS_Multiply_Flag { LEFT_MULTIPLY = 0, RIGHT_MULTIPLY, LEFT_MULTIPLY_TRANSPOSE, RIGHT_MULTIPLY_TRANSPOSE }
 
typedef PetscBLASInt pivot_index_t
 
typedef int pivot_index_t
 

Private Member Functions

void _lu_decompose ()
 
void _lu_back_substitute (const DenseVector< T > &b, DenseVector< T > &x) const
 
void _cholesky_decompose ()
 
template<typename T2 >
void _cholesky_back_substitute (const DenseVector< T2 > &b, DenseVector< T2 > &x) const
 
void _multiply_blas (const DenseMatrixBase< T > &other, _BLAS_Multiply_Flag flag)
 
void _lu_decompose_lapack ()
 
void _svd_lapack (DenseVector< Real > &sigma)
 
void _svd_lapack (DenseVector< Real > &sigma, DenseMatrix< Number > &U, DenseMatrix< Number > &VT)
 
void _svd_solve_lapack (const DenseVector< T > &rhs, DenseVector< T > &x, Real rcond) const
 
void _svd_helper (char JOBU, char JOBVT, std::vector< Real > &sigma_val, std::vector< Number > &U_val, std::vector< Number > &VT_val)
 
void _evd_lapack (DenseVector< T > &lambda_real, DenseVector< T > &lambda_imag, DenseMatrix< T > *VL=libmesh_nullptr, DenseMatrix< T > *VR=libmesh_nullptr)
 
void _lu_back_substitute_lapack (const DenseVector< T > &b, DenseVector< T > &x)
 
void _matvec_blas (T alpha, T beta, DenseVector< T > &dest, const DenseVector< T > &arg, bool trans=false) const
 

Private Attributes

std::vector< T > _val
 
DecompositionType _decomposition_type
 
std::vector< pivot_index_t_pivots
 

Detailed Description

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

A matrix object used for finite element assembly and numerics.

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

Definition at line 54 of file dense_matrix.h.

Member Typedef Documentation

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

template<typename T>
typedef int libMesh::DenseMatrix< T >::pivot_index_t
private

Definition at line 654 of file dense_matrix.h.

Member Enumeration Documentation

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

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

Constructor & Destructor Documentation

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

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

729  :
730  DenseMatrixBase<T>(new_m,new_n),
731 #if defined(LIBMESH_HAVE_PETSC) && defined(LIBMESH_USE_REAL_NUMBERS) && defined(LIBMESH_DEFAULT_DOUBLE_PRECISION)
732  use_blas_lapack(true),
733 #else
734  use_blas_lapack(false),
735 #endif
736  _val(),
738 {
739  this->resize(new_m,new_n);
740 }
DecompositionType _decomposition_type
Definition: dense_matrix.h:566
void resize(const unsigned int new_m, const unsigned int new_n)
Definition: dense_matrix.h:776
std::vector< T > _val
Definition: dense_matrix.h:516
template<typename T>
virtual libMesh::DenseMatrix< T >::~DenseMatrix ( )
inlinevirtual

Destructor. Empty.

Definition at line 67 of file dense_matrix.h.

67 {}

Member Function Documentation

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

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

994 {
995  // Shorthand notation for number of rows and columns.
996  const unsigned int
997  n_rows = this->m(),
998  n_cols = this->n();
999 
1000  // Just to be really sure...
1001  libmesh_assert_equal_to (n_rows, n_cols);
1002 
1003  // A convenient reference to *this
1004  const DenseMatrix<T> & A = *this;
1005 
1006  // Now compute the solution to Ax =b using the factorization.
1007  x.resize(n_rows);
1008 
1009  // Solve for Ly=b
1010  for (unsigned int i=0; i<n_cols; ++i)
1011  {
1012  T2 temp = b(i);
1013 
1014  for (unsigned int k=0; k<i; ++k)
1015  temp -= A(i,k)*x(k);
1016 
1017  x(i) = temp / A(i,i);
1018  }
1019 
1020  // Solve for L^T x = y
1021  for (unsigned int i=0; i<n_cols; ++i)
1022  {
1023  const unsigned int ib = (n_cols-1)-i;
1024 
1025  for (unsigned int k=(ib+1); k<n_cols; ++k)
1026  x(ib) -= A(k,ib) * x(k);
1027 
1028  x(ib) /= A(ib,ib);
1029  }
1030 }
unsigned int n() const
unsigned int m() const
static PetscErrorCode Mat * A
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 946 of file dense_matrix_impl.h.

References A.

947 {
948  // If we called this function, there better not be any
949  // previous decomposition of the matrix.
950  libmesh_assert_equal_to (this->_decomposition_type, NONE);
951 
952  // Shorthand notation for number of rows and columns.
953  const unsigned int
954  n_rows = this->m(),
955  n_cols = this->n();
956 
957  // Just to be really sure...
958  libmesh_assert_equal_to (n_rows, n_cols);
959 
960  // A convenient reference to *this
961  DenseMatrix<T> & A = *this;
962 
963  for (unsigned int i=0; i<n_rows; ++i)
964  {
965  for (unsigned int j=i; j<n_cols; ++j)
966  {
967  for (unsigned int k=0; k<i; ++k)
968  A(i,j) -= A(i,k) * A(j,k);
969 
970  if (i == j)
971  {
972 #ifndef LIBMESH_USE_COMPLEX_NUMBERS
973  if (A(i,j) <= 0.0)
974  libmesh_error_msg("Error! Can only use Cholesky decomposition with symmetric positive definite matrices.");
975 #endif
976 
977  A(i,i) = std::sqrt(A(i,j));
978  }
979  else
980  A(j,i) = A(i,j) / A(i,i);
981  }
982  }
983 
984  // Set the flag for CHOLESKY decomposition
986 }
unsigned int n() const
DecompositionType _decomposition_type
Definition: dense_matrix.h:566
unsigned int m() const
static PetscErrorCode Mat * A
template<typename T>
template void libMesh::DenseMatrix< T >::_evd_lapack ( DenseVector< T > &  lambda_real,
DenseVector< T > &  lambda_imag,
DenseMatrix< T > *  VL = libmesh_nullptr,
DenseMatrix< T > *  VR = libmesh_nullptr 
)
private

Computes the eigenvalues of the matrix using the Lapack routine "DGEEV". If VR and/or VL are non-NULL, 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 704 of file dense_matrix_blas_lapack.C.

References libMesh::DenseVector< T >::get_values(), libMesh::DenseMatrix< T >::get_values(), libmesh_nullptr, libMesh::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::resize(), and swap().

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

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

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

653 {
654  const unsigned int
655  n_cols = this->n();
656 
657  libmesh_assert_equal_to (this->m(), n_cols);
658  libmesh_assert_equal_to (this->m(), b.size());
659 
660  x.resize (n_cols);
661 
662  // A convenient reference to *this
663  const DenseMatrix<T> & A = *this;
664 
665  // Temporary vector storage. We use this instead of
666  // modifying the RHS.
667  DenseVector<T> z = b;
668 
669  // Lower-triangular "top to bottom" solve step, taking into account pivots
670  for (unsigned int i=0; i<n_cols; ++i)
671  {
672  // Swap
673  if (_pivots[i] != static_cast<pivot_index_t>(i))
674  std::swap( z(i), z(_pivots[i]) );
675 
676  x(i) = z(i);
677 
678  for (unsigned int j=0; j<i; ++j)
679  x(i) -= A(i,j)*x(j);
680 
681  x(i) /= A(i,i);
682  }
683 
684  // Upper-triangular "bottom to top" solve step
685  const unsigned int last_row = n_cols-1;
686 
687  for (int i=last_row; i>=0; --i)
688  {
689  for (int j=i+1; j<static_cast<int>(n_cols); ++j)
690  x(i) -= A(i,j)*x(j);
691  }
692 }
unsigned int n() const
unsigned int m() const
static PetscErrorCode Mat * A
void swap(Iterator &lhs, Iterator &rhs)
std::vector< pivot_index_t > _pivots
Definition: dense_matrix.h:656
template<typename T>
template 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 913 of file dense_matrix_blas_lapack.C.

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

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

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

References A, std::abs(), libMesh::Real, swap(), and libMesh::zero.

703 {
704  // If this function was called, there better not be any
705  // previous decomposition of the matrix.
706  libmesh_assert_equal_to (this->_decomposition_type, NONE);
707 
708  // Get the matrix size and make sure it is square
709  const unsigned int
710  n_rows = this->m();
711 
712  // A convenient reference to *this
713  DenseMatrix<T> & A = *this;
714 
715  _pivots.resize(n_rows);
716 
717  for (unsigned int i=0; i<n_rows; ++i)
718  {
719  // Find the pivot row by searching down the i'th column
720  _pivots[i] = i;
721 
722  // std::abs(complex) must return a Real!
723  Real the_max = std::abs( A(i,i) );
724  for (unsigned int j=i+1; j<n_rows; ++j)
725  {
726  Real candidate_max = std::abs( A(j,i) );
727  if (the_max < candidate_max)
728  {
729  the_max = candidate_max;
730  _pivots[i] = j;
731  }
732  }
733 
734  // libMesh::out << "the_max=" << the_max << " found at row " << _pivots[i] << std::endl;
735 
736  // If the max was found in a different row, interchange rows.
737  // Here we interchange the *entire* row, in Gaussian elimination
738  // you would only interchange the subrows A(i,j) and A(p(i),j), for j>i
739  if (_pivots[i] != static_cast<pivot_index_t>(i))
740  {
741  for (unsigned int j=0; j<n_rows; ++j)
742  std::swap( A(i,j), A(_pivots[i], j) );
743  }
744 
745 
746  // If the max abs entry found is zero, the matrix is singular
747  if (A(i,i) == libMesh::zero)
748  libmesh_error_msg("Matrix A is singular!");
749 
750  // Scale upper triangle entries of row i by the diagonal entry
751  // Note: don't scale the diagonal entry itself!
752  const T diag_inv = 1. / A(i,i);
753  for (unsigned int j=i+1; j<n_rows; ++j)
754  A(i,j) *= diag_inv;
755 
756  // Update the remaining sub-matrix A[i+1:m][i+1:m]
757  // by subtracting off (the diagonal-scaled)
758  // upper-triangular part of row i, scaled by the
759  // i'th column entry of each row. In terms of
760  // row operations, this is:
761  // for each r > i
762  // SubRow(r) = SubRow(r) - A(r,i)*SubRow(i)
763  //
764  // If we were scaling the i'th column as well, like
765  // in Gaussian elimination, this would 'zero' the
766  // entry in the i'th column.
767  for (unsigned int row=i+1; row<n_rows; ++row)
768  for (unsigned int col=i+1; col<n_rows; ++col)
769  A(row,col) -= A(row,i) * A(i,col);
770 
771  } // end i loop
772 
773  // Set the flag for LU decomposition
774  this->_decomposition_type = LU;
775 }
double abs(double a)
DecompositionType _decomposition_type
Definition: dense_matrix.h:566
const Number zero
Definition: libmesh.h:178
unsigned int m() const
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static PetscErrorCode Mat * A
void swap(Iterator &lhs, Iterator &rhs)
std::vector< pivot_index_t > _pivots
Definition: dense_matrix.h:656
template<typename T >
template 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 205 of file dense_matrix_blas_lapack.C.

References std::min().

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

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

References libMesh::DenseMatrix< T >::_evd_lapack(), libMesh::DenseMatrix< T >::_lu_back_substitute_lapack(), libMesh::DenseMatrix< T >::_lu_decompose_lapack(), libMesh::DenseMatrix< T >::_multiply_blas(), libMesh::DenseMatrix< T >::_svd_helper(), libMesh::DenseMatrix< T >::_svd_lapack(), libMesh::DenseMatrix< T >::_svd_solve_lapack(), libMesh::DenseVector< T >::get_values(), libMesh::DenseMatrix< T >::get_values(), libMesh::Real, and libMesh::DenseVector< T >::size().

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

References libMesh::DenseMatrix< T >::_val, A, libMesh::DenseMatrixBase< T >::m(), libMesh::DenseMatrixBase< T >::n(), and swap().

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

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

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

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

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

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

274 {
275  // The calling sequence for dgetrf is:
276  // DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, INFO )
277 
278  // JOBU (input)
279  // Specifies options for computing all or part of the matrix U:
280  // = 'A': all M columns of U are returned in array U:
281  // = 'S': the first min(m,n) columns of U (the left singular
282  // vectors) are returned in the array U;
283  // = 'O': the first min(m,n) columns of U (the left singular
284  // vectors) are overwritten on the array A;
285  // = 'N': no columns of U (no left singular vectors) are
286  // computed.
287  char JOBU = 'N';
288 
289  // JOBVT (input)
290  // Specifies options for computing all or part of the matrix
291  // V**T:
292  // = 'A': all N rows of V**T are returned in the array VT;
293  // = 'S': the first min(m,n) rows of V**T (the right singular
294  // vectors) are returned in the array VT;
295  // = 'O': the first min(m,n) rows of V**T (the right singular
296  // vectors) are overwritten on the array A;
297  // = 'N': no rows of V**T (no right singular vectors) are
298  // computed.
299  char JOBVT = 'N';
300 
301  std::vector<Real> sigma_val;
302  std::vector<Number> U_val;
303  std::vector<Number> VT_val;
304 
305  _svd_helper(JOBU, JOBVT, sigma_val, U_val, VT_val);
306 
307  // Copy the singular values into sigma, ignore U_val and VT_val
308  sigma.resize(cast_int<unsigned int>(sigma_val.size()));
309  for (unsigned int i=0; i<sigma.size(); i++)
310  sigma(i) = sigma_val[i];
311 
312 }
void _svd_helper(char JOBU, char JOBVT, std::vector< Real > &sigma_val, std::vector< Number > &U_val, std::vector< Number > &VT_val)
template<typename T >
template 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 315 of file dense_matrix_blas_lapack.C.

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

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

References libMesh::DenseVector< T >::get_values(), libMesh::DenseMatrix< T >::get_values(), libMesh::DenseMatrixBase< T >::m(), std::max(), std::min(), libMesh::DenseMatrixBase< T >::n(), libMesh::Real, libMesh::DenseVector< T >::resize(), and libMesh::DenseVector< T >::size().

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

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

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

185 {
186  libmesh_assert_equal_to (this->m(), mat.m());
187  libmesh_assert_equal_to (this->n(), mat.n());
188 
189  for (unsigned int j=0; j<this->n(); j++)
190  for (unsigned int i=0; i<this->m(); i++)
191  this->el(i,j) += factor*mat.el(i,j);
192 }
unsigned int n() const
unsigned int m() const
virtual T el(const unsigned int i, const unsigned int j) const =0
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 883 of file dense_matrix.h.

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

Referenced by libMesh::FEMSystem::assembly(), and libMesh::DenseMatrix< Number >::el().

885 {
886  libmesh_assert_equal_to (this->m(), mat.m());
887  libmesh_assert_equal_to (this->n(), mat.n());
888 
889  for (unsigned int i=0; i<this->m(); i++)
890  for (unsigned int j=0; j<this->n(); j++)
891  (*this)(i,j) += factor * mat(i,j);
892 }
unsigned int n() const
unsigned int m() const
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 914 of file dense_matrix_impl.h.

Referenced by libMesh::FEGenericBase< OutputType >::coarsened_dof_values(), libMesh::FEGenericBase< OutputType >::compute_periodic_constraints(), libMesh::FEGenericBase< OutputType >::compute_proj_constraints(), libMesh::DenseMatrix< Number >::condense(), libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction >::operator()(), libMesh::BoundaryProjectSolution::operator()(), and libMesh::HPCoarsenTest::select_refinement().

916 {
917  // Check for a previous decomposition
918  switch(this->_decomposition_type)
919  {
920  case NONE:
921  {
922  this->_cholesky_decompose ();
923  break;
924  }
925 
926  case CHOLESKY:
927  {
928  // Already factored, just need to call back_substitute.
929  break;
930  }
931 
932  default:
933  libmesh_error_msg("Error! This matrix already has a different decomposition...");
934  }
935 
936  // Perform back substitution
937  this->_cholesky_back_substitute (b, x);
938 }
DecompositionType _decomposition_type
Definition: dense_matrix.h:566
void _cholesky_back_substitute(const DenseVector< T2 > &b, DenseVector< T2 > &x) const
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 58 of file dense_matrix_base.C.

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

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

62 {
63  libmesh_assert_equal_to (this->_m, rhs.size());
64  libmesh_assert_equal_to (iv, jv);
65 
66 
67  // move the known value into the RHS
68  // and zero the column
69  for (unsigned int i=0; i<this->m(); i++)
70  {
71  rhs.el(i) -= this->el(i,jv)*val;
72  this->el(i,jv) = 0.;
73  }
74 
75  // zero the row
76  for (unsigned int j=0; j<this->n(); j++)
77  this->el(iv,j) = 0.;
78 
79  this->el(iv,jv) = 1.;
80  rhs.el(iv) = val;
81 
82 }
unsigned int n() const
unsigned int m() const
virtual T el(const unsigned int i, const unsigned int j) const =0
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 338 of file dense_matrix.h.

342  { DenseMatrixBase<T>::condense (i, j, val, rhs); }
void condense(const unsigned int i, const unsigned int j, const T val, DenseVectorBase< T > &rhs)
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 853 of file dense_matrix_impl.h.

References libMesh::Real.

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

854 {
855  switch(this->_decomposition_type)
856  {
857  case NONE:
858  {
859  // First LU decompose the matrix.
860  // Note that the lu_decompose routine will check to see if the
861  // matrix is square so we don't worry about it.
862  if (this->use_blas_lapack)
863  this->_lu_decompose_lapack();
864  else
865  this->_lu_decompose ();
866  }
867  case LU:
868  case LU_BLAS_LAPACK:
869  {
870  // Already decomposed, don't do anything
871  break;
872  }
873  default:
874  libmesh_error_msg("Error! Can't compute the determinant under the current decomposition.");
875  }
876 
877  // A variable to keep track of the running product of diagonal terms.
878  T determinant = 1.;
879 
880  // Loop over diagonal terms, computing the product. In practice,
881  // be careful because this value could easily become too large to
882  // fit in a double or float. To be safe, one should keep track of
883  // the power (of 10) of the determinant in a separate variable
884  // and maintain an order 1 value for the determinant itself.
885  unsigned int n_interchanges = 0;
886  for (unsigned int i=0; i<this->m(); i++)
887  {
888  if (this->_decomposition_type==LU)
889  if (_pivots[i] != static_cast<pivot_index_t>(i))
890  n_interchanges++;
891 
892  // Lapack pivots are 1-based!
894  if (_pivots[i] != static_cast<pivot_index_t>(i+1))
895  n_interchanges++;
896 
897  determinant *= (*this)(i,i);
898  }
899 
900  // Compute sign of determinant, depends on number of row interchanges!
901  // The sign should be (-1)^{n}, where n is the number of interchanges.
902  Real sign = n_interchanges % 2 == 0 ? 1. : -1.;
903 
904  return sign*determinant;
905 }
DecompositionType _decomposition_type
Definition: dense_matrix.h:566
unsigned int m() const
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
std::vector< pivot_index_t > _pivots
Definition: dense_matrix.h:656
template<typename T>
virtual T libMesh::DenseMatrix< T >::el ( const unsigned int  i,
const unsigned int  j 
) const
inlinevirtual
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 84 of file dense_matrix.h.

86  { return (*this)(i,j); }
template<typename T>
virtual T& libMesh::DenseMatrix< T >::el ( const unsigned int  i,
const unsigned int  j 
)
inlinevirtual
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 88 of file dense_matrix.h.

90  { return (*this)(i,j); }
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 809 of file dense_matrix_impl.h.

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

811 {
812  // We use the LAPACK eigenvalue problem implementation
813  _evd_lapack(lambda_real, lambda_imag);
814 }
void _evd_lapack(DenseVector< T > &lambda_real, DenseVector< T > &lambda_imag, DenseMatrix< T > *VL=libmesh_nullptr, DenseMatrix< T > *VR=libmesh_nullptr)
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 819 of file dense_matrix_impl.h.

References libmesh_nullptr.

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

822 {
823  // We use the LAPACK eigenvalue problem implementation
824  _evd_lapack(lambda_real, lambda_imag, &VL, libmesh_nullptr);
825 }
const class libmesh_nullptr_t libmesh_nullptr
void _evd_lapack(DenseVector< T > &lambda_real, DenseVector< T > &lambda_imag, DenseMatrix< T > *VL=libmesh_nullptr, DenseMatrix< T > *VR=libmesh_nullptr)
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 841 of file dense_matrix_impl.h.

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

845 {
846  // We use the LAPACK eigenvalue problem implementation
847  _evd_lapack(lambda_real, lambda_imag, &VL, &VR);
848 }
void _evd_lapack(DenseVector< T > &lambda_real, DenseVector< T > &lambda_imag, DenseMatrix< T > *VL=libmesh_nullptr, DenseMatrix< T > *VR=libmesh_nullptr)
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 830 of file dense_matrix_impl.h.

References libmesh_nullptr.

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

833 {
834  // We use the LAPACK eigenvalue problem implementation
835  _evd_lapack(lambda_real, lambda_imag, libmesh_nullptr, &VR);
836 }
const class libmesh_nullptr_t libmesh_nullptr
void _evd_lapack(DenseVector< T > &lambda_real, DenseVector< T > &lambda_imag, DenseMatrix< T > *VL=libmesh_nullptr, DenseMatrix< T > *VR=libmesh_nullptr)
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.

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

Referenced by libMesh::DenseMatrix< Number >::el().

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 n() const
unsigned int m() const
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 571 of file dense_matrix_impl.h.

572 {
573  get_principal_submatrix(sub_m, sub_m, dest);
574 }
void get_principal_submatrix(unsigned int sub_m, unsigned int sub_n, DenseMatrix< T > &dest) const
template<typename T>
void libMesh::DenseMatrix< T >::get_transpose ( DenseMatrix< T > &  dest) const

Put the tranposed matrix into dest.

Definition at line 579 of file dense_matrix_impl.h.

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

Referenced by libMesh::DenseMatrix< Number >::el().

580 {
581  dest.resize(this->n(), this->m());
582 
583  for (unsigned int i=0; i<dest.m(); i++)
584  for (unsigned int j=0; j<dest.n(); j++)
585  dest(i,j) = (*this)(j,i);
586 }
unsigned int n() const
unsigned int m() const
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 325 of file dense_matrix.h.

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

325 { return _val; }
std::vector< T > _val
Definition: dense_matrix.h:516
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 330 of file dense_matrix.h.

330 { return _val; }
std::vector< T > _val
Definition: dense_matrix.h:516
template<typename T >
Real libMesh::DenseMatrix< T >::l1_norm ( ) const
inline
Returns
The l1-norm of the matrix, that is, the max column sum:

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

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

Definition at line 990 of file dense_matrix.h.

References libMesh::DenseMatrixBase< T >::_m, libMesh::DenseMatrixBase< T >::_n, std::abs(), and libMesh::Real.

Referenced by libMesh::FEMSystem::assembly(), and libMesh::DenseMatrix< Number >::el().

991 {
992  libmesh_assert (this->_m);
993  libmesh_assert (this->_n);
994 
995  Real columnsum = 0.;
996  for (unsigned int i=0; i!=this->_m; i++)
997  {
998  columnsum += std::abs((*this)(i,0));
999  }
1000  Real my_max = columnsum;
1001  for (unsigned int j=1; j!=this->_n; j++)
1002  {
1003  columnsum = 0.;
1004  for (unsigned int i=0; i!=this->_m; i++)
1005  {
1006  columnsum += std::abs((*this)(i,j));
1007  }
1008  my_max = (my_max > columnsum? my_max : columnsum);
1009  }
1010  return my_max;
1011 }
double abs(double a)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
template<typename T>
void libMesh::DenseMatrix< T >::left_multiply ( const DenseMatrixBase< T > &  M2)
virtual

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

Implements libMesh::DenseMatrixBase< T >.

Definition at line 40 of file dense_matrix_impl.h.

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

Referenced by libMesh::DenseMatrix< Number >::el().

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 }
void _multiply_blas(const DenseMatrixBase< T > &other, _BLAS_Multiply_Flag flag)
void resize(const unsigned int new_m, const unsigned int new_n)
Definition: dense_matrix.h:776
static void multiply(DenseMatrixBase< T > &M1, const DenseMatrixBase< T > &M2, const DenseMatrixBase< T > &M3)
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.

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

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 }
void resize(const unsigned int new_m, const unsigned int new_n)
Definition: dense_matrix.h:776
static void multiply(DenseMatrixBase< T > &M1, const DenseMatrixBase< T > &M2, const DenseMatrixBase< T > &M3)
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.

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

Referenced by libMesh::DofMap::constrain_element_matrix(), libMesh::DofMap::constrain_element_matrix_and_vector(), libMesh::DenseMatrix< Number >::el(), and libMesh::DofMap::heterogenously_constrain_element_matrix_and_vector().

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 }
unsigned int n() const
void _multiply_blas(const DenseMatrixBase< T > &other, _BLAS_Multiply_Flag flag)
unsigned int m() const
static PetscErrorCode Mat * A
void resize(const unsigned int new_m, const unsigned int new_n)
Definition: dense_matrix.h:776
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.

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

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 n() const
unsigned int m() const
static PetscErrorCode Mat * A
void resize(const unsigned int new_m, const unsigned int new_n)
Definition: dense_matrix.h:776
template<typename T >
Real libMesh::DenseMatrix< T >::linfty_norm ( ) const
inline
Returns
The linfty-norm of the matrix, that is, the max row sum:

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

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

Definition at line 1017 of file dense_matrix.h.

References libMesh::DenseMatrixBase< T >::_m, libMesh::DenseMatrixBase< T >::_n, std::abs(), and libMesh::Real.

Referenced by libMesh::DenseMatrix< Number >::el().

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

Referenced by libMesh::DenseMatrix< Number >::condense(), libMesh::WeightedPatchRecoveryErrorEstimator::EstimateError::operator()(), and libMesh::PatchRecoveryErrorEstimator::EstimateError::operator()().

594 {
595  // Check to be sure that the matrix is square before attempting
596  // an LU-solve. In general, one can compute the LU factorization of
597  // a non-square matrix, but:
598  //
599  // Overdetermined systems (m>n) have a solution only if enough of
600  // the equations are linearly-dependent.
601  //
602  // Underdetermined systems (m<n) typically have infinitely many
603  // solutions.
604  //
605  // We don't want to deal with either of these ambiguous cases here...
606  libmesh_assert_equal_to (this->m(), this->n());
607 
608  switch(this->_decomposition_type)
609  {
610  case NONE:
611  {
612  if (this->use_blas_lapack)
613  this->_lu_decompose_lapack();
614  else
615  this->_lu_decompose ();
616  break;
617  }
618 
619  case LU_BLAS_LAPACK:
620  {
621  // Already factored, just need to call back_substitute.
622  if (this->use_blas_lapack)
623  break;
624  }
625  libmesh_fallthrough();
626 
627  case LU:
628  {
629  // Already factored, just need to call back_substitute.
630  if (!(this->use_blas_lapack))
631  break;
632  }
633  libmesh_fallthrough();
634 
635  default:
636  libmesh_error_msg("Error! This matrix already has a different decomposition...");
637  }
638 
639  if (this->use_blas_lapack)
640  this->_lu_back_substitute_lapack (b, x);
641  else
642  this->_lu_back_substitute (b, x);
643 }
void _lu_back_substitute(const DenseVector< T > &b, DenseVector< T > &x) const
unsigned int n() const
DecompositionType _decomposition_type
Definition: dense_matrix.h:566
void _lu_back_substitute_lapack(const DenseVector< T > &b, DenseVector< T > &x)
unsigned int m() const
template<typename T>
unsigned int libMesh::DenseMatrixBase< T >::m ( ) const
inlineinherited
Returns
The row-dimension of the matrix.

Definition at line 97 of file dense_matrix_base.h.

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

Referenced by libMesh::DenseMatrix< T >::_multiply_blas(), libMesh::DenseMatrix< T >::_svd_solve_lapack(), libMesh::DenseMatrixBase< T >::add(), libMesh::DenseMatrix< T >::add(), libMesh::PetscMatrix< T >::add_block_matrix(), libMesh::SparseMatrix< T >::add_block_matrix(), libMesh::EigenSparseMatrix< T >::add_matrix(), libMesh::LaspackMatrix< T >::add_matrix(), libMesh::EpetraMatrix< T >::add_matrix(), libMesh::PetscMatrix< T >::add_matrix(), libMesh::DofMap::build_constraint_matrix(), libMesh::DofMap::build_constraint_matrix_and_vector(), libMesh::DofMap::constrain_element_dyad_matrix(), libMesh::DofMap::constrain_element_matrix(), libMesh::DofMap::constrain_element_matrix_and_vector(), libMesh::DofMap::constrain_element_vector(), libMesh::DofMap::extract_local_vector(), libMesh::DenseMatrix< T >::get_transpose(), libMesh::DofMap::heterogenously_constrain_element_matrix_and_vector(), libMesh::DofMap::heterogenously_constrain_element_vector(), libMesh::DenseMatrix< T >::left_multiply(), libMesh::DenseMatrix< T >::left_multiply_transpose(), libMesh::DofMap::max_constraint_error(), libMesh::DenseMatrixBase< T >::multiply(), libMesh::DenseSubMatrix< T >::operator()(), libMesh::WeightedPatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::PatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::DenseMatrix< T >::operator=(), libMesh::DenseSubMatrix< T >::reposition(), libMesh::DenseMatrix< T >::right_multiply(), libMesh::DenseMatrix< T >::right_multiply_transpose(), libMesh::DenseMatrix< T >::scale_column(), and libMesh::DenseSubMatrix< T >::zero().

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

Definition at line 969 of file dense_matrix.h.

References libMesh::DenseMatrixBase< T >::_m, libMesh::DenseMatrixBase< T >::_n, libMesh::libmesh_real(), and libMesh::Real.

Referenced by libMesh::DenseMatrix< Number >::el().

970 {
971  libmesh_assert (this->_m);
972  libmesh_assert (this->_n);
973  Real my_max = libmesh_real((*this)(0,0));
974 
975  for (unsigned int i=0; i!=this->_m; i++)
976  {
977  for (unsigned int j=0; j!=this->_n; j++)
978  {
979  Real current = libmesh_real((*this)(i,j));
980  my_max = (my_max > current? my_max : current);
981  }
982  }
983  return my_max;
984 }
T libmesh_real(T a)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
template<typename T >
Real libMesh::DenseMatrix< T >::min ( ) const
inline
Returns
The minimum entry in the matrix, or the minimum real part in the case of complex numbers.

Definition at line 948 of file dense_matrix.h.

References libMesh::DenseMatrixBase< T >::_m, libMesh::DenseMatrixBase< T >::_n, libMesh::libmesh_real(), and libMesh::Real.

Referenced by libMesh::DenseMatrix< Number >::el().

949 {
950  libmesh_assert (this->_m);
951  libmesh_assert (this->_n);
952  Real my_min = libmesh_real((*this)(0,0));
953 
954  for (unsigned int i=0; i!=this->_m; i++)
955  {
956  for (unsigned int j=0; j!=this->_n; j++)
957  {
958  Real current = libmesh_real((*this)(i,j));
959  my_min = (my_min < current? my_min : current);
960  }
961  }
962  return my_min;
963 }
T libmesh_real(T a)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
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 31 of file dense_matrix_base.C.

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

34 {
35  // Assertions to make sure we have been
36  // passed matrices of the correct dimension.
37  libmesh_assert_equal_to (M1.m(), M2.m());
38  libmesh_assert_equal_to (M1.n(), M3.n());
39  libmesh_assert_equal_to (M2.n(), M3.m());
40 
41  const unsigned int m_s = M2.m();
42  const unsigned int p_s = M2.n();
43  const unsigned int n_s = M1.n();
44 
45  // Do it this way because there is a
46  // decent chance (at least for constraint matrices)
47  // that M3(k,j) = 0. when right-multiplying.
48  for (unsigned int k=0; k<p_s; k++)
49  for (unsigned int j=0; j<n_s; j++)
50  if (M3.el(k,j) != 0.)
51  for (unsigned int i=0; i<m_s; i++)
52  M1.el(i,j) += M2.el(i,k) * M3.el(k,j);
53 }
template<typename T>
unsigned int libMesh::DenseMatrixBase< T >::n ( ) const
inlineinherited
Returns
The column-dimension of the matrix.

Definition at line 102 of file dense_matrix_base.h.

References libMesh::DenseMatrixBase< T >::_n, libMesh::out, and libMesh::DenseMatrixBase< T >::print().

Referenced by libMesh::DenseMatrix< T >::_multiply_blas(), libMesh::DenseMatrix< T >::_svd_solve_lapack(), libMesh::DenseMatrixBase< T >::add(), libMesh::DenseMatrix< T >::add(), libMesh::PetscMatrix< T >::add_block_matrix(), libMesh::SparseMatrix< T >::add_block_matrix(), libMesh::EigenSparseMatrix< T >::add_matrix(), libMesh::LaspackMatrix< T >::add_matrix(), libMesh::EpetraMatrix< T >::add_matrix(), libMesh::PetscMatrix< T >::add_matrix(), libMesh::DofMap::build_constraint_matrix(), libMesh::DofMap::build_constraint_matrix_and_vector(), libMesh::DofMap::constrain_element_dyad_matrix(), libMesh::DofMap::constrain_element_matrix(), libMesh::DofMap::constrain_element_matrix_and_vector(), libMesh::DofMap::constrain_element_vector(), libMesh::DofMap::extract_local_vector(), libMesh::DenseMatrix< T >::get_transpose(), libMesh::DofMap::heterogenously_constrain_element_matrix_and_vector(), libMesh::DofMap::heterogenously_constrain_element_vector(), libMesh::DenseMatrix< T >::left_multiply(), libMesh::DenseMatrix< T >::left_multiply_transpose(), libMesh::DofMap::max_constraint_error(), libMesh::DenseMatrixBase< T >::multiply(), libMesh::DenseSubMatrix< T >::operator()(), libMesh::WeightedPatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::PatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::DenseMatrix< T >::operator=(), libMesh::DenseSubMatrix< T >::reposition(), libMesh::DenseMatrix< T >::right_multiply(), libMesh::DenseMatrix< T >::right_multiply_transpose(), and libMesh::DenseSubMatrix< T >::zero().

102 { return _n; }
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 911 of file dense_matrix.h.

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

Referenced by libMesh::DenseMatrix< Number >::el().

912 {
913  for (std::size_t i=0; i<_val.size(); i++)
914  if (_val[i] != mat._val[i])
915  return true;
916 
917  return false;
918 }
std::vector< T > _val
Definition: dense_matrix.h:516
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 818 of file dense_matrix.h.

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

Referenced by libMesh::DenseMatrix< Number >::~DenseMatrix().

820 {
821  libmesh_assert_less (i*j, _val.size());
822  libmesh_assert_less (i, this->_m);
823  libmesh_assert_less (j, this->_n);
824 
825 
826  // return _val[(i) + (this->_m)*(j)]; // col-major
827  return _val[(i)*(this->_n) + (j)]; // row-major
828 }
std::vector< T > _val
Definition: dense_matrix.h:516
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 834 of file dense_matrix.h.

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

836 {
837  libmesh_assert_less (i*j, _val.size());
838  libmesh_assert_less (i, this->_m);
839  libmesh_assert_less (j, this->_n);
840 
841  //return _val[(i) + (this->_m)*(j)]; // col-major
842  return _val[(i)*(this->_n) + (j)]; // row-major
843 }
std::vector< T > _val
Definition: dense_matrix.h:516
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 870 of file dense_matrix.h.

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

Referenced by libMesh::DenseMatrix< Number >::el().

871 {
872  this->scale(factor);
873  return *this;
874 }
void scale(const T factor)
Definition: dense_matrix.h:851
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 924 of file dense_matrix.h.

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

Referenced by libMesh::DenseMatrix< Number >::el().

925 {
926  for (std::size_t i=0; i<_val.size(); i++)
927  _val[i] += mat._val[i];
928 
929  return *this;
930 }
std::vector< T > _val
Definition: dense_matrix.h:516
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 936 of file dense_matrix.h.

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

Referenced by libMesh::DenseMatrix< Number >::el().

937 {
938  for (std::size_t i=0; i<_val.size(); i++)
939  _val[i] -= mat._val[i];
940 
941  return *this;
942 }
std::vector< T > _val
Definition: dense_matrix.h:516
template<typename T>
DenseMatrix< T > & libMesh::DenseMatrix< T >::operator= ( const DenseMatrix< T > &  other_matrix)
inline

Assignment operator.

Returns
A reference to *this.

Definition at line 803 of file dense_matrix.h.

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

Referenced by libMesh::DenseMatrix< Number >::el().

804 {
805  this->_m = other_matrix._m;
806  this->_n = other_matrix._n;
807 
808  _val = other_matrix._val;
809  _decomposition_type = other_matrix._decomposition_type;
810 
811  return *this;
812 }
DecompositionType _decomposition_type
Definition: dense_matrix.h:566
std::vector< T > _val
Definition: dense_matrix.h:516
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 761 of file dense_matrix.h.

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

762 {
763  unsigned int mat_m = mat.m(), mat_n = mat.n();
764  this->resize(mat_m, mat_n);
765  for (unsigned int i=0; i<mat_m; i++)
766  for (unsigned int j=0; j<mat_n; j++)
767  (*this)(i,j) = mat(i,j);
768 
769  return *this;
770 }
void resize(const unsigned int new_m, const unsigned int new_n)
Definition: dense_matrix.h:776
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 898 of file dense_matrix.h.

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

Referenced by libMesh::DenseMatrix< Number >::el().

899 {
900  for (std::size_t i=0; i<_val.size(); i++)
901  if (_val[i] != mat._val[i])
902  return false;
903 
904  return true;
905 }
std::vector< T > _val
Definition: dense_matrix.h:516
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 128 of file dense_matrix_base.C.

Referenced by libMesh::DenseMatrixBase< T >::n().

129 {
130  for (unsigned int i=0; i<this->m(); i++)
131  {
132  for (unsigned int j=0; j<this->n(); j++)
133  os << std::setw(8)
134  << this->el(i,j) << " ";
135 
136  os << std::endl;
137  }
138 
139  return;
140 }
unsigned int n() const
unsigned int m() const
virtual T el(const unsigned int i, const unsigned int j) const =0
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 86 of file dense_matrix_base.C.

87 {
88 #ifndef LIBMESH_BROKEN_IOSTREAM
89 
90  // save the initial format flags
91  std::ios_base::fmtflags os_flags = os.flags();
92 
93  // Print the matrix entries.
94  for (unsigned int i=0; i<this->m(); i++)
95  {
96  for (unsigned int j=0; j<this->n(); j++)
97  os << std::setw(15)
98  << std::scientific
99  << std::setprecision(precision)
100  << this->el(i,j) << " ";
101 
102  os << std::endl;
103  }
104 
105  // reset the original format flags
106  os.flags(os_flags);
107 
108 #else
109 
110  // Print the matrix entries.
111  for (unsigned int i=0; i<this->m(); i++)
112  {
113  for (unsigned int j=0; j<this->n(); j++)
114  os << std::setprecision(precision)
115  << this->el(i,j)
116  << " ";
117 
118  os << std::endl;
119  }
120 
121 
122 #endif
123 }
unsigned int n() const
unsigned int m() const
virtual T el(const unsigned int i, const unsigned int j) const =0
template<typename T >
void libMesh::DenseMatrix< T >::resize ( const unsigned int  new_m,
const unsigned int  new_n 
)
inline

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

Definition at line 776 of file dense_matrix.h.

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

Referenced by libMesh::DenseMatrix< T >::_evd_lapack(), libMesh::DenseMatrix< T >::_svd_lapack(), libMesh::HPCoarsenTest::add_projection(), libMesh::DofMap::build_constraint_matrix(), libMesh::DofMap::build_constraint_matrix_and_vector(), libMesh::FEGenericBase< OutputType >::coarsened_dof_values(), libMesh::FEGenericBase< OutputType >::compute_periodic_constraints(), libMesh::FEGenericBase< OutputType >::compute_proj_constraints(), libMesh::DenseMatrix< T >::DenseMatrix(), libMesh::DenseMatrix< Number >::el(), libMesh::DenseMatrix< T >::get_principal_submatrix(), libMesh::DenseMatrix< T >::get_transpose(), libMesh::FESubdivision::init_subdivision_matrix(), libMesh::DGFEMContext::neighbor_side_fe_reinit(), libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction >::operator()(), libMesh::BoundaryProjectSolution::operator()(), libMesh::DenseMatrix< T >::operator=(), libMesh::FEMContext::pre_fe_reinit(), and libMesh::HPCoarsenTest::select_refinement().

778 {
779  _val.resize(new_m*new_n);
780 
781  this->_m = new_m;
782  this->_n = new_n;
783 
784  // zero and set decomposition_type to NONE
785  this->zero();
786 }
virtual void zero() libmesh_override
Definition: dense_matrix.h:792
std::vector< T > _val
Definition: dense_matrix.h:516
template<typename T>
void libMesh::DenseMatrix< T >::right_multiply ( const DenseMatrixBase< T > &  M3)
virtual

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

Implements libMesh::DenseMatrixBase< T >.

Definition at line 213 of file dense_matrix_impl.h.

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

Referenced by libMesh::DofMap::build_constraint_matrix(), libMesh::DofMap::build_constraint_matrix_and_vector(), libMesh::DofMap::constrain_element_matrix(), libMesh::DofMap::constrain_element_matrix_and_vector(), libMesh::DenseMatrix< Number >::el(), libMesh::DofMap::heterogenously_constrain_element_matrix_and_vector(), and libMesh::FESubdivision::init_shape_functions().

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 }
void _multiply_blas(const DenseMatrixBase< T > &other, _BLAS_Multiply_Flag flag)
void resize(const unsigned int new_m, const unsigned int new_n)
Definition: dense_matrix.h:776
static void multiply(DenseMatrixBase< T > &M1, const DenseMatrixBase< T > &M2, const DenseMatrixBase< T > &M3)
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.

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

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 }
void resize(const unsigned int new_m, const unsigned int new_n)
Definition: dense_matrix.h:776
static void multiply(DenseMatrixBase< T > &M1, const DenseMatrixBase< T > &M2, const DenseMatrixBase< T > &M3)
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.

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

Referenced by libMesh::DenseMatrix< Number >::el().

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 }
unsigned int n() const
void _multiply_blas(const DenseMatrixBase< T > &other, _BLAS_Multiply_Flag flag)
unsigned int m() const
static PetscErrorCode Mat * A
void resize(const unsigned int new_m, const unsigned int new_n)
Definition: dense_matrix.h:776
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.

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

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 n() const
unsigned int m() const
static PetscErrorCode Mat * A
void resize(const unsigned int new_m, const unsigned int new_n)
Definition: dense_matrix.h:776
template<typename T>
void libMesh::DenseMatrix< T >::scale ( const T  factor)
inline

Multiplies every element in the matrix by factor.

Definition at line 851 of file dense_matrix.h.

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

Referenced by libMesh::DenseMatrix< Number >::el(), and libMesh::DenseMatrix< T >::operator*=().

852 {
853  for (std::size_t i=0; i<_val.size(); i++)
854  _val[i] *= factor;
855 }
std::vector< T > _val
Definition: dense_matrix.h:516
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 860 of file dense_matrix.h.

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

Referenced by libMesh::DenseMatrix< Number >::el().

861 {
862  for (unsigned int i=0; i<this->m(); i++)
863  (*this)(i, col) *= factor;
864 }
unsigned int m() const
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 780 of file dense_matrix_impl.h.

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

781 {
782  // We use the LAPACK svd implementation
783  _svd_lapack(sigma);
784 }
void _svd_lapack(DenseVector< Real > &sigma)
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 788 of file dense_matrix_impl.h.

791 {
792  // We use the LAPACK svd implementation
793  _svd_lapack(sigma, U, VT);
794 }
void _svd_lapack(DenseVector< Real > &sigma)
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 799 of file dense_matrix_impl.h.

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

802 {
803  _svd_solve_lapack(rhs, x, rcond);
804 }
void _svd_solve_lapack(const DenseVector< T > &rhs, DenseVector< T > &x, Real rcond) const
template<typename T>
void libMesh::DenseMatrix< T >::swap ( DenseMatrix< T > &  other_matrix)
inline

STL-like swap method

Definition at line 746 of file dense_matrix.h.

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

Referenced by libMesh::EulerSolver::_general_residual(), libMesh::Euler2Solver::_general_residual(), libMesh::NewmarkSolver::_general_residual(), and libMesh::DenseMatrix< Number >::el().

747 {
748  std::swap(this->_m, other_matrix._m);
749  std::swap(this->_n, other_matrix._n);
750  _val.swap(other_matrix._val);
752  _decomposition_type = other_matrix._decomposition_type;
753  other_matrix._decomposition_type = _temp;
754 }
DecompositionType _decomposition_type
Definition: dense_matrix.h:566
void swap(Iterator &lhs, Iterator &rhs)
std::vector< T > _val
Definition: dense_matrix.h:516
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 1044 of file dense_matrix.h.

Referenced by libMesh::DenseMatrix< Number >::el(), libMesh::DenseMatrix< T >::left_multiply_transpose(), and libMesh::DenseMatrix< T >::right_multiply_transpose().

1046 {
1047  // Implement in terms of operator()
1048  return (*this)(j,i);
1049 }
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.

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

Referenced by libMesh::FEMap::compute_single_point_map(), libMesh::DenseMatrix< Number >::el(), libMesh::DofMap::heterogenously_constrain_element_matrix_and_vector(), and libMesh::DofMap::heterogenously_constrain_element_vector().

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 n() const
void _matvec_blas(T alpha, T beta, DenseVector< T > &dest, const DenseVector< T > &arg, bool trans=false) const
unsigned int m() const
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.

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

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 n() const
unsigned int m() const
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.

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

Referenced by libMesh::DofMap::build_constraint_matrix_and_vector(), and libMesh::DenseMatrix< Number >::el().

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 }
void _matvec_blas(T alpha, T beta, DenseVector< T > &dest, const DenseVector< T > &arg, bool trans=false) const
unsigned int m() const
void vector_mult(DenseVector< T > &dest, const DenseVector< T > &arg) const
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.

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

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

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

Referenced by libMesh::DofMap::constrain_element_dyad_matrix(), libMesh::DofMap::constrain_element_matrix_and_vector(), libMesh::DofMap::constrain_element_vector(), libMesh::DenseMatrix< Number >::el(), libMesh::DofMap::heterogenously_constrain_element_matrix_and_vector(), and libMesh::DofMap::heterogenously_constrain_element_vector().

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 n() const
void _matvec_blas(T alpha, T beta, DenseVector< T > &dest, const DenseVector< T > &arg, bool trans=false) const
unsigned int m() const
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.

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

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 n() const
unsigned int m() const

Member Data Documentation

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

Referenced by libMesh::DenseMatrix< T >::operator=(), libMesh::DenseMatrix< T >::swap(), and libMesh::DenseMatrix< T >::zero().

template<typename T>
std::vector<pivot_index_t> libMesh::DenseMatrix< T >::_pivots
private

Definition at line 656 of file dense_matrix.h.

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


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