libMesh::DenseMatrix< T > Class Template Reference

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

#include <dense_matrix.h>

Inheritance diagram for libMesh::DenseMatrix< T >:

Public Member Functions

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

Public Attributes

bool use_blas_lapack
 

Protected Member Functions

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

Static Protected Member Functions

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

Protected Attributes

unsigned int _m
 
unsigned int _n
 

Private Types

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

Private Member Functions

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

Private Attributes

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

Detailed Description

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

A matrix object used for finite element assembly and numerics.

Defines a dense matrix for use in Finite Element-type computations. Useful for storing element stiffness matrices before summation into a global matrix. All overridden virtual functions are documented in dense_matrix_base.h.

Author
Benjamin S. Kirk
Date
2002

Definition at line 54 of file dense_matrix.h.

Member Typedef Documentation

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

Array used to store pivot indices. May be used by whatever factorization is currently active, clients of the class should not rely on it for any reason.

Definition at line 652 of file dense_matrix.h.

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

Definition at line 654 of file dense_matrix.h.

Member Enumeration Documentation

template<typename T>
enum libMesh::DenseMatrix::_BLAS_Multiply_Flag
private

Enumeration used to determine the behavior of the _multiply_blas function.

Enumerator
LEFT_MULTIPLY 
RIGHT_MULTIPLY 
LEFT_MULTIPLY_TRANSPOSE 
RIGHT_MULTIPLY_TRANSPOSE 

Definition at line 572 of file dense_matrix.h.

template<typename T>
enum libMesh::DenseMatrix::DecompositionType
private

The decomposition schemes above change the entries of the matrix A. It is therefore an error to call A.lu_solve() and subsequently call A.cholesky_solve() since the result will probably not match any desired outcome. This typedef keeps track of which decomposition has been called for this matrix.

Enumerator
LU 
CHOLESKY 
LU_BLAS_LAPACK 
NONE 

Definition at line 560 of file dense_matrix.h.

Constructor & Destructor Documentation

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

Constructor. Creates a dense matrix of dimension m by n.

Definition at line 728 of file dense_matrix.h.

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

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

Destructor. Empty.

Definition at line 67 of file dense_matrix.h.

67 {}

Member Function Documentation

template<typename T >
template<typename T2 >
template void libMesh::DenseMatrix< T >::_cholesky_back_substitute ( const DenseVector< T2 > &  b,
DenseVector< T2 > &  x 
) const
private

Solves the equation Ax=b for the unknown value x and rhs b based on the Cholesky factorization of A.

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

Definition at line 989 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.

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

Decomposes a symmetric positive definite matrix into a product of two lower triangular matrices according to A = LL^T.

Note
This program generates an error if the matrix is not SPD.

Definition at line 943 of file dense_matrix.C.

References A.

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

Computes the eigenvalues of the matrix using the Lapack routine "DGEEV". If VR and/or VL are non-NULL, then the matrix of right and/or left eigenvectors is also computed and returned by this function.

[ Implementation in dense_matrix_blas_lapack.C ]

Definition at line 704 of file dense_matrix_blas_lapack.C.

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

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

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

Solves the system Ax=b through back substitution. This function is private since it is only called as part of the implementation of the lu_solve(...) function.

Definition at line 648 of file dense_matrix.C.

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

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

Companion function to _lu_decompose_lapack(). Do not use directly, called through the public lu_solve() interface. This function is logically const in that it does not modify the matrix, but since we are just calling LAPACK routines, it's less const_cast hassle to just declare the function non-const. [ Implementation in dense_matrix_blas_lapack.C ]

Definition at line 913 of file dense_matrix_blas_lapack.C.

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

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

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

Form the LU decomposition of the matrix. This function is private since it is only called as part of the implementation of the lu_solve(...) function.

Definition at line 699 of file dense_matrix.C.

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

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

Computes an LU factorization of the matrix using the Lapack routine "getrf". This routine should only be used by the "use_blas_lapack" branch of the lu_solve() function. After the call to this function, the matrix is replaced by its factorized version, and the DecompositionType is set to LU_BLAS_LAPACK. [ Implementation in dense_matrix_blas_lapack.C ]

Definition at line 205 of file dense_matrix_blas_lapack.C.

References std::min().

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

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

Uses the BLAS GEMV function (through PETSc) to compute

dest := alpha*A*arg + beta*dest

where alpha and beta are scalars, A is this matrix, and arg and dest are input vectors of appropriate size. If trans is true, the transpose matvec is computed instead. By default, trans==false.

[ Implementation in dense_matrix_blas_lapack.C ]

Definition at line 1003 of file dense_matrix_blas_lapack.C.

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

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

Definition at line 35 of file dense_matrix_blas_lapack.C.

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

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

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

