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 >
template 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 987 of file dense_matrix.C.

References A, libMesh::DenseMatrix< T >::cholesky_solve(), libMesh::LIBMESH_VMA_INSTANTIATE(), libMesh::Real, libMesh::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::vector_mult(), libMesh::DenseMatrix< T >::vector_mult_transpose(), and libMesh::x.

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

References A.

942 {
943  // If we called this function, there better not be any
944  // previous decomposition of the matrix.
945  libmesh_assert_equal_to (this->_decomposition_type, NONE);
946 
947  // Shorthand notation for number of rows and columns.
948  const unsigned int
949  n_rows = this->m(),
950  n_cols = this->n();
951 
952  // Just to be really sure...
953  libmesh_assert_equal_to (n_rows, n_cols);
954 
955  // A convenient reference to *this
956  DenseMatrix<T> & A = *this;
957 
958  for (unsigned int i=0; i<n_rows; ++i)
959  {
960  for (unsigned int j=i; j<n_cols; ++j)
961  {
962  for (unsigned int k=0; k<i; ++k)
963  A(i,j) -= A(i,k) * A(j,k);
964 
965  if (i == j)
966  {
967 #ifndef LIBMESH_USE_COMPLEX_NUMBERS
968  if (A(i,j) <= 0.0)
969  libmesh_error_msg("Error! Can only use Cholesky decomposition with symmetric positive definite matrices.");
970 #endif
971 
972  A(i,i) = std::sqrt(A(i,j));
973  }
974  else
975  A(j,i) = A(i,j) / A(i,i);
976  }
977  }
978 
979  // Set the flag for CHOLESKY decomposition
981 }
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 700 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().

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

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

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

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

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

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

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

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

References std::min().

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

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

