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.

Author
Benjamin S. Kirk
Date
2002

Definition at line 53 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 660 of file dense_matrix.h.

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

Definition at line 662 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 580 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 568 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 738 of file dense_matrix.h.

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

739  :
740  DenseMatrixBase<T>(new_m,new_n),
741 #if defined(LIBMESH_HAVE_PETSC) && defined(LIBMESH_USE_REAL_NUMBERS) && defined(LIBMESH_DEFAULT_DOUBLE_PRECISION)
742  use_blas_lapack(true),
743 #else
744  use_blas_lapack(false),
745 #endif
746  _val(),
748 {
749  this->resize(new_m,new_n);
750 }
DecompositionType _decomposition_type
Definition: dense_matrix.h:574
void resize(const unsigned int new_m, const unsigned int new_n)
Definition: dense_matrix.h:799
std::vector< T > _val
Definition: dense_matrix.h:526
template<typename T>
virtual libMesh::DenseMatrix< T >::~DenseMatrix ( )
inlinevirtual

Copy-constructor. Destructor. Empty.

Definition at line 71 of file dense_matrix.h.

71 {}

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 that 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 that 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:574
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:526
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:664
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:526
std::vector< pivot_index_t > _pivots
Definition: dense_matrix.h:664
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:574
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:664
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:574
unsigned int m() const
std::vector< T > _val
Definition: dense_matrix.h:526
std::vector< pivot_index_t > _pivots
Definition: dense_matrix.h:664
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:526
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:526
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:338
void resize(const unsigned int new_m, const unsigned int new_n)
Definition: dense_matrix.h:799
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.

Definition at line 906 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().

908 {
909  libmesh_assert_equal_to (this->m(), mat.m());
910  libmesh_assert_equal_to (this->n(), mat.n());
911 
912  for (unsigned int i=0; i<this->m(); i++)
913  for (unsigned int j=0; j<this->n(); j++)
914  (*this)(i,j) += factor * mat(i,j);
915 }
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 that 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:574
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 351 of file dense_matrix.h.