Helper function that actually performs the SVD. [ Implementation in dense_matrix_blas_lapack.C ]

Definition at line 381 of file dense_matrix_blas_lapack.C.

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

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

Computes an SVD of the matrix using the Lapack routine "getsvd". [ Implementation in dense_matrix_blas_lapack.C ]

Definition at line 273 of file dense_matrix_blas_lapack.C.

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

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

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

Computes a "reduced" SVD of the matrix using the Lapack routine "getsvd". [ Implementation in dense_matrix_blas_lapack.C ]

Definition at line 315 of file dense_matrix_blas_lapack.C.

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

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

Called by svd_solve(rhs).

Definition at line 527 of file dense_matrix_blas_lapack.C.

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

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

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

Adds factor to every element in the matrix. This should only work if T += T2 * T3 is valid C++ and if T2 is scalar. Return type is void

Definition at line 183 of file dense_matrix_base.h.

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

185 {
186  libmesh_assert_equal_to (this->m(), mat.m());
187  libmesh_assert_equal_to (this->n(), mat.n());
188 
189  for (unsigned int j=0; j<this->n(); j++)
190  for (unsigned int i=0; i<this->m(); i++)
191  this->el(i,j) += factor*mat.el(i,j);
192 }
unsigned int n() const
unsigned int m() const
virtual T el(const unsigned int i, const unsigned int j) const =0
template<typename T >
template<typename T2 , typename T3 >
boostcopy::enable_if_c< ScalarTraits< T2 >::value, void >::type libMesh::DenseMatrix< T >::add ( const T2  factor,
const DenseMatrix< T3 > &  mat 
)
inline

Adds factor times mat to this matrix.

Returns
A reference to *this.

Definition at line 883 of file dense_matrix.h.

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

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

885 {
886  libmesh_assert_equal_to (this->m(), mat.m());
887  libmesh_assert_equal_to (this->n(), mat.n());
888 
889  for (unsigned int i=0; i<this->m(); i++)
890  for (unsigned int j=0; j<this->n(); j++)
891  (*this)(i,j) += factor * mat(i,j);
892 }
unsigned int n() const
unsigned int m() const
template<typename T >
template<typename T2 >
template void libMesh::DenseMatrix< T >::cholesky_solve ( const DenseVector< T2 > &  b,
DenseVector< T2 > &  x 
)

For symmetric positive definite (SPD) matrices. A Cholesky factorization of A such that A = L L^T is about twice as fast as a standard LU factorization. Therefore you can use this method if you know a-priori that the matrix is SPD. If the matrix is not SPD, an error is generated. One nice property of Cholesky decompositions is that they do not require pivoting for stability.

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

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

913 {
914  // Check for a previous decomposition
915  switch(this->_decomposition_type)
916  {
917  case NONE:
918  {
919  this->_cholesky_decompose ();
920  break;
921  }
922 
923  case CHOLESKY:
924  {
925  // Already factored, just need to call back_substitute.
926  break;
927  }
928 
929  default:
930  libmesh_error_msg("Error! This matrix already has a different decomposition...");
931  }
932 
933  // Perform back substitution
934  this->_cholesky_back_substitute (b, x);
935 }
void _cholesky_back_substitute(const DenseVector< T2 > &b, DenseVector< T2 > &x) const
Definition: dense_matrix.C:989
DecompositionType _decomposition_type
Definition: dense_matrix.h:566
PetscErrorCode Vec x
template<typename T >
void libMesh::DenseMatrixBase< T >::condense ( const unsigned int  i,
const unsigned int  j,
const T  val,
DenseVectorBase< T > &  rhs 
)
protectedinherited

Condense-out the (i,j) entry of the matrix, forcing it to take on the value val. This is useful in numerical simulations for applying boundary conditions. Preserves the symmetry of the matrix.

Definition at line 58 of file dense_matrix_base.C.

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

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

62 {
63  libmesh_assert_equal_to (this->_m, rhs.size());
64  libmesh_assert_equal_to (iv, jv);
65 
66 
67  // move the known value into the RHS
68  // and zero the column
69  for (unsigned int i=0; i<this->m(); i++)
70  {
71  rhs.el(i) -= this->el(i,jv)*val;
72  this->el(i,jv) = 0.;
73  }
74 
75  // zero the row
76  for (unsigned int j=0; j<this->n(); j++)
77  this->el(iv,j) = 0.;
78 
79  this->el(iv,jv) = 1.;
80  rhs.el(iv) = val;
81 
82 }
unsigned int n() const
unsigned int m() const
virtual T el(const unsigned int i, const unsigned int j) const =0
template<typename T>
void libMesh::DenseMatrix< T >::condense ( const unsigned int  i,
const unsigned int  j,
const T  val,
DenseVector< T > &  rhs 
)
inline