1004 {
1005  // Ensure that dest and arg sizes are compatible
1006  if (!trans)
1007  {
1008  // dest ~ A * arg
1009  // (mx1) (mxn) * (nx1)
1010  if ((dest.size() != this->m()) || (arg.size() != this->n()))
1011  libmesh_error_msg("Improper input argument sizes!");
1012  }
1013 
1014  else // trans == true
1015  {
1016  // Ensure that dest and arg are proper size
1017  // dest ~ A^T * arg
1018  // (nx1) (nxm) * (mx1)
1019  if ((dest.size() != this->n()) || (arg.size() != this->m()))
1020  libmesh_error_msg("Improper input argument sizes!");
1021  }
1022 
1023  // Calling sequence for dgemv:
1024  //
1025  // dgemv(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
1026 
1027  // TRANS (input)
1028  // 't' for transpose, 'n' for non-transpose multiply
1029  // We store everything in row-major order, so pass the transpose flag for
1030  // non-transposed matvecs and the 'n' flag for transposed matvecs
1031  char TRANS[] = "t";
1032  if (trans)
1033  TRANS[0] = 'n';
1034 
1035  // M (input)
1036  // On entry, M specifies the number of rows of the matrix A.
1037  // In C/C++, pass the number of *cols* of A
1038  PetscBLASInt M = this->n();
1039 
1040  // N (input)
1041  // On entry, N specifies the number of columns of the matrix A.
1042  // In C/C++, pass the number of *rows* of A
1043  PetscBLASInt N = this->m();
1044 
1045  // ALPHA (input)
1046  // The scalar constant passed to this function
1047 
1048  // A (input) double precision array of DIMENSION ( LDA, n ).
1049  // Before entry, the leading m by n part of the array A must
1050  // contain the matrix of coefficients.
1051  // The matrix, *this. Note that _matvec_blas is called from
1052  // a const function, vector_mult(), and so we have made this function const
1053  // as well. Since BLAS knows nothing about const, we have to cast it away
1054  // now.
1055  DenseMatrix<T> & a_ref = const_cast< DenseMatrix<T> &> ( *this );
1056  std::vector<T> & a = a_ref.get_values();
1057 
1058  // LDA (input)
1059  // On entry, LDA specifies the first dimension of A as declared
1060  // in the calling (sub) program. LDA must be at least
1061  // max( 1, m ).
1062  PetscBLASInt LDA = M;
1063 
1064  // X (input) double precision array of DIMENSION at least
1065  // ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
1066  // and at least
1067  // ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
1068  // Before entry, the incremented array X must contain the
1069  // vector x.
1070  // Here, we must cast away the const-ness of "arg" since BLAS knows
1071  // nothing about const
1072  DenseVector<T> & x_ref = const_cast< DenseVector<T> &> ( arg );
1073  std::vector<T> & x = x_ref.get_values();
1074 
1075  // INCX (input)
1076  // On entry, INCX specifies the increment for the elements of
1077  // X. INCX must not be zero.
1078  PetscBLASInt INCX = 1;
1079 
1080  // BETA (input)
1081  // On entry, BETA specifies the scalar beta. When BETA is
1082  // supplied as zero then Y need not be set on input.
1083  // The second scalar constant passed to this function
1084 
1085  // Y (input) double precision array of DIMENSION at least
1086  // ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
1087  // and at least
1088  // ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
1089  // Before entry with BETA non-zero, the incremented array Y
1090  // must contain the vector y. On exit, Y is overwritten by the
1091  // updated vector y.
1092  // The input vector "dest"
1093  std::vector<T> & y = dest.get_values();
1094 
1095  // INCY (input)
1096  // On entry, INCY specifies the increment for the elements of
1097  // Y. INCY must not be zero.
1098  PetscBLASInt INCY = 1;
1099 
1100  // Finally, ready to call the BLAS function
1101  BLASgemv_(TRANS, &M, &N, &alpha, &(a[0]), &LDA, &(x[0]), &INCX, &beta, &(y[0]), &INCY);
1102 }
unsigned int n() const
PetscErrorCode Vec x
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_tranpose() 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  case RIGHT_MULTIPLY:
51  {
52  result_size = other.n() * this->m();
53  if (other.m() == this->n())
54  break;
55  }
57  {
58  result_size = other.n() * this->n();
59  if (other.m() == this->m())
60  break;
61  }
63  {
64  result_size = other.m() * this->m();
65  if (other.n() == this->n())
66  break;
67  }
68  default:
69  libmesh_error_msg("Unknown flag selected or matrices are incompatible for multiplication.");
70  }
71 
72  // For this to work, the passed arg. must actually be a DenseMatrix<T>
73  const DenseMatrix<T> * const_that = cast_ptr<const DenseMatrix<T> *>(&other);
74 
75  // Also, although 'that' is logically const in this BLAS routine,
76  // the PETSc BLAS interface does not specify that any of the inputs are
77  // const. To use it, I must cast away const-ness.
78  DenseMatrix<T> * that = const_cast< DenseMatrix<T> * > (const_that);
79 
80  // Initialize A, B pointers for LEFT_MULTIPLY* cases
81  DenseMatrix<T> * A = this;
82  DenseMatrix<T> * B = that;
83 
84  // For RIGHT_MULTIPLY* cases, swap the meaning of A and B.
85  // Here is a full table of combinations we can pass to BLASgemm, and what the answer is when finished:
86  // pass A B -> (Fortran) -> A^T B^T -> (C++) -> (A^T B^T)^T -> (identity) -> B A "lt multiply"
87  // pass B A -> (Fortran) -> B^T A^T -> (C++) -> (B^T A^T)^T -> (identity) -> A B "rt multiply"
88  // pass A B^T -> (Fortran) -> A^T B -> (C++) -> (A^T B)^T -> (identity) -> B^T A "lt multiply t"
89  // pass B^T A -> (Fortran) -> B A^T -> (C++) -> (B A^T)^T -> (identity) -> A B^T "rt multiply t"
90  if (flag==RIGHT_MULTIPLY || flag==RIGHT_MULTIPLY_TRANSPOSE)
91  std::swap(A,B);
92 
93  // transa, transb values to pass to blas
94  char
95  transa[] = "n",
96  transb[] = "n";
97 
98  // Integer values to pass to BLAS:
99  //
100  // M
101  // In Fortran, the number of rows of op(A),
102  // In the BLAS documentation, typically known as 'M'.
103  //
104  // In C/C++, we set:
105  // M = n_cols(A) if (transa='n')
106  // n_rows(A) if (transa='t')
107  PetscBLASInt M = static_cast<PetscBLASInt>( A->n() );
108 
109  // N
110  // In Fortran, the number of cols of op(B), and also the number of cols of C.
111  // In the BLAS documentation, typically known as 'N'.
112  //
113  // In C/C++, we set:
114  // N = n_rows(B) if (transb='n')
115  // n_cols(B) if (transb='t')
116  PetscBLASInt N = static_cast<PetscBLASInt>( B->m() );
117 
118  // K
119  // In Fortran, the number of cols of op(A), and also
120  // the number of rows of op(B). In the BLAS documentation,
121  // typically known as 'K'.
122  //
123  // In C/C++, we set:
124  // K = n_rows(A) if (transa='n')
125  // n_cols(A) if (transa='t')
126  PetscBLASInt K = static_cast<PetscBLASInt>( A->m() );
127 
128  // LDA (leading dimension of A). In our cases,
129  // LDA is always the number of columns of A.
130  PetscBLASInt LDA = static_cast<PetscBLASInt>( A->n() );
131 
132  // LDB (leading dimension of B). In our cases,
133  // LDB is always the number of columns of B.
134  PetscBLASInt LDB = static_cast<PetscBLASInt>( B->n() );
135 
136  if (flag == LEFT_MULTIPLY_TRANSPOSE)
137  {
138  transb[0] = 't';
139  N = static_cast<PetscBLASInt>( B->n() );
140  }
141 
142  else if (flag == RIGHT_MULTIPLY_TRANSPOSE)
143  {
144  transa[0] = 't';
145  std::swap(M,K);
146  }
147 
148  // LDC (leading dimension of C). LDC is the
149  // number of columns in the solution matrix.
150  PetscBLASInt LDC = M;
151 
152  // Scalar values to pass to BLAS
153  //
154  // scalar multiplying the whole product AB
155  T alpha = 1.;
156 
157  // scalar multiplying C, which is the original matrix.
158  T beta = 0.;
159 
160  // Storage for the result
161  std::vector<T> result (result_size);
162 
163  // Finally ready to call the BLAS
164  BLASgemm_(transa, transb, &M, &N, &K, &alpha, &(A->_val[0]), &LDA, &(B->_val[0]), &LDB, &beta, &result[0], &LDC);
165 
166  // Update the relevant dimension for this matrix.
167  switch (flag)
168  {
169  case LEFT_MULTIPLY: { this->_m = other.m(); break; }
170  case RIGHT_MULTIPLY: { this->_n = other.n(); break; }
171  case LEFT_MULTIPLY_TRANSPOSE: { this->_m = other.n(); break; }
172  case RIGHT_MULTIPLY_TRANSPOSE: { this->_n = other.m(); break; }
173  default:
174  libmesh_error_msg("Unknown flag selected.");
175  }
176 
177  // Swap my data vector with the result
178  this->_val.swap(result);
179 }
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 377 of file dense_matrix_blas_lapack.C.

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

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

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

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

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

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

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

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