355  { 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 that this means doing an LU decomposition and then computing the product of the diagonal terms. Therefore this is a non-const method.

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:574
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:664
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 as a writeable reference.

Implements libMesh::DenseMatrixBase< T >.

Definition at line 94 of file dense_matrix.h.

96  { 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 writeable reference.

Implements libMesh::DenseMatrixBase< T >.

Definition at line 101 of file dense_matrix.h.

103  { 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 A are overwritten by this function!

The implementation requires the LAPACKgeev_ 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.

Warning: the contents of A 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_ which is wrapped 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 A are overwritten by this function!

See the documentation of the evd_left() and evd_right() functions for more information. The implementation requires the LAPACKgeev_ which is wrapped 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.

Warning: the contents of A 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_ which is wrapped 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
template<typename T>
const std::vector<T>& libMesh::DenseMatrix< T >::get_values ( ) const
inline

Return a constant reference to the matrix values.

Definition at line 343 of file dense_matrix.h.

343 { return _val; }
std::vector< T > _val
Definition: dense_matrix.h:526
template<typename T >
Real libMesh::DenseMatrix< T >::l1_norm ( ) const
inline

Return the l1-norm of the matrix, that is $|M|_1=max_{all columns j}\sum_{all rows i} |M_ij|$, (max. sum of columns). 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 1013 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().

1014 {
1015  libmesh_assert (this->_m);
1016  libmesh_assert (this->_n);
1017 
1018  Real columnsum = 0.;
1019  for (unsigned int i=0; i!=this->_m; i++)
1020  {
1021  columnsum += std::abs((*this)(i,0));
1022  }
1023  Real my_max = columnsum;
1024  for (unsigned int j=1; j!=this->_n; j++)
1025  {
1026  columnsum = 0.;
1027  for (unsigned int i=0; i!=this->_m; i++)
1028  {
1029  columnsum += std::abs((*this)(i,j));
1030  }
1031  my_max = (my_max > columnsum? my_max : columnsum);
1032  }
1033  return my_max;
1034 }
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

Left multipliess by the matrix M2.

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:799
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 multipliess 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:799
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:799
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 of different 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:799
template<typename T >
Real libMesh::DenseMatrix< T >::linfty_norm ( ) const
inline

Return the linfty-norm of the matrix, that is $|M|_\infty=max_{all rows i}\sum_{all columns j} |M_ij|$, (max. sum of rows). 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 1040 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().

1041 {
1042  libmesh_assert (this->_m);
1043  libmesh_assert (this->_n);
1044 
1045  Real rowsum = 0.;
1046  for (unsigned int j=0; j!=this->_n; j++)
1047  {
1048  rowsum += std::abs((*this)(0,j));
1049  }
1050  Real my_max = rowsum;
1051  for (unsigned int i=1; i!=this->_m; i++)
1052  {
1053  rowsum = 0.;
1054  for (unsigned int j=0; j!=this->_n; j++)
1055  {
1056  rowsum += std::abs((*this)(i,j));
1057  }
1058  my_max = (my_max > rowsum? my_max : rowsum);
1059  }
1060  return my_max;
1061 }
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:574
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::SparseMatrix< T >::add_block_matrix(), libMesh::PetscMatrix< 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 element in the matrix. In case of complex numbers, this returns the maximum Real part.

Definition at line 992 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().

993 {
994  libmesh_assert (this->_m);
995  libmesh_assert (this->_n);
996  Real my_max = libmesh_real((*this)(0,0));
997 
998  for (unsigned int i=0; i!=this->_m; i++)
999  {
1000  for (unsigned int j=0; j!=this->_n; j++)
1001  {
1002  Real current = libmesh_real((*this)(i,j));
1003  my_max = (my_max > current? my_max : current);
1004  }
1005  }
1006  return my_max;
1007 }
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 element in the matrix. In case of complex numbers, this returns the minimum Real part.

Definition at line 971 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().

972 {
973  libmesh_assert (this->_m);
974  libmesh_assert (this->_n);
975  Real my_min = libmesh_real((*this)(0,0));
976 
977  for (unsigned int i=0; i!=this->_m; i++)
978  {
979  for (unsigned int j=0; j!=this->_n; j++)
980  {
981  Real current = libmesh_real((*this)(i,j));
982  my_min = (my_min < current? my_min : current);
983  }
984  }
985  return my_min;
986 }
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::SparseMatrix< T >::add_block_matrix(), libMesh::PetscMatrix< 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

Tests if mat is not exactly equal to this matrix.

Definition at line 934 of file dense_matrix.h.

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

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

935 {
936  for (std::size_t i=0; i<_val.size(); i++)
937  if (_val[i] != mat._val[i])
938  return true;
939 
940  return false;
941 }
std::vector< T > _val
Definition: dense_matrix.h:526
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 841 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().

843 {
844  libmesh_assert_less (i*j, _val.size());
845  libmesh_assert_less (i, this->_m);
846  libmesh_assert_less (j, this->_n);
847 
848 
849  // return _val[(i) + (this->_m)*(j)]; // col-major
850  return _val[(i)*(this->_n) + (j)]; // row-major
851 }
std::vector< T > _val
Definition: dense_matrix.h:526
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 writeable reference.

Definition at line 857 of file dense_matrix.h.

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

859 {
860  libmesh_assert_less (i*j, _val.size());
861  libmesh_assert_less (i, this->_m);
862  libmesh_assert_less (j, this->_n);
863 
864  //return _val[(i) + (this->_m)*(j)]; // col-major
865  return _val[(i)*(this->_n) + (j)]; // row-major
866 }
std::vector< T > _val
Definition: dense_matrix.h:526
template<typename T>
DenseMatrix< T > & libMesh::DenseMatrix< T >::operator*= ( const T  factor)
inline

Multiplies every element in the matrix by factor.

Definition at line 893 of file dense_matrix.h.

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

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

894 {
895  this->scale(factor);
896  return *this;
897 }
void scale(const T factor)
Definition: dense_matrix.h:874
template<typename T>
DenseMatrix< T > & libMesh::DenseMatrix< T >::operator+= ( const DenseMatrix< T > &  mat)
inline

Adds mat to this matrix.

Definition at line 947 of file dense_matrix.h.

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

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

948 {
949  for (std::size_t i=0; i<_val.size(); i++)
950  _val[i] += mat._val[i];
951 
952  return *this;
953 }
std::vector< T > _val
Definition: dense_matrix.h:526
template<typename T>
DenseMatrix< T > & libMesh::DenseMatrix< T >::operator-= ( const DenseMatrix< T > &  mat)
inline

Subtracts mat from this matrix.

Definition at line 959 of file dense_matrix.h.

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

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

960 {
961  for (std::size_t i=0; i<_val.size(); i++)
962  _val[i] -= mat._val[i];
963 
964  return *this;
965 }
std::vector< T > _val
Definition: dense_matrix.h:526
template<typename T>
DenseMatrix< T > & libMesh::DenseMatrix< T >::operator= ( const DenseMatrix< T > &  other_matrix)
inline

Assignment operator.

Definition at line 826 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().

827 {
828  this->_m = other_matrix._m;
829  this->_n = other_matrix._n;
830 
831  _val = other_matrix._val;
832  _decomposition_type = other_matrix._decomposition_type;
833 
834  return *this;
835 }
DecompositionType _decomposition_type
Definition: dense_matrix.h:574
std::vector< T > _val
Definition: dense_matrix.h:526
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

Definition at line 784 of file dense_matrix.h.

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

785 {
786  unsigned int mat_m = mat.m(), mat_n = mat.n();
787  this->resize(mat_m, mat_n);
788  for (unsigned int i=0; i<mat_m; i++)
789  for (unsigned int j=0; j<mat_n; j++)
790  (*this)(i,j) = mat(i,j);
791 
792  return *this;
793 }
void resize(const unsigned int new_m, const unsigned int new_n)
Definition: dense_matrix.h:799
template<typename T>
bool libMesh::DenseMatrix< T >::operator== ( const DenseMatrix< T > &  mat) const
inline

Tests if mat is exactly equal to this matrix.

Definition at line 921 of file dense_matrix.h.

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

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

922 {
923  for (std::size_t i=0; i<_val.size(); i++)
924  if (_val[i] != mat._val[i])
925  return false;
926 
927  return true;
928 }
std::vector< T > _val
Definition: dense_matrix.h:526
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 799 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().

801 {
802  _val.resize(new_m*new_n);
803 
804  this->_m = new_m;
805  this->_n = new_n;
806 
807  // zero and set decomposition_type to NONE
808  this->zero();
809 }
virtual void zero() libmesh_override
Definition: dense_matrix.h:815
std::vector< T > _val
Definition: dense_matrix.h:526
template<typename T>
void libMesh::DenseMatrix< T >::right_multiply ( const DenseMatrixBase< T > &  M2)
virtual

Right multiplies by the matrix M2.

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:799
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:799
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:799
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 of different 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:799
template<typename T>
void libMesh::DenseMatrix< T >::scale ( const T  factor)
inline

Multiplies every element in the matrix by factor.

Definition at line 874 of file dense_matrix.h.

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

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

875 {
876  for (std::size_t i=0; i<_val.size(); i++)
877  _val[i] *= factor;
878 }
std::vector< T > _val
Definition: dense_matrix.h:526
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 883 of file dense_matrix.h.

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

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

884 {
885  for (unsigned int i=0; i<this->m(); i++)
886  (*this)(i, col) *= factor;
887 }
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 769 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().

770 {
771  std::swap(this->_m, other_matrix._m);
772  std::swap(this->_n, other_matrix._n);
773  _val.swap(other_matrix._val);
775  _decomposition_type = other_matrix._decomposition_type;
776  other_matrix._decomposition_type = _temp;
777 }
DecompositionType _decomposition_type
Definition: dense_matrix.h:574
void swap(Iterator &lhs, Iterator &rhs)
std::vector< T > _val
Definition: dense_matrix.h:526
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 1067 of file dense_matrix.h.

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

1069 {
1070  // Implement in terms of operator()
1071  return (*this)(j,i);
1072 }
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 574 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 664 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 available on the web. This is not the most memory efficient routine since the inverse is not computed "in place" but is instead placed into a the matrix inv passed in by the user. Run-time selectable option to turn on/off blas support. This was primarily used for testing purposes, and could be removed...

Definition at line 519 of file dense_matrix.h.


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