Condense-out the (i,j) entry of the matrix, forcing it to take on the value val. This is useful in numerical simulations for applying boundary conditions. Preserves the symmetry of the matrix.

Definition at line 338 of file dense_matrix.h.

342  { DenseMatrixBase<T>::condense (i, j, val, rhs); }
void condense(const unsigned int i, const unsigned int j, const T val, DenseVectorBase< T > &rhs)
template<typename T >
T libMesh::DenseMatrix< T >::det ( )
Returns
The determinant of the matrix.
Note
Implemented by computing an LU decomposition and then taking the product of the diagonal terms. Therefore this is a non-const method which modifies the entries of the matrix.

Definition at line 850 of file dense_matrix.C.

References libMesh::Real.

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

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

Implements libMesh::DenseMatrixBase< T >.

Definition at line 84 of file dense_matrix.h.

86  { return (*this)(i,j); }
template<typename T>
virtual T& libMesh::DenseMatrix< T >::el ( const unsigned int  i,
const unsigned int  j 
)
inlinevirtual
Returns
The (i,j) element of the matrix as a writable reference. Since internal data representations may differ, you must redefine this function.

Implements libMesh::DenseMatrixBase< T >.

Definition at line 88 of file dense_matrix.h.

90  { return (*this)(i,j); }
template<typename T>
void libMesh::DenseMatrix< T >::evd ( DenseVector< T > &  lambda_real,
DenseVector< T > &  lambda_imag 
)

Compute the eigenvalues (both real and imaginary parts) of a general matrix.

Warning: the contents of *this are overwritten by this function!

The implementation requires the LAPACKgeev_ function which is wrapped by PETSc.

Definition at line 806 of file dense_matrix.C.

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

808 {
809  // We use the LAPACK eigenvalue problem implementation
810  _evd_lapack(lambda_real, lambda_imag);
811 }
void _evd_lapack(DenseVector< T > &lambda_real, DenseVector< T > &lambda_imag, DenseMatrix< T > *VL=libmesh_nullptr, DenseMatrix< T > *VR=libmesh_nullptr)
template<typename T>
void libMesh::DenseMatrix< T >::evd_left ( DenseVector< T > &  lambda_real,
DenseVector< T > &  lambda_imag,
DenseMatrix< T > &  VL 
)

Compute the eigenvalues (both real and imaginary parts) and left eigenvectors of a general matrix, $ A $.

Warning: the contents of *this are overwritten by this function!

The left eigenvector $ u_j $ of $ A $ satisfies: $ u_j^H A = lambda_j u_j^H $ where $ u_j^H $ denotes the conjugate-transpose of $ u_j $.

If the j-th and (j+1)-st eigenvalues form a complex conjugate pair, then the j-th and (j+1)-st columns of VL "share" their real-valued storage in the following way: u_j = VL(:,j) + i*VL(:,j+1) and u_{j+1} = VL(:,j) - i*VL(:,j+1).

The implementation requires the LAPACKgeev_ routine which is provided by PETSc.

Definition at line 816 of file dense_matrix.C.

References libmesh_nullptr.

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

819 {
820  // We use the LAPACK eigenvalue problem implementation
821  _evd_lapack(lambda_real, lambda_imag, &VL, libmesh_nullptr);
822 }
const class libmesh_nullptr_t libmesh_nullptr
void _evd_lapack(DenseVector< T > &lambda_real, DenseVector< T > &lambda_imag, DenseMatrix< T > *VL=libmesh_nullptr, DenseMatrix< T > *VR=libmesh_nullptr)
template<typename T>
void libMesh::DenseMatrix< T >::evd_left_and_right ( DenseVector< T > &  lambda_real,
DenseVector< T > &  lambda_imag,
DenseMatrix< T > &  VL,
DenseMatrix< T > &  VR 
)

Compute the eigenvalues (both real and imaginary parts) as well as the left and right eigenvectors of a general matrix.

Warning: the contents of *this are overwritten by this function!

See the documentation of the evd_left() and evd_right() functions for more information. The implementation requires the LAPACKgeev_ routine which is provided by PETSc.

Definition at line 838 of file dense_matrix.C.

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

842 {
843  // We use the LAPACK eigenvalue problem implementation
844  _evd_lapack(lambda_real, lambda_imag, &VL, &VR);
845 }
void _evd_lapack(DenseVector< T > &lambda_real, DenseVector< T > &lambda_imag, DenseMatrix< T > *VL=libmesh_nullptr, DenseMatrix< T > *VR=libmesh_nullptr)
template<typename T>
void libMesh::DenseMatrix< T >::evd_right ( DenseVector< T > &  lambda_real,
DenseVector< T > &  lambda_imag,
DenseMatrix< T > &  VR 
)