Referenced by libMesh::DenseMatrix< T >::_cholesky_back_substitute(), 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().

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

References libMesh::Real.

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

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

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

806 {
807  // We use the LAPACK eigenvalue problem implementation
808  _evd_lapack(lambda_real, lambda_imag);
809 }
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 814 of file dense_matrix.C.

References libmesh_nullptr.

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

817 {
818  // We use the LAPACK eigenvalue problem implementation
819  _evd_lapack(lambda_real, lambda_imag, &VL, libmesh_nullptr);
820 }
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 836 of file dense_matrix.C.

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

840 {
841  // We use the LAPACK eigenvalue problem implementation
842  _evd_lapack(lambda_real, lambda_imag, &VL, &VR);
843 }
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 825 of file dense_matrix.C.

References libmesh_nullptr.

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

828 {
829  // We use the LAPACK eigenvalue problem implementation
830  _evd_lapack(lambda_real, lambda_imag, libmesh_nullptr, &VR);
831 }
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 553 of file dense_matrix.C.

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

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

556 {
557  libmesh_assert( (sub_m <= this->m()) && (sub_n <= this->n()) );
558 
559  dest.resize(sub_m, sub_n);
560  for (unsigned int i=0; i<sub_m; i++)
561  for (unsigned int j=0; j<sub_n; j++)
562  dest(i,j) = (*this)(i,j);
563 }
unsigned int n() const
libmesh_assert(j)
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 568 of file dense_matrix.C.