Compute the eigenvalues (both real and imaginary parts) and right eigenvectors of a general matrix, $ A $.

Warning: the contents of *this are overwritten by this function!

The right eigenvector $ v_j $ of $ A $ satisfies: $ A v_j = lambda_j v_j $ where $ lambda_j $ is its corresponding eigenvalue.

Note
If the j-th and (j+1)-st eigenvalues form a complex conjugate pair, then the j-th and (j+1)-st columns of VR "share" their real-valued storage in the following way: v_j = VR(:,j) + i*VR(:,j+1) and v_{j+1} = VR(:,j) - i*VR(:,j+1).

The implementation requires the LAPACKgeev_ routine which is provided by PETSc.

Definition at line 827 of file dense_matrix.C.

References libmesh_nullptr.

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

830 {
831  // We use the LAPACK eigenvalue problem implementation
832  _evd_lapack(lambda_real, lambda_imag, libmesh_nullptr, &VR);
833 }
const class libmesh_nullptr_t libmesh_nullptr
void _evd_lapack(DenseVector< T > &lambda_real, DenseVector< T > &lambda_imag, DenseMatrix< T > *VL=libmesh_nullptr, DenseMatrix< T > *VR=libmesh_nullptr)
template<typename T>
void libMesh::DenseMatrix< T >::get_principal_submatrix ( unsigned int  sub_m,
unsigned int  sub_n,
DenseMatrix< T > &  dest 
) const

Put the sub_m x sub_n principal submatrix into dest.

Definition at line 553 of file dense_matrix.C.

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

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

556 {
557  libmesh_assert( (sub_m <= this->m()) && (sub_n <= this->n()) );
558 
559  dest.resize(sub_m, sub_n);
560  for (unsigned int i=0; i<sub_m; i++)
561  for (unsigned int j=0; j<sub_n; j++)
562  dest(i,j) = (*this)(i,j);
563 }
unsigned int n() const
libmesh_assert(j)
unsigned int m() const
template<typename T>
void libMesh::DenseMatrix< T >::get_principal_submatrix ( unsigned int  sub_m,
DenseMatrix< T > &  dest 
) const

Put the sub_m x sub_m principal submatrix into dest.

Definition at line 568 of file dense_matrix.C.

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

Put the tranposed matrix into dest.

Definition at line 576 of file dense_matrix.C.

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

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

577 {
578  dest.resize(this->n(), this->m());
579 
580  for (unsigned int i=0; i<dest.m(); i++)
581  for (unsigned int j=0; j<dest.n(); j++)
582  dest(i,j) = (*this)(j,i);
583 }
unsigned int n() const
unsigned int m() const
template<typename T>
std::vector<T>& libMesh::DenseMatrix< T >::get_values ( )
inline
Returns
A reference to the underlying data storage vector.

This should be used with caution (i.e. one should not change the size of the vector, etc.) but is useful for interoperating with low level BLAS routines which expect a simple array.

Definition at line 325 of file dense_matrix.h.

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

325 { return _val; }
std::vector< T > _val
Definition: dense_matrix.h:516
template<typename T>
const std::vector<T>& libMesh::DenseMatrix< T >::get_values ( ) const
inline
Returns
A constant reference to the underlying data storage vector.

Definition at line 330 of file dense_matrix.h.

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

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

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

Definition at line 990 of file dense_matrix.h.

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

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

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

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

Implements libMesh::DenseMatrixBase< T >.

Definition at line 37 of file dense_matrix.C.

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

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

38 {
39  if (this->use_blas_lapack)
40  this->_multiply_blas(M2, LEFT_MULTIPLY);
41  else
42  {
43  // (*this) <- M2 * (*this)
44  // Where:
45  // (*this) = (m x n),
46  // M2 = (m x p),
47  // M3 = (p x n)
48 
49  // M3 is a copy of *this before it gets resize()d
50  DenseMatrix<T> M3(*this);
51 
52  // Resize *this so that the result can fit
53  this->resize (M2.m(), M3.n());
54 
55  // Call the multiply function in the base class
56  this->multiply(*this, M2, M3);
57  }
58 }
void _multiply_blas(const DenseMatrixBase< T > &other, _BLAS_Multiply_Flag flag)
void resize(const unsigned int new_m, const unsigned int new_n)
Definition: dense_matrix.h:776
static void multiply(DenseMatrixBase< T > &M1, const DenseMatrixBase< T > &M2, const DenseMatrixBase< T > &M3)
template<typename T >
template<typename T2 >
void libMesh::DenseMatrix< T >::left_multiply ( const DenseMatrixBase< T2 > &  M2)

Left multiplies by the matrix M2 of different type

Definition at line 64 of file dense_matrix.C.

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

65 {
66  // (*this) <- M2 * (*this)
67  // Where:
68  // (*this) = (m x n),
69  // M2 = (m x p),
70  // M3 = (p x n)
71 
72  // M3 is a copy of *this before it gets resize()d
73  DenseMatrix<T> M3(*this);
74 
75  // Resize *this so that the result can fit
76  this->resize (M2.m(), M3.n());
77 
78  // Call the multiply function in the base class
79  this->multiply(*this, M2, M3);
80 }
void resize(const unsigned int new_m, const unsigned int new_n)
Definition: dense_matrix.h:776
static void multiply(DenseMatrixBase< T > &M1, const DenseMatrixBase< T > &M2, const DenseMatrixBase< T > &M3)
template<typename T>
void libMesh::DenseMatrix< T >::left_multiply_transpose ( const DenseMatrix< T > &  A)

Left multiplies by the transpose of the matrix A.

Definition at line 85 of file dense_matrix.C.

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

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

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

Left multiplies by the transpose of the matrix A which contains a different numerical type.

Definition at line 151 of file dense_matrix.C.

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

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

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

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

Definition at line 1017 of file dense_matrix.h.

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

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

1018 {
1019  libmesh_assert (this->_m);
1020  libmesh_assert (this->_n);
1021 
1022  Real rowsum = 0.;
1023  for (unsigned int j=0; j!=this->_n; j++)
1024  {
1025  rowsum += std::abs((*this)(0,j));
1026  }
1027  Real my_max = rowsum;
1028  for (unsigned int i=1; i!=this->_m; i++)
1029  {
1030  rowsum = 0.;
1031  for (unsigned int j=0; j!=this->_n; j++)
1032  {
1033  rowsum += std::abs((*this)(i,j));
1034  }
1035  my_max = (my_max > rowsum? my_max : rowsum);
1036  }
1037  return my_max;
1038 }
double abs(double a)
libmesh_assert(j)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
template<typename T>
void libMesh::DenseMatrix< T >::lu_solve ( const DenseVector< T > &  b,
DenseVector< T > &  x 
)

Solve the system Ax=b given the input vector b. Partial pivoting is performed by default in order to keep the algorithm stable to the effects of round-off error.

Definition at line 589 of file dense_matrix.C.

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

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

Definition at line 97 of file dense_matrix_base.h.

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

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

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

Definition at line 969 of file dense_matrix.h.

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

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

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

Definition at line 948 of file dense_matrix.h.

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

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

949 {
950  libmesh_assert (this->_m);
951  libmesh_assert (this->_n);
952  Real my_min = libmesh_real((*this)(0,0));
953 
954  for (unsigned int i=0; i!=this->_m; i++)
955  {
956  for (unsigned int j=0; j!=this->_n; j++)
957  {
958  Real current = libmesh_real((*this)(i,j));
959  my_min = (my_min < current? my_min : current);
960  }
961  }
962  return my_min;
963 }
T libmesh_real(T a)
libmesh_assert(j)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
template<typename T >
void libMesh::DenseMatrixBase< T >::multiply ( DenseMatrixBase< T > &  M1,
const DenseMatrixBase< T > &  M2,
const DenseMatrixBase< T > &  M3 
)
staticprotectedinherited

Helper function - Performs the computation M1 = M2 * M3 where: M1 = (m x n) M2 = (m x p) M3 = (p x n)

Definition at line 31 of file dense_matrix_base.C.

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

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

Definition at line 102 of file dense_matrix_base.h.

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

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

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

Definition at line 911 of file dense_matrix.h.

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

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

912 {
913  for (std::size_t i=0; i<_val.size(); i++)
914  if (_val[i] != mat._val[i])
915  return true;
916 
917  return false;
918 }
std::vector< T > _val
Definition: dense_matrix.h:516
template<typename T >
T libMesh::DenseMatrix< T >::operator() ( const unsigned int  i,
const unsigned int  j 
) const
inline
Returns
The (i,j) element of the matrix.

Definition at line 818 of file dense_matrix.h.

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

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

820 {
821  libmesh_assert_less (i*j, _val.size());
822  libmesh_assert_less (i, this->_m);
823  libmesh_assert_less (j, this->_n);
824 
825 
826  // return _val[(i) + (this->_m)*(j)]; // col-major
827  return _val[(i)*(this->_n) + (j)]; // row-major
828 }
std::vector< T > _val
Definition: dense_matrix.h:516
template<typename T >
T & libMesh::DenseMatrix< T >::operator() ( const unsigned int  i,
const unsigned int  j 
)
inline
Returns
The (i,j) element of the matrix as a writable reference.

Definition at line 834 of file dense_matrix.h.

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