569 {
570  get_principal_submatrix(sub_m, sub_m, dest);
571 }
void get_principal_submatrix(unsigned int sub_m, unsigned int sub_n, DenseMatrix< T > &dest) const
Definition: dense_matrix.C:553
template<typename T>
void libMesh::DenseMatrix< T >::get_transpose ( DenseMatrix< T > &  dest) const

Put the tranposed matrix into dest.

Definition at line 576 of file dense_matrix.C.

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

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

577 {
578  dest.resize(this->n(), this->m());
579 
580  for (unsigned int i=0; i<dest.m(); i++)
581  for (unsigned int j=0; j<dest.n(); j++)
582  dest(i,j) = (*this)(j,i);
583 }
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(), libMesh::libmesh_assert(), 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)
libmesh_assert(j)
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 37 of file dense_matrix.C.

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

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

38 {
39  if (this->use_blas_lapack)
40  this->_multiply_blas(M2, LEFT_MULTIPLY);
41  else
42  {
43  // (*this) <- M2 * (*this)
44  // Where:
45  // (*this) = (m x n),
46  // M2 = (m x p),
47  // M3 = (p x n)
48 
49  // M3 is a copy of *this before it gets resize()d
50  DenseMatrix<T> M3(*this);
51 
52  // Resize *this so that the result can fit
53  this->resize (M2.m(), M3.n());
54 
55  // Call the multiply function in the base class
56  this->multiply(*this, M2, M3);
57  }
58 }
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 64 of file dense_matrix.C.

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

65 {
66  // (*this) <- M2 * (*this)
67  // Where:
68  // (*this) = (m x n),
69  // M2 = (m x p),
70  // M3 = (p x n)
71 
72  // M3 is a copy of *this before it gets resize()d
73  DenseMatrix<T> M3(*this);
74 
75  // Resize *this so that the result can fit
76  this->resize (M2.m(), M3.n());
77 
78  // Call the multiply function in the base class
79  this->multiply(*this, M2, M3);
80 }
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 85 of file dense_matrix.C.

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

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

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

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

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

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

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

211 {
212  if (this->use_blas_lapack)
213  this->_multiply_blas(M3, RIGHT_MULTIPLY);
214  else
215  {
216  // (*this) <- M3 * (*this)
217  // Where:
218  // (*this) = (m x n),
219  // M2 = (m x p),
220  // M3 = (p x n)
221 
222  // M2 is a copy of *this before it gets resize()d
223  DenseMatrix<T> M2(*this);
224 
225  // Resize *this so that the result can fit
226  this->resize (M2.m(), M3.n());
227 
228  this->multiply(*this, M2, M3);
229  }
230 }
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 236 of file dense_matrix.C.

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

237 {
238  // (*this) <- M3 * (*this)
239  // Where:
240  // (*this) = (m x n),
241  // M2 = (m x p),
242  // M3 = (p x n)
243 
244  // M2 is a copy of *this before it gets resize()d
245  DenseMatrix<T> M2(*this);
246 
247  // Resize *this so that the result can fit
248  this->resize (M2.m(), M3.n());
249 
250  this->multiply(*this, M2, M3);
251 }
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 257 of file dense_matrix.C.

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

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

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

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

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

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

776 {
777  // We use the LAPACK svd implementation
778  _svd_lapack(sigma);
779 }
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 783 of file dense_matrix.C.

786 {
787  // We use the LAPACK svd implementation
788  _svd_lapack(sigma, U, VT);
789 }
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 794 of file dense_matrix.C.

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

797 {
798  _svd_solve_lapack(rhs, x, rcond);
799 }
PetscErrorCode Vec x
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 382 of file dense_matrix.C.

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

Referenced by libMesh::DenseMatrix< T >::_cholesky_back_substitute(), 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().

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

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

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

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

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

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

536 {
537  // Short-circuit if the matrix is empty
538  if (this->m() == 0)
539  {
540  dest.resize(0);
541  return;
542  }
543 
544  DenseVector<typename CompareTypes<T,T3>::supertype>
545  temp(arg.size());
546  this->vector_mult(temp, arg);
547  dest.add(factor, temp);
548 }
unsigned int m() const
void vector_mult(DenseVector< T > &dest, const DenseVector< T > &arg) const
Definition: dense_matrix.C:382
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 438 of file dense_matrix.C.

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

Referenced by libMesh::DenseMatrix< T >::_cholesky_back_substitute(), 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().

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

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

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