836 {
837  libmesh_assert_less (i*j, _val.size());
838  libmesh_assert_less (i, this->_m);
839  libmesh_assert_less (j, this->_n);
840 
841  //return _val[(i) + (this->_m)*(j)]; // col-major
842  return _val[(i)*(this->_n) + (j)]; // row-major
843 }
std::vector< T > _val
Definition: dense_matrix.h:516
template<typename T>
DenseMatrix< T > & libMesh::DenseMatrix< T >::operator*= ( const T  factor)
inline

Multiplies every element in the matrix by factor.

Returns
A reference to *this.

Definition at line 870 of file dense_matrix.h.

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

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

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

Adds mat to this matrix.

Returns
A reference to *this.

Definition at line 924 of file dense_matrix.h.

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

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

925 {
926  for (std::size_t i=0; i<_val.size(); i++)
927  _val[i] += mat._val[i];
928 
929  return *this;
930 }
std::vector< T > _val
Definition: dense_matrix.h:516
template<typename T>
DenseMatrix< T > & libMesh::DenseMatrix< T >::operator-= ( const DenseMatrix< T > &  mat)
inline

Subtracts mat from this matrix.

Returns
A reference to *this.

Definition at line 936 of file dense_matrix.h.

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

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

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

Assignment operator.

Returns
A reference to *this.

Definition at line 803 of file dense_matrix.h.

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

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

804 {
805  this->_m = other_matrix._m;
806  this->_n = other_matrix._n;
807 
808  _val = other_matrix._val;
809  _decomposition_type = other_matrix._decomposition_type;
810 
811  return *this;
812 }
DecompositionType _decomposition_type
Definition: dense_matrix.h:566
std::vector< T > _val
Definition: dense_matrix.h:516
template<typename T >
template<typename T2 >
DenseMatrix< T > & libMesh::DenseMatrix< T >::operator= ( const DenseMatrix< T2 > &  other_matrix)
inline

Assignment-from-other-matrix-type operator.

Copies the dense matrix of type T2 into the present matrix. This is useful for copying real matrices into complex ones for further operations.

Returns
A reference to *this.

Definition at line 761 of file dense_matrix.h.

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

762 {
763  unsigned int mat_m = mat.m(), mat_n = mat.n();
764  this->resize(mat_m, mat_n);
765  for (unsigned int i=0; i<mat_m; i++)
766  for (unsigned int j=0; j<mat_n; j++)
767  (*this)(i,j) = mat(i,j);
768 
769  return *this;
770 }
void resize(const unsigned int new_m, const unsigned int new_n)
Definition: dense_matrix.h:776
template<typename T>
bool libMesh::DenseMatrix< T >::operator== ( const DenseMatrix< T > &  mat) const
inline
Returns
true if mat is exactly equal to this matrix, false otherwise.

Definition at line 898 of file dense_matrix.h.

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

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

899 {
900  for (std::size_t i=0; i<_val.size(); i++)
901  if (_val[i] != mat._val[i])
902  return false;
903 
904  return true;
905 }
std::vector< T > _val
Definition: dense_matrix.h:516
template<typename T >
void libMesh::DenseMatrixBase< T >::print ( std::ostream &  os = libMesh::out) const
inherited

Pretty-print the matrix, by default to libMesh::out.

Definition at line 128 of file dense_matrix_base.C.

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

129 {
130  for (unsigned int i=0; i<this->m(); i++)
131  {
132  for (unsigned int j=0; j<this->n(); j++)
133  os << std::setw(8)
134  << this->el(i,j) << " ";
135 
136  os << std::endl;
137  }
138 
139  return;
140 }
unsigned int n() const
unsigned int m() const
virtual T el(const unsigned int i, const unsigned int j) const =0
template<typename T >
void libMesh::DenseMatrixBase< T >::print_scientific ( std::ostream &  os,
unsigned  precision = 8 
) const
inherited

Prints the matrix entries with more decimal places in scientific notation.

Definition at line 86 of file dense_matrix_base.C.

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

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

Definition at line 776 of file dense_matrix.h.

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

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

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

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

Implements libMesh::DenseMatrixBase< T >.

Definition at line 210 of file dense_matrix.C.

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

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

211 {
212  if (this->use_blas_lapack)
213  this->_multiply_blas(M3, RIGHT_MULTIPLY);
214  else
215  {
216  // (*this) <- M3 * (*this)
217  // Where:
218  // (*this) = (m x n),
219  // M2 = (m x p),
220  // M3 = (p x n)
221 
222  // M2 is a copy of *this before it gets resize()d
223  DenseMatrix<T> M2(*this);
224 
225  // Resize *this so that the result can fit
226  this->resize (M2.m(), M3.n());
227 
228  this->multiply(*this, M2, M3);
229  }
230 }
void _multiply_blas(const DenseMatrixBase< T > &other, _BLAS_Multiply_Flag flag)
void resize(const unsigned int new_m, const unsigned int new_n)
Definition: dense_matrix.h:776
static void multiply(DenseMatrixBase< T > &M1, const DenseMatrixBase< T > &M2, const DenseMatrixBase< T > &M3)
template<typename T >
template<typename T2 >
void libMesh::DenseMatrix< T >::right_multiply ( const DenseMatrixBase< T2 > &  M2)

Right multiplies by the matrix M2 of different type

Definition at line 236 of file dense_matrix.C.

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

237 {
238  // (*this) <- M3 * (*this)
239  // Where:
240  // (*this) = (m x n),
241  // M2 = (m x p),
242  // M3 = (p x n)
243 
244  // M2 is a copy of *this before it gets resize()d
245  DenseMatrix<T> M2(*this);
246 
247  // Resize *this so that the result can fit
248  this->resize (M2.m(), M3.n());
249 
250  this->multiply(*this, M2, M3);
251 }
void resize(const unsigned int new_m, const unsigned int new_n)
Definition: dense_matrix.h:776
static void multiply(DenseMatrixBase< T > &M1, const DenseMatrixBase< T > &M2, const DenseMatrixBase< T > &M3)
template<typename T>
void libMesh::DenseMatrix< T >::right_multiply_transpose ( const DenseMatrix< T > &  A)

Right multiplies by the transpose of the matrix A

Definition at line 257 of file dense_matrix.C.

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

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

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

Right multiplies by the transpose of the matrix A which contains a different numerical type.

Definition at line 322 of file dense_matrix.C.

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

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

Multiplies every element in the matrix by factor.

Definition at line 851 of file dense_matrix.h.

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

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

852 {
853  for (std::size_t i=0; i<_val.size(); i++)
854  _val[i] *= factor;
855 }
std::vector< T > _val
Definition: dense_matrix.h:516
template<typename T>
void libMesh::DenseMatrix< T >::scale_column ( const unsigned int  col,
const T  factor 
)
inline

Multiplies every element in the column col matrix by factor.

Definition at line 860 of file dense_matrix.h.

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

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

861 {
862  for (unsigned int i=0; i<this->m(); i++)
863  (*this)(i, col) *= factor;
864 }
unsigned int m() const
template<typename T >
void libMesh::DenseMatrix< T >::svd ( DenseVector< Real > &  sigma)

Compute the singular value decomposition of the matrix. On exit, sigma holds all of the singular values (in descending order).

The implementation uses PETSc's interface to BLAS/LAPACK. If this is not available, this function throws an error.

Definition at line 777 of file dense_matrix.C.

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

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

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

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

799 {
800  _svd_solve_lapack(rhs, x, rcond);
801 }
PetscErrorCode Vec x
void _svd_solve_lapack(const DenseVector< T > &rhs, DenseVector< T > &x, Real rcond) const
template<typename T>
void libMesh::DenseMatrix< T >::swap ( DenseMatrix< T > &  other_matrix)
inline

STL-like swap method

Definition at line 746 of file dense_matrix.h.

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

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

747 {
748  std::swap(this->_m, other_matrix._m);
749  std::swap(this->_n, other_matrix._n);
750  _val.swap(other_matrix._val);
752  _decomposition_type = other_matrix._decomposition_type;
753  other_matrix._decomposition_type = _temp;
754 }
DecompositionType _decomposition_type
Definition: dense_matrix.h:566
void swap(Iterator &lhs, Iterator &rhs)
std::vector< T > _val
Definition: dense_matrix.h:516
template<typename T >
T libMesh::DenseMatrix< T >::transpose ( const unsigned int  i,
const unsigned int  j 
) const
inline
Returns
The (i,j) element of the transposed matrix.

Definition at line 1044 of file dense_matrix.h.

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

1046 {
1047  // Implement in terms of operator()
1048  return (*this)(j,i);
1049 }
template<typename T>
void libMesh::DenseMatrix< T >::vector_mult ( DenseVector< T > &  dest,
const DenseVector< T > &  arg 
) const

Performs the matrix-vector multiplication, dest := (*this) * arg.

Definition at line 382 of file dense_matrix.C.

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

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

384 {
385  // Make sure the input sizes are compatible
386  libmesh_assert_equal_to (this->n(), arg.size());
387 
388  // Resize and clear dest.
389  // Note: DenseVector::resize() also zeros the vector.
390  dest.resize(this->m());
391 
392  // Short-circuit if the matrix is empty
393  if(this->m() == 0 || this->n() == 0)
394  return;
395 
396  if (this->use_blas_lapack)
397  this->_matvec_blas(1., 0., dest, arg);
398  else
399  {
400  const unsigned int n_rows = this->m();
401  const unsigned int n_cols = this->n();
402 
403  for (unsigned int i=0; i<n_rows; i++)
404  for (unsigned int j=0; j<n_cols; j++)
405  dest(i) += (*this)(i,j)*arg(j);
406  }
407 }
unsigned int n() const
void _matvec_blas(T alpha, T beta, DenseVector< T > &dest, const DenseVector< T > &arg, bool trans=false) const
unsigned int m() const
template<typename T>
template<typename T2 >
void libMesh::DenseMatrix< T >::vector_mult ( DenseVector< typename CompareTypes< T, T2 >::supertype > &  dest,
const DenseVector< T2 > &  arg 
) const

Performs the matrix-vector multiplication, dest := (*this) * arg on mixed types

Definition at line 413 of file dense_matrix.C.

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

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

Performs the scaled matrix-vector multiplication, dest += factor * (*this) * arg.

Definition at line 508 of file dense_matrix.C.

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

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

511 {
512  // Short-circuit if the matrix is empty
513  if (this->m() == 0)
514  {
515  dest.resize(0);
516  return;
517  }
518 
519  if (this->use_blas_lapack)
520  this->_matvec_blas(factor, 1., dest, arg);
521  else
522  {
523  DenseVector<T> temp(arg.size());
524  this->vector_mult(temp, arg);
525  dest.add(factor, temp);
526  }
527 }
void _matvec_blas(T alpha, T beta, DenseVector< T > &dest, const DenseVector< T > &arg, bool trans=false) const
unsigned int m() const
void vector_mult(DenseVector< T > &dest, const DenseVector< T > &arg) const
Definition: dense_matrix.C:382
template<typename T>
template<typename T2 , typename T3 >
void libMesh::DenseMatrix< T >::vector_mult_add ( DenseVector< typename CompareTypes< T, typename CompareTypes< T2, T3 >::supertype >::supertype > &  dest,
const T2  factor,
const DenseVector< T3 > &  arg 
) const

Performs the scaled matrix-vector multiplication, dest += factor * (*this) * arg. on mixed types

Definition at line 533 of file dense_matrix.C.

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

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

Performs the matrix-vector multiplication, dest := (*this)^T * arg.

Definition at line 438 of file dense_matrix.C.

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

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

440 {
441  // Make sure the input sizes are compatible
442  libmesh_assert_equal_to (this->m(), arg.size());
443 
444  // Resize and clear dest.
445  // Note: DenseVector::resize() also zeros the vector.
446  dest.resize(this->n());
447 
448  // Short-circuit if the matrix is empty
449  if (this->m() == 0)
450  return;
451 
452  if (this->use_blas_lapack)
453  {
454  this->_matvec_blas(1., 0., dest, arg, /*trans=*/true);
455  }
456  else
457  {
458  const unsigned int n_rows = this->m();
459  const unsigned int n_cols = this->n();
460 
461  // WORKS
462  // for (unsigned int j=0; j<n_cols; j++)
463  // for (unsigned int i=0; i<n_rows; i++)
464  // dest(j) += (*this)(i,j)*arg(i);
465 
466  // ALSO WORKS, (i,j) just swapped
467  for (unsigned int i=0; i<n_cols; i++)
468  for (unsigned int j=0; j<n_rows; j++)
469  dest(i) += (*this)(j,i)*arg(j);
470  }
471 }
unsigned int n() const
void _matvec_blas(T alpha, T beta, DenseVector< T > &dest, const DenseVector< T > &arg, bool trans=false) const
unsigned int m() const
template<typename T>
template<typename T2 >
void libMesh::DenseMatrix< T >::vector_mult_transpose ( DenseVector< typename CompareTypes< T, T2 >::supertype > &  dest,
const DenseVector< T2 > &  arg 
) const

Performs the matrix-vector multiplication, dest := (*this)^T * arg. on mixed types

Definition at line 477 of file dense_matrix.C.

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

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

Member Data Documentation

template<typename T>
DecompositionType libMesh::DenseMatrix< T >::_decomposition_type
private

This flag keeps track of which type of decomposition has been performed on the matrix.

Definition at line 566 of file dense_matrix.h.

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

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

Definition at line 656 of file dense_matrix.h.

template<typename T>
bool libMesh::DenseMatrix< T >::use_blas_lapack

Computes the inverse of the dense matrix (assuming it is invertible) by first computing the LU decomposition and then performing multiple back substitution steps. Follows the algorithm from Numerical Recipes in C that is available on the web.

This routine is commented out since it is not really a memory- or computationally- efficient implementation. Also, you typically don't need the actual inverse for anything, and can use something like lu_solve() instead. Run-time selectable option to turn on/off BLAS support. This was primarily used for testing purposes, and could be removed...

Definition at line 509 of file dense_matrix.h.


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