libMesh::DenseMatrix< T > Class Template Reference

#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 ()
T	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
virtual T &	el (const unsigned int i, const unsigned int j)
virtual void	left_multiply (const DenseMatrixBase< T > &M2)
virtual void	right_multiply (const DenseMatrixBase< T > &M3)
void	vector_mult (DenseVector< T > &dest, const DenseVector< T > &arg) const
void	vector_mult_transpose (DenseVector< T > &dest, const DenseVector< T > &arg) const
void	vector_mult_add (DenseVector< T > &dest, const T factor, const DenseVector< T > &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)
void	swap (DenseMatrix< T > &other_matrix)
void	resize (const unsigned int new_m, const unsigned int new_n)
void	scale (const T factor)
DenseMatrix< T > &	operator*= (const T factor)
void	add (const T factor, const DenseMatrix< T > &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)
void	right_multiply_transpose (const DenseMatrix< T > &A)
T	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< T > &sigma)
void	svd (DenseVector< T > &sigma, DenseMatrix< T > &U, DenseMatrix< T > &VT)
T	det ()
virtual void	left_multiply (const DenseMatrixBase< T > &M2)=0
virtual void	right_multiply (const DenseMatrixBase< T > &M3)=0
unsigned int	m () const
unsigned int	n () const
void	print (std::ostream &os=libMesh::out) const
void	print_scientific (std::ostream &os) 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	multiply (DenseMatrixBase< T > &M1, const DenseMatrixBase< T > &M2, const DenseMatrixBase< T > &M3)
void	condense (const unsigned int i, const unsigned int j, const T val, DenseVectorBase< T > &rhs)
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 }
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< T > &sigma)
void	_svd_lapack (DenseVector< T > &sigma, DenseMatrix< T > &U, DenseMatrix< T > &VT)
void	_svd_helper (char JOBU, char JOBVT, std::vector< T > &sigma_val, std::vector< T > &U_val, std::vector< T > &VT_val)
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< int >	_pivots
Friends
std::ostream &	operator<< (std::ostream &os, const DenseMatrixBase< T > &m)

Detailed Description

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

Defines a dense matrix for use in Finite Element-type computations. Useful for storing element stiffness matrices before summation into a global matrix.

Author:: Benjamin S. Kirk, 2002

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

00411                            {
00412     LEFT_MULTIPLY = 0,
00413     RIGHT_MULTIPLY,
00414     LEFT_MULTIPLY_TRANSPOSE,
00415     RIGHT_MULTIPLY_TRANSPOSE
00416   };

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

00399 {LU=0, CHOLESKY=1, LU_BLAS_LAPACK, NONE};

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

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

00546   : DenseMatrixBase<T>(new_m,new_n),
00547 #if defined(LIBMESH_HAVE_PETSC) && defined(LIBMESH_USE_REAL_NUMBERS) && defined(LIBMESH_DEFAULT_DOUBLE_PRECISION)
00548     use_blas_lapack(true),
00549 #else
00550     use_blas_lapack(false),
00551 #endif
00552     _val(),
00553     _decomposition_type(NONE),
00554     _pivots()
00555 {
00556   this->resize(new_m,new_n);
00557 }

template<typename T>

virtual libMesh::DenseMatrix< T >::~DenseMatrix ( ) [inline, virtual]

Copy-constructor. Destructor. Empty.

Definition at line 69 of file dense_matrix.h.

00069 {}

Member Function Documentation

template<typename T >

template<typename T2 >

void libMesh::DenseMatrix< T >::_cholesky_back_substitute	(	const DenseVector< T2 > &	b,
		DenseVector< T2 > &	x
	)			const `[inline, private]`

Solves the equation Ax=b for the unknown value x and rhs b based on the Cholesky factorization of A. Note that this method may be used when A is real-valued and b and x are complex-valued.

Definition at line 716 of file dense_matrix.C.

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

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

00718 {
00719   // Shorthand notation for number of rows and columns.
00720   const unsigned int
00721     n_rows = this->m(),
00722     n_cols = this->n();
00723 
00724   // Just to be really sure...
00725   libmesh_assert_equal_to (n_rows, n_cols);
00726 
00727   // A convenient reference to *this
00728   const DenseMatrix<T>& A = *this;
00729 
00730   // Now compute the solution to Ax =b using the factorization.
00731   x.resize(n_rows);
00732 
00733   // Solve for Ly=b
00734   for (unsigned int i=0; i<n_cols; ++i)
00735     {
00736       T2 temp = b(i);
00737 
00738       for (unsigned int k=0; k<i; ++k)
00739         temp -= A(i,k)*x(k);
00740 
00741       x(i) = temp / A(i,i);
00742     }
00743 
00744   // Solve for L^T x = y
00745   for (unsigned int i=0; i<n_cols; ++i)
00746     {
00747       const unsigned int ib = (n_cols-1)-i;
00748 
00749       for (unsigned int k=(ib+1); k<n_cols; ++k)
00750         x(ib) -= A(k,ib) * x(k);
00751 
00752       x(ib) /= A(ib,ib);
00753     }
00754 }

template<typename T >

void libMesh::DenseMatrix< T >::_cholesky_decompose ( ) [inline, private]

Decomposes a symmetric positive definite matrix into a product of two lower triangular matrices according to A = LL^T. Note that this program generates an error if the matrix is not SPD.

Definition at line 665 of file dense_matrix.C.

References libMesh::DenseMatrix< T >::_decomposition_type, libMesh::DenseMatrix< T >::CHOLESKY, libMesh::err, libMesh::DenseMatrixBase< T >::m(), libMesh::DenseMatrixBase< T >::n(), and libMesh::DenseMatrix< T >::NONE.

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

00666 {
00667   // If we called this function, there better not be any
00668   // previous decomposition of the matrix.
00669   libmesh_assert_equal_to (this->_decomposition_type, NONE);
00670 
00671   // Shorthand notation for number of rows and columns.
00672   const unsigned int
00673     n_rows = this->m(),
00674     n_cols = this->n();
00675 
00676   // Just to be really sure...
00677   libmesh_assert_equal_to (n_rows, n_cols);
00678 
00679   // A convenient reference to *this
00680   DenseMatrix<T>& A = *this;
00681 
00682   for (unsigned int i=0; i<n_rows; ++i)
00683     {
00684       for (unsigned int j=i; j<n_cols; ++j)
00685         {
00686           for (unsigned int k=0; k<i; ++k)
00687             A(i,j) -= A(i,k) * A(j,k);
00688 
00689           if (i == j)
00690             {
00691 #ifndef LIBMESH_USE_COMPLEX_NUMBERS
00692               if (A(i,j) <= 0.0)
00693                 {
00694                   libMesh::err << "Error! Can only use Cholesky decomposition "
00695                                 << "with symmetric positive definite matrices."
00696                                 << std::endl;
00697                   libmesh_error();
00698                 }
00699 #endif
00700 
00701               A(i,i) = std::sqrt(A(i,j));
00702             }
00703           else
00704             A(j,i) = A(i,j) / A(i,i);
00705         }
00706     }
00707 
00708   // Set the flag for CHOLESKY decomposition
00709   this->_decomposition_type = CHOLESKY;
00710 }

template<typename T>

void libMesh::DenseMatrix< T >::_lu_back_substitute	(	const DenseVector< T > &	b,
		DenseVector< T > &	x
	)			const `[inline, 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 413 of file dense_matrix.C.

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

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

00415 {
00416   const unsigned int
00417     n_cols = this->n();
00418 
00419   libmesh_assert_equal_to (this->m(), n_cols);
00420   libmesh_assert_equal_to (this->m(), b.size());
00421 
00422   x.resize (n_cols);
00423 
00424   // A convenient reference to *this
00425   const DenseMatrix<T>& A = *this;
00426 
00427   // Temporary vector storage.  We use this instead of
00428   // modifying the RHS.
00429   DenseVector<T> z = b;
00430 
00431   // Lower-triangular "top to bottom" solve step, taking into account pivots
00432   for (unsigned int i=0; i<n_cols; ++i)
00433     {
00434       // Swap
00435       if (_pivots[i] != static_cast<int>(i))
00436         std::swap( z(i), z(_pivots[i]) );
00437 
00438       x(i) = z(i);
00439 
00440       for (unsigned int j=0; j<i; ++j)
00441         x(i) -= A(i,j)*x(j);
00442 
00443       x(i) /= A(i,i);
00444     }
00445 
00446   // Upper-triangular "bottom to top" solve step
00447   const unsigned int last_row = n_cols-1;
00448 
00449   for (int i=last_row; i>=0; --i)
00450     {
00451       for (int j=i+1; j<static_cast<int>(n_cols); ++j)
00452         x(i) -= A(i,j)*x(j);
00453     }
00454 }

template<typename T>

void libMesh::DenseMatrix< T >::_lu_back_substitute_lapack	(	const DenseVector< T > &	b,
		DenseVector< T > &	x
	)			`[inline, 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 535 of file dense_matrix_blas_lapack.C.

References libMesh::DenseMatrix< T >::_pivots, libMesh::DenseMatrix< T >::_val, libMesh::DenseVector< T >::get_values(), libMesh::DenseMatrixBase< T >::m(), and libMesh::out.

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

00537 {
00538   // The calling sequence for getrs is:
00539   // dgetrs(TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
00540 
00541   //    trans   (input) char*
00542   //            'n' for no tranpose, 't' for transpose
00543   char TRANS[] = "t";
00544 
00545   //    N       (input) int*
00546   //            The order of the matrix A.  N >= 0.
00547   int N = this->m();
00548 
00549 
00550   //    NRHS    (input) int*
00551   //            The number of right hand sides, i.e., the number of columns
00552   //            of the matrix B.  NRHS >= 0.
00553   int NRHS = 1;
00554 
00555   //    A       (input) DOUBLE PRECISION array, dimension (LDA,N)
00556   //            The factors L and U from the factorization A = P*L*U
00557   //            as computed by dgetrf.
00558   // Here, we pass &(_val[0])
00559 
00560   //    LDA     (input) int*
00561   //            The leading dimension of the array A.  LDA >= max(1,N).
00562   int LDA = N;
00563 
00564   //    ipiv    (input) int array, dimension (N)
00565   //            The pivot indices from DGETRF; for 1<=i<=N, row i of the
00566   //            matrix was interchanged with row IPIV(i).
00567   // Here, we pass &(_pivots[0]) which was computed in _lu_decompose_lapack
00568 
00569   //    B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
00570   //            On entry, the right hand side matrix B.
00571   //            On exit, the solution matrix X.
00572   // Here, we pass a copy of the rhs vector's data array in x, so that the
00573   // passed right-hand side b is unmodified.  I don't see a way around this
00574   // copy if we want to maintain an unmodified rhs in LibMesh.
00575   x = b;
00576   std::vector<T>& x_vec = x.get_values();
00577 
00578   // We can avoid the copy if we don't care about overwriting the RHS: just
00579   // pass b to the Lapack routine and then swap with x before exiting
00580   // std::vector<T>& x_vec = b.get_values();
00581 
00582   //    LDB     (input) int*
00583   //            The leading dimension of the array B.  LDB >= max(1,N).
00584   int LDB = N;
00585 
00586   //    INFO    (output) int*
00587   //            = 0:  successful exit
00588   //            < 0:  if INFO = -i, the i-th argument had an illegal value
00589   int INFO = 0;
00590 
00591   // Finally, ready to call the Lapack getrs function
00592   LAPACKgetrs_(TRANS, &N, &NRHS, &(_val[0]), &LDA, &(_pivots[0]), &(x_vec[0]), &LDB, &INFO);
00593 
00594   // Check return value for errors
00595   if (INFO != 0)
00596     {
00597       libMesh::out << "INFO="
00598                     << INFO
00599                     << ", Error during Lapack LU solve!" << std::endl;
00600       libmesh_error();
00601     }
00602 
00603   // Don't do this if you already made a copy of b above
00604   // Swap b and x.  The solution will then be in x, and whatever was originally
00605   // in x, maybe garbage, maybe nothing, will be in b.
00606   // FIXME: Rewrite the LU and Cholesky solves to just take one input, and overwrite
00607   // the input.  This *should* make user code simpler, as they don't have to create
00608   // an extra vector just to pass it in to the solve function!
00609   // b.swap(x);
00610 }

template<typename T >

void libMesh::DenseMatrix< T >::_lu_decompose ( ) [inline, 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 464 of file dense_matrix.C.

References libMesh::DenseMatrix< T >::_decomposition_type, libMesh::DenseMatrix< T >::_pivots, std::abs(), libMesh::DenseMatrix< T >::LU, libMesh::DenseMatrixBase< T >::m(), libMesh::DenseMatrix< T >::NONE, libMesh::out, libMesh::Real, swap(), and libMesh::zero.

Referenced by libMesh::DenseMatrix< T >::det(), and libMesh::DenseMatrix< T >::lu_solve().

00466 00467 00468 00469 00470 00471 00472 00473 00474 00475 00476 00477 00478 00479 00480     { 00481 00482 00483 00484 00485 00486 00487         { 00488 00489 00490 00491 00492 00493 00494         } 00495 00496 00497 00498 00499 00500 00501 00502         { 00503 00504 00505         } 00506 00507 00508 00509 00510         { 00511 00512 00513         } 00514 00515 00516 00517 00518 00519 00520 00521 00522 00523 00524 00525 00526 00527 00528 00529 00530 00531 00532 00533 00534 00535 00536 00537 00538 00539 00540 }

class="fragment">00465 { // If this function was called, there better not be any // previous decomposition of the matrix. libmesh_assert_equal_to (this->_decomposition_type, NONE); // Get the matrix size and make sure it is square const unsigned int n_rows = this->m(); // A convenient reference to *this DenseMatrix<T>& A = *this; _pivots.resize(n_rows); for (unsigned int i=0; i<n_rows; ++i) // Find the pivot row by searching down the i'th column _pivots[i] = i; // std::abs(complex) must return a Real! Real the_max = std::abs( A(i,i) ); for (unsigned int j=i+1; j<n_rows; ++j) Real candidate_max = std::abs( A(j,i) ); if (the_max < candidate_max) { the_max = candidate_max; _pivots[i] = j; } // libMesh::out << "the_max=" << the_max << " found at row " << _pivots[i] << std::endl; // If the max was found in a different row, interchange rows. // Here we interchange the *entire* row, in Gaussian elimination // you would only interchange the subrows A(i,j) and A(p(i),j), for j>i if (_pivots[i] != static_cast<int>(i)) for (unsigned int j=0; j<n_rows; ++j) std::swap( A(i,j), A(_pivots[i], j) ); // If the max abs entry found is zero, the matrix is singular if (A(i,i) == libMesh::zero) libMesh::out << "Matrix A is singular!" << std::endl; libmesh_error(); // Scale upper triangle entries of row i by the diagonal entry // Note: don't scale the diagonal entry itself! const T diag_inv = 1. / A(i,i); for (unsigned int j=i+1; j<n_rows; ++j) A(i,j) *= diag_inv; // Update the remaining sub-matrix A[i+1:m][i+1:m] // by subtracting off (the diagonal-scaled) // upper-triangular part of row i, scaled by the // i'th column entry of each row. In terms of // row operations, this is: // for each r > i // SubRow(r) = SubRow(r) - A(r,i)*SubRow(i) // // If we were scaling the i'th column as well, like // in Gaussian elimination, this would 'zero' the // entry in the i'th column. for (unsigned int row=i+1; row<n_rows; ++row) for (unsigned int col=i+1; col<n_rows; ++col) A(row,col) -= A(row,i) * A(i,col); } // end i loop // Set the flag for LU decomposition this->_decomposition_type = LU;

template<typename T >

template void libMesh::DenseMatrix< T >::_lu_decompose_lapack ( ) [inline, 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 215 of file dense_matrix_blas_lapack.C.

References libMesh::DenseMatrix< T >::_decomposition_type, libMesh::DenseMatrix< T >::_pivots, libMesh::DenseMatrix< T >::_val, libMesh::DenseMatrix< T >::LU_BLAS_LAPACK, libMesh::DenseMatrixBase< T >::m(), std::min(), libMesh::DenseMatrixBase< T >::n(), libMesh::DenseMatrix< T >::NONE, and libMesh::out.

Referenced by libMesh::DenseMatrix< T >::det(), and libMesh::DenseMatrix< T >::lu_solve().

00216 {
00217   // If this function was called, there better not be any
00218   // previous decomposition of the matrix.
00219   libmesh_assert_equal_to (this->_decomposition_type, NONE);
00220 
00221   // The calling sequence for dgetrf is:
00222   // dgetrf(M, N, A, lda, ipiv, info)
00223 
00224   //    M       (input) int*
00225   //            The number of rows of the matrix A.  M >= 0.
00226   // In C/C++, pass the number of *cols* of A
00227   int M = this->n();
00228 
00229   //    N       (input) int*
00230   //            The number of columns of the matrix A.  N >= 0.
00231   // In C/C++, pass the number of *rows* of A
00232   int N = this->m();
00233 
00234   //    A (input/output) double precision array, dimension (LDA,N)
00235   //      On entry, the M-by-N matrix to be factored.
00236   //      On exit, the factors L and U from the factorization
00237   //      A = P*L*U; the unit diagonal elements of L are not stored.
00238   // Here, we pass &(_val[0]).
00239 
00240   //    LDA     (input) int*
00241   //            The leading dimension of the array A.  LDA >= max(1,M).
00242   int LDA = M;
00243 
00244   //    ipiv    (output) integer array, dimension (min(m,n))
00245   //            The pivot indices; for 1 <= i <= min(m,n), row i of the
00246   //            matrix was interchanged with row IPIV(i).
00247   // Here, we pass &(_pivots[0]), a private class member used to store pivots
00248   this->_pivots.resize( std::min(M,N) );
00249 
00250   //    info    (output) int*
00251   //            = 0:  successful exit
00252   //            < 0:  if INFO = -i, the i-th argument had an illegal value
00253   //            > 0:  if INFO = i, U(i,i) is exactly zero. The factorization
00254   //                  has been completed, but the factor U is exactly
00255   //                  singular, and division by zero will occur if it is used
00256   //                  to solve a system of equations.
00257   int INFO = 0;
00258 
00259   // Ready to call the actual factorization routine through PETSc's interface
00260   LAPACKgetrf_(&M, &N, &(this->_val[0]), &LDA, &(_pivots[0]), &INFO);
00261 
00262   // Check return value for errors
00263   if (INFO != 0)
00264     {
00265       libMesh::out << "INFO="
00266                     << INFO
00267                     << ", Error during Lapack LU factorization!" << std::endl;
00268       libmesh_error();
00269     }
00270 
00271   // Set the flag for LU decomposition
00272   this->_decomposition_type = LU_BLAS_LAPACK;
00273 }

template<typename T>

void libMesh::DenseMatrix< T >::_matvec_blas	(	T	alpha,
		T	beta,
		DenseVector< T > &	dest,
		const DenseVector< T > &	arg,
		bool	trans = `false`
	)			const `[inline, 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 631 of file dense_matrix_blas_lapack.C.

References libMesh::DenseVector< T >::get_values(), libMesh::DenseMatrix< T >::get_values(), libMesh::DenseMatrixBase< T >::m(), libMesh::DenseMatrixBase< T >::n(), libMesh::out, and libMesh::DenseVector< T >::size().

Referenced by libMesh::DenseMatrix< T >::vector_mult(), libMesh::DenseMatrix< T >::vector_mult_add(), and libMesh::DenseMatrix< T >::vector_mult_transpose().

00635 {
00636   // Ensure that dest and arg sizes are compatible
00637   if (!trans)
00638     {
00639       // dest  ~ A     * arg
00640       // (mx1)   (mxn) * (nx1)
00641       if ((dest.size() != this->m()) || (arg.size() != this->n()))
00642         {
00643           libMesh::out << "Improper input argument sizes!" << std::endl;
00644           libmesh_error();
00645         }
00646     }
00647 
00648   else // trans == true
00649     {
00650       // Ensure that dest and arg are proper size
00651       // dest  ~ A^T   * arg
00652       // (nx1)   (nxm) * (mx1)
00653       if ((dest.size() != this->n()) || (arg.size() != this->m()))
00654         {
00655           libMesh::out << "Improper input argument sizes!" << std::endl;
00656           libmesh_error();
00657         }
00658     }
00659 
00660   // Calling sequence for dgemv:
00661   //
00662   // dgemv(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
00663 
00664   //   TRANS  - CHARACTER*1, 't' for transpose, 'n' for non-transpose multiply
00665   // We store everything in row-major order, so pass the transpose flag for
00666   // non-transposed matvecs and the 'n' flag for transposed matvecs
00667   char TRANS[] = "t";
00668   if (trans)
00669     TRANS[0] = 'n';
00670 
00671   //   M      - INTEGER.
00672   //            On entry, M specifies the number of rows of the matrix A.
00673   // In C/C++, pass the number of *cols* of A
00674   int M = this->n();
00675 
00676   //   N      - INTEGER.
00677   //            On entry, N specifies the number of columns of the matrix A.
00678   // In C/C++, pass the number of *rows* of A
00679   int N = this->m();
00680 
00681   //   ALPHA  - DOUBLE PRECISION.
00682   // The scalar constant passed to this function
00683 
00684   //   A      - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
00685   //            Before entry, the leading m by n part of the array A must
00686   //            contain the matrix of coefficients.
00687   // The matrix, *this.  Note that _matvec_blas is called from
00688   // a const function, vector_mult(), and so we have made this function const
00689   // as well.  Since BLAS knows nothing about const, we have to cast it away
00690   // now.
00691   DenseMatrix<T>& a_ref = const_cast< DenseMatrix<T>& > ( *this );
00692   std::vector<T>& a = a_ref.get_values();
00693 
00694   //   LDA    - INTEGER.
00695   //            On entry, LDA specifies the first dimension of A as declared
00696   //            in the calling (sub) program. LDA must be at least
00697   //            max( 1, m ).
00698   int LDA = M;
00699 
00700   //   X      - DOUBLE PRECISION array of DIMENSION at least
00701   //            ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
00702   //            and at least
00703   //            ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
00704   //            Before entry, the incremented array X must contain the
00705   //            vector x.
00706   // Here, we must cast away the const-ness of "arg" since BLAS knows
00707   // nothing about const
00708   DenseVector<T>& x_ref = const_cast< DenseVector<T>& > ( arg );
00709   std::vector<T>& x = x_ref.get_values();
00710 
00711   //   INCX   - INTEGER.
00712   //            On entry, INCX specifies the increment for the elements of
00713   //            X. INCX must not be zero.
00714   int INCX = 1;
00715 
00716   //   BETA   - DOUBLE PRECISION.
00717   //            On entry, BETA specifies the scalar beta. When BETA is
00718   //            supplied as zero then Y need not be set on input.
00719   // The second scalar constant passed to this function
00720 
00721   //   Y      - DOUBLE PRECISION array of DIMENSION at least
00722   //            ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
00723   //            and at least
00724   //            ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
00725   //            Before entry with BETA non-zero, the incremented array Y
00726   //            must contain the vector y. On exit, Y is overwritten by the
00727   //            updated vector y.
00728   // The input vector "dest"
00729   std::vector<T>& y = dest.get_values();
00730 
00731   //   INCY   - INTEGER.
00732   //            On entry, INCY specifies the increment for the elements of
00733   //            Y. INCY must not be zero.
00734   int INCY = 1;
00735 
00736   // Finally, ready to call the BLAS function
00737   BLASgemv_(TRANS, &M, &N, &alpha, &(a[0]), &LDA, &(x[0]), &INCX, &beta, &(y[0]), &INCY);
00738 }

template<typename T>

void libMesh::DenseMatrix< T >::_multiply_blas	(	const DenseMatrixBase< T > &	other,
		_BLAS_Multiply_Flag	flag
	)			`[inline, private]`

The _multiply_blas function computes A <- op(A) * op(B) using BLAS gemm function. Used in the right_multiply(), left_multiply(), right_multiply_transpose(), and left_multiply_tranpose() routines. [ Implementation in dense_matrix_blas_lapack.C ]

Definition at line 40 of file dense_matrix_blas_lapack.C.

References libMesh::DenseMatrixBase< T >::_m, libMesh::DenseMatrixBase< T >::_n, libMesh::DenseMatrix< T >::_val, libMesh::DenseMatrix< T >::LEFT_MULTIPLY, libMesh::DenseMatrix< T >::LEFT_MULTIPLY_TRANSPOSE, libMesh::DenseMatrixBase< T >::m(), libMesh::DenseMatrixBase< T >::n(), libMesh::out, libMesh::DenseMatrix< T >::RIGHT_MULTIPLY, libMesh::DenseMatrix< T >::RIGHT_MULTIPLY_TRANSPOSE, and libMesh::DenseMatrix< T >::swap().

Referenced by libMesh::DenseMatrix< T >::left_multiply(), libMesh::DenseMatrix< T >::left_multiply_transpose(), libMesh::DenseMatrix< T >::right_multiply(), and libMesh::DenseMatrix< T >::right_multiply_transpose().

00042 {
00043   int result_size = 0;
00044 
00045   // For each case, determine the size of the final result make sure
00046   // that the inner dimensions match
00047   switch (flag)
00048     {
00049     case LEFT_MULTIPLY:
00050       {
00051         result_size = other.m() * this->n();
00052         if (other.n() == this->m())
00053           break;
00054       }
00055     case RIGHT_MULTIPLY:
00056       {
00057         result_size = other.n() * this->m();
00058         if (other.m() == this->n())
00059           break;
00060       }
00061     case LEFT_MULTIPLY_TRANSPOSE:
00062       {
00063         result_size = other.n() * this->n();
00064         if (other.m() == this->m())
00065           break;
00066       }
00067     case RIGHT_MULTIPLY_TRANSPOSE:
00068       {
00069         result_size = other.m() * this->m();
00070         if (other.n() == this->n())
00071           break;
00072       }
00073     default:
00074       {
00075         libMesh::out << "Unknown flag selected or matrices are ";
00076         libMesh::out << "incompatible for multiplication." << std::endl;
00077         libmesh_error();
00078       }
00079     }
00080 
00081   // For this to work, the passed arg. must actually be a DenseMatrix<T>
00082   const DenseMatrix<T>* const_that = libmesh_cast_ptr< const DenseMatrix<T>* >(&other);
00083 
00084   // Also, although 'that' is logically const in this BLAS routine,
00085   // the PETSc BLAS interface does not specify that any of the inputs are
00086   // const.  To use it, I must cast away const-ness.
00087   DenseMatrix<T>* that = const_cast< DenseMatrix<T>* > (const_that);
00088 
00089   // Initialize A, B pointers for LEFT_MULTIPLY* cases
00090   DenseMatrix<T>
00091     *A = this,
00092     *B = that;
00093 
00094   // For RIGHT_MULTIPLY* cases, swap the meaning of A and B.
00095   // Here is a full table of combinations we can pass to BLASgemm, and what the answer is when finished:
00096   // pass A B   -> (Fortran) -> A^T B^T -> (C++) -> (A^T B^T)^T -> (identity) -> B A   "lt multiply"
00097   // pass B A   -> (Fortran) -> B^T A^T -> (C++) -> (B^T A^T)^T -> (identity) -> A B   "rt multiply"
00098   // pass A B^T -> (Fortran) -> A^T B   -> (C++) -> (A^T B)^T   -> (identity) -> B^T A "lt multiply t"
00099   // pass B^T A -> (Fortran) -> B A^T   -> (C++) -> (B A^T)^T   -> (identity) -> A B^T "rt multiply t"
00100   if (flag==RIGHT_MULTIPLY || flag==RIGHT_MULTIPLY_TRANSPOSE)
00101     std::swap(A,B);
00102 
00103   // transa, transb values to pass to blas
00104   char
00105     transa[] = "n",
00106     transb[] = "n";
00107 
00108   // Integer values to pass to BLAS:
00109   //
00110   // M
00111   // In Fortran, the number of rows of op(A),
00112   // In the BLAS documentation, typically known as 'M'.
00113   //
00114   // In C/C++, we set:
00115   // M = n_cols(A) if (transa='n')
00116   //     n_rows(A) if (transa='t')
00117   int M = static_cast<int>( A->n() );
00118 
00119   // N
00120   // In Fortran, the number of cols of op(B), and also the number of cols of C.
00121   // In the BLAS documentation, typically known as 'N'.
00122   //
00123   // In C/C++, we set:
00124   // N = n_rows(B) if (transb='n')
00125   //     n_cols(B) if (transb='t')
00126   int N = static_cast<int>( B->m() );
00127 
00128   // K
00129   // In Fortran, the number of cols of op(A), and also
00130   // the number of rows of op(B). In the BLAS documentation,
00131   // typically known as 'K'.
00132   //
00133   // In C/C++, we set:
00134   // K = n_rows(A) if (transa='n')
00135   //     n_cols(A) if (transa='t')
00136   int K = static_cast<int>( A->m() );
00137 
00138   // LDA (leading dimension of A). In our cases,
00139   // LDA is always the number of columns of A.
00140   int LDA = static_cast<int>( A->n() );
00141 
00142   // LDB (leading dimension of B).  In our cases,
00143   // LDB is always the number of columns of B.
00144   int LDB = static_cast<int>( B->n() );
00145 
00146   if (flag == LEFT_MULTIPLY_TRANSPOSE)
00147     {
00148       transb[0] = 't';
00149       N = static_cast<int>( B->n() );
00150     }
00151 
00152   else if (flag == RIGHT_MULTIPLY_TRANSPOSE)
00153     {
00154       transa[0] = 't';
00155       std::swap(M,K);
00156     }
00157 
00158   // LDC (leading dimension of C).  LDC is the
00159   // number of columns in the solution matrix.
00160   int LDC = M;
00161 
00162   // Scalar values to pass to BLAS
00163   //
00164   // scalar multiplying the whole product AB
00165   T alpha = 1.;
00166 
00167   // scalar multiplying C, which is the original matrix.
00168   T beta  = 0.;
00169 
00170   // Storage for the result
00171   std::vector<T> result (result_size);
00172 
00173   // Finally ready to call the BLAS
00174   BLASgemm_(transa, transb, &M, &N, &K, &alpha, &(A->_val[0]), &LDA, &(B->_val[0]), &LDB, &beta, &result[0], &LDC);
00175 
00176   // Update the relevant dimension for this matrix.
00177   switch (flag)
00178     {
00179     case LEFT_MULTIPLY:            { this->_m = other.m(); break; }
00180     case RIGHT_MULTIPLY:           { this->_n = other.n(); break; }
00181     case LEFT_MULTIPLY_TRANSPOSE:  { this->_m = other.n(); break; }
00182     case RIGHT_MULTIPLY_TRANSPOSE: { this->_n = other.m(); break; }
00183     default:
00184       {
00185         libMesh::out << "Unknown flag selected." << std::endl;
00186         libmesh_error();
00187       }
00188     }
00189 
00190   // Swap my data vector with the result
00191   this->_val.swap(result);
00192 }

template<typename T>

void libMesh::DenseMatrix< T >::_svd_helper	(	char	JOBU,
		char	JOBVT,
		std::vector< T > &	sigma_val,
		std::vector< T > &	U_val,
		std::vector< T > &	VT_val
	)			`[inline, private]`

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

Definition at line 400 of file dense_matrix_blas_lapack.C.

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

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

00405 {
00406 
00407   //    M       (input) int*
00408   //            The number of rows of the matrix A.  M >= 0.
00409   // In C/C++, pass the number of *cols* of A
00410   int M = this->n();
00411 
00412   //    N       (input) int*
00413   //            The number of columns of the matrix A.  N >= 0.
00414   // In C/C++, pass the number of *rows* of A
00415   int N = this->m();
00416 
00417   int min_MN = (M < N) ? M : N;
00418   int max_MN = (M > N) ? M : N;
00419 
00420   //  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
00421   //          On entry, the M-by-N matrix A.
00422   //          On exit,
00423   //          if JOBU = 'O',  A is overwritten with the first min(m,n)
00424   //                          columns of U (the left singular vectors,
00425   //                          stored columnwise);
00426   //          if JOBVT = 'O', A is overwritten with the first min(m,n)
00427   //                          rows of V**T (the right singular vectors,
00428   //                          stored rowwise);
00429   //          if JOBU .ne. 'O' and JOBVT .ne. 'O', the contents of A
00430   //                          are destroyed.
00431   // Here, we pass &(_val[0]).
00432 
00433   //    LDA     (input) int*
00434   //            The leading dimension of the array A.  LDA >= max(1,M).
00435   int LDA = M;
00436 
00437   //  S       (output) DOUBLE PRECISION array, dimension (min(M,N))
00438   //          The singular values of A, sorted so that S(i) >= S(i+1).
00439   sigma_val.resize( min_MN );
00440 
00441   //  LDU     (input) INTEGER
00442   //          The leading dimension of the array U.  LDU >= 1; if
00443   //          JOBU = 'S' or 'A', LDU >= M.
00444   int LDU = M;
00445 
00446   //  U       (output) DOUBLE PRECISION array, dimension (LDU,UCOL)
00447   //          (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'.
00448   //          If JOBU = 'A', U contains the M-by-M orthogonal matrix U;
00449   //          if JOBU = 'S', U contains the first min(m,n) columns of U
00450   //          (the left singular vectors, stored columnwise);
00451   //          if JOBU = 'N' or 'O', U is not referenced.
00452   U_val.resize( LDU*M );
00453 
00454   //  LDVT    (input) INTEGER
00455   //          The leading dimension of the array VT.  LDVT >= 1; if
00456   //          JOBVT = 'A', LDVT >= N; if JOBVT = 'S', LDVT >= min(M,N).
00457   int LDVT = N;
00458 
00459   //  VT      (output) DOUBLE PRECISION array, dimension (LDVT,N)
00460   //          If JOBVT = 'A', VT contains the N-by-N orthogonal matrix
00461   //          V**T;
00462   //          if JOBVT = 'S', VT contains the first min(m,n) rows of
00463   //          V**T (the right singular vectors, stored rowwise);
00464   //          if JOBVT = 'N' or 'O', VT is not referenced.
00465   VT_val.resize( LDVT*N );
00466 
00467   //  LWORK   (input) INTEGER
00468   //          The dimension of the array WORK.
00469   //          LWORK >= MAX(1,3*MIN(M,N)+MAX(M,N),5*MIN(M,N)).
00470   //          For good performance, LWORK should generally be larger.
00471   //
00472   //          If LWORK = -1, then a workspace query is assumed; the routine
00473   //          only calculates the optimal size of the WORK array, returns
00474   //          this value as the first entry of the WORK array, and no error
00475   //          message related to LWORK is issued by XERBLA.
00476   int larger = (3*min_MN+max_MN > 5*min_MN) ? 3*min_MN+max_MN : 5*min_MN;
00477   int LWORK  = (larger > 1) ? larger : 1;
00478 
00479 
00480   //  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
00481   //          On exit, if INFO = 0, WORK(1) returns the optimal LWORK;
00482   //          if INFO > 0, WORK(2:MIN(M,N)) contains the unconverged
00483   //          superdiagonal elements of an upper bidiagonal matrix B
00484   //          whose diagonal is in S (not necessarily sorted). B
00485   //          satisfies A = U * B * VT, so it has the same singular values
00486   //          as A, and singular vectors related by U and VT.
00487   std::vector<T> WORK( LWORK );
00488 
00489   //  INFO    (output) INTEGER
00490   //          = 0:  successful exit.
00491   //          < 0:  if INFO = -i, the i-th argument had an illegal value.
00492   //          > 0:  if DBDSQR did not converge, INFO specifies how many
00493   //                superdiagonals of an intermediate bidiagonal form B
00494   //                did not converge to zero. See the description of WORK
00495   //                above for details.
00496   int INFO = 0;
00497 
00498   // Ready to call the actual factorization routine through PETSc's interface
00499   LAPACKgesvd_(&JOBU, &JOBVT, &M, &N, &(_val[0]), &LDA, &(sigma_val[0]), &(U_val[0]),
00500                &LDU, &(VT_val[0]), &LDVT, &(WORK[0]), &LWORK, &INFO);
00501 
00502   // Check return value for errors
00503   if (INFO != 0)
00504     {
00505       libMesh::out << "INFO="
00506                     << INFO
00507                     << ", Error during Lapack SVD calculation!" << std::endl;
00508       libmesh_error();
00509     }
00510 }

template<typename T>

void libMesh::DenseMatrix< T >::_svd_lapack	(	DenseVector< T > &	sigma,
		DenseMatrix< T > &	U,
		DenseMatrix< T > &	VT
	)			`[inline, private]`

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

Definition at line 334 of file dense_matrix_blas_lapack.C.

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

00335 {
00336   // The calling sequence for dgetrf is:
00337   // DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, INFO )
00338 
00339 
00340   //  JOBU    (input) CHARACTER*1
00341   //          Specifies options for computing all or part of the matrix U:
00342   //          = 'A':  all M columns of U are returned in array U:
00343   //          = 'S':  the first min(m,n) columns of U (the left singular
00344   //                  vectors) are returned in the array U;
00345   //          = 'O':  the first min(m,n) columns of U (the left singular
00346   //                  vectors) are overwritten on the array A;
00347   //          = 'N':  no columns of U (no left singular vectors) are
00348   //                  computed.
00349   char JOBU = 'S';
00350 
00351   //  JOBVT   (input) CHARACTER*1
00352   //          Specifies options for computing all or part of the matrix
00353   //          V**T:
00354   //          = 'A':  all N rows of V**T are returned in the array VT;
00355   //          = 'S':  the first min(m,n) rows of V**T (the right singular
00356   //                  vectors) are returned in the array VT;
00357   //          = 'O':  the first min(m,n) rows of V**T (the right singular
00358   //                  vectors) are overwritten on the array A;
00359   //          = 'N':  no rows of V**T (no right singular vectors) are
00360   //                  computed.
00361   char JOBVT = 'S';
00362 
00363   std::vector<T> sigma_val;
00364   std::vector<T> U_val;
00365   std::vector<T> VT_val;
00366 
00367   _svd_helper(JOBU, JOBVT, sigma_val, U_val, VT_val);
00368 
00369   // Load the singular values into sigma, ignore U_val and VT_val
00370   const unsigned int n_sigma_vals =
00371     libmesh_cast_int<unsigned int>(sigma_val.size());
00372   sigma.resize(n_sigma_vals);
00373   for(unsigned int i=0; i<n_sigma_vals; i++)
00374     sigma(i) = sigma_val[i];
00375 
00376   int M = this->n();
00377   int N = this->m();
00378   int min_MN = (M < N) ? M : N;
00379   U.resize(M,min_MN);
00380   for(unsigned int i=0; i<U.m(); i++)
00381     for(unsigned int j=0; j<U.n(); j++)
00382     {
00383       unsigned int index = i + j*U.n();  // Column major storage
00384       U(i,j) = U_val[index];
00385     }
00386 
00387   VT.resize(min_MN,N);
00388   for(unsigned int i=0; i<VT.m(); i++)
00389     for(unsigned int j=0; j<VT.n(); j++)
00390     {
00391       unsigned int index = i + j*U.n(); // Column major storage
00392       VT(i,j) = VT_val[index];
00393     }
00394 
00395 }

template<typename T>

void libMesh::DenseMatrix< T >::_svd_lapack ( DenseVector< T > & sigma ) [inline, private]

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

Definition at line 289 of file dense_matrix_blas_lapack.C.

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

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

00290 {
00291   // The calling sequence for dgetrf is:
00292   // DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, INFO )
00293 
00294 
00295   //  JOBU    (input) CHARACTER*1
00296   //          Specifies options for computing all or part of the matrix U:
00297   //          = 'A':  all M columns of U are returned in array U:
00298   //          = 'S':  the first min(m,n) columns of U (the left singular
00299   //                  vectors) are returned in the array U;
00300   //          = 'O':  the first min(m,n) columns of U (the left singular
00301   //                  vectors) are overwritten on the array A;
00302   //          = 'N':  no columns of U (no left singular vectors) are
00303   //                  computed.
00304   char JOBU = 'N';
00305 
00306   //  JOBVT   (input) CHARACTER*1
00307   //          Specifies options for computing all or part of the matrix
00308   //          V**T:
00309   //          = 'A':  all N rows of V**T are returned in the array VT;
00310   //          = 'S':  the first min(m,n) rows of V**T (the right singular
00311   //                  vectors) are returned in the array VT;
00312   //          = 'O':  the first min(m,n) rows of V**T (the right singular
00313   //                  vectors) are overwritten on the array A;
00314   //          = 'N':  no rows of V**T (no right singular vectors) are
00315   //                  computed.
00316   char JOBVT = 'N';
00317 
00318   std::vector<T> sigma_val;
00319   std::vector<T> U_val;
00320   std::vector<T> VT_val;
00321 
00322   _svd_helper(JOBU, JOBVT, sigma_val, U_val, VT_val);
00323 
00324   // Load the singular values into sigma, ignore U_val and VT_val
00325   const unsigned int n_sigma_vals =
00326     libmesh_cast_int<unsigned int>(sigma_val.size());
00327   sigma.resize(n_sigma_vals);
00328   for(unsigned int i=0; i<n_sigma_vals; i++)
00329     sigma(i) = sigma_val[i];
00330 
00331 }

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
	)			`[inline, inherited]`

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

Definition at line 184 of file dense_matrix_base.h.

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

00186 {
00187   libmesh_assert_equal_to (this->m(), mat.m());
00188   libmesh_assert_equal_to (this->n(), mat.n());
00189 
00190   for (unsigned int j=0; j<this->n(); j++)
00191     for (unsigned int i=0; i<this->m(); i++)
00192       this->el(i,j) += factor*mat.el(i,j);
00193 }

template<typename T>

void libMesh::DenseMatrix< T >::add	(	const T	factor,
		const DenseMatrix< T > &	mat
	)			`[inline]`

Adds factor times mat to this matrix.

Definition at line 691 of file dense_matrix.h.

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

00692 {
00693   for (unsigned int i=0; i<_val.size(); i++)
00694     _val[i] += factor * mat._val[i];
00695 }

template<typename T >

template<typename T2 >

void libMesh::DenseMatrix< T >::cholesky_solve	(	const DenseVector< T2 > &	b,
		DenseVector< T2 > &	x
	)			`[inline]`

For symmetric positive definite (SPD) matrices. A Cholesky factorization of A such that A = L L^T is about twice as fast as a standard LU factorization. Therefore you can use this method if you know a-priori that the matrix is SPD. If the matrix is not SPD, an error is generated. One nice property of cholesky decompositions is that they do not require pivoting for stability. Note that this method may also be used when A is real-valued and x and b are complex-valued.

Definition at line 628 of file dense_matrix.C.

References libMesh::DenseMatrix< T >::_cholesky_back_substitute(), libMesh::DenseMatrix< T >::_cholesky_decompose(), libMesh::DenseMatrix< T >::_decomposition_type, libMesh::DenseMatrix< T >::CHOLESKY, libMesh::err, and libMesh::DenseMatrix< T >::NONE.

Referenced by libMesh::FEGenericBase< OutputType >::coarsened_dof_values(), libMesh::FEGenericBase< OutputType >::compute_periodic_constraints(), libMesh::FEGenericBase< OutputType >::compute_proj_constraints(), libMesh::BoundaryProjectSolution::operator()(), libMesh::ProjectFEMSolution::operator()(), libMesh::ProjectSolution::operator()(), and libMesh::HPCoarsenTest::select_refinement().

00630 {
00631   // Check for a previous decomposition
00632   switch(this->_decomposition_type)
00633     {
00634     case NONE:
00635       {
00636         this->_cholesky_decompose ();
00637         break;
00638       }
00639 
00640     case CHOLESKY:
00641       {
00642         // Already factored, just need to call back_substitute.
00643         break;
00644       }
00645 
00646     default:
00647       {
00648         libMesh::err << "Error! This matrix already has a "
00649                       << "different decomposition..."
00650                       << std::endl;
00651         libmesh_error();
00652       }
00653     }
00654 
00655   // Perform back substitution
00656   this->_cholesky_back_substitute (b, x);
00657 }

template<typename T>

void libMesh::DenseMatrixBase< T >::condense	(	const unsigned int	i,
		const unsigned int	j,
		const T	val,
		DenseVectorBase< T > &	rhs
	)			`[inline, protected, inherited]`

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::DenseMatrixBase< T >::_m, libMesh::DenseMatrixBase< T >::el(), libMesh::DenseVectorBase< T >::el(), libMesh::DenseMatrixBase< T >::m(), libMesh::DenseMatrixBase< T >::n(), and libMesh::DenseVectorBase< T >::size().

00062 {
00063   libmesh_assert_equal_to (this->_m, rhs.size());
00064   libmesh_assert_equal_to (iv, jv);
00065 
00066 
00067   // move the known value into the RHS
00068   // and zero the column
00069   for (unsigned int i=0; i<this->m(); i++)
00070     {
00071       rhs.el(i) -= this->el(i,jv)*val;
00072       this->el(i,jv) = 0.;
00073     }
00074 
00075   // zero the row
00076   for (unsigned int j=0; j<this->n(); j++)
00077     this->el(iv,j) = 0.;
00078 
00079   this->el(iv,jv) = 1.;
00080   rhs.el(iv) = val;
00081 
00082 }

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

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

00277   { DenseMatrixBase<T>::condense (i, j, val, rhs); }

template<typename T >

T libMesh::DenseMatrix< T >::det ( ) [inline]

Returns:: the determinant of the matrix. Note that this means doing an LU decomposition and then computing the product of the diagonal terms. Therefore this is a non-const method.

Definition at line 562 of file dense_matrix.C.

References libMesh::DenseMatrix< T >::_decomposition_type, libMesh::DenseMatrix< T >::_lu_decompose(), libMesh::DenseMatrix< T >::_lu_decompose_lapack(), libMesh::DenseMatrix< T >::_pivots, libMesh::err, libMesh::DenseMatrix< T >::LU, libMesh::DenseMatrix< T >::LU_BLAS_LAPACK, libMesh::DenseMatrixBase< T >::m(), libMesh::DenseMatrix< T >::NONE, libMesh::Real, and libMesh::DenseMatrix< T >::use_blas_lapack.

00563 {
00564   switch(this->_decomposition_type)
00565     {
00566     case NONE:
00567       {
00568         // First LU decompose the matrix.
00569         // Note that the lu_decompose routine will check to see if the
00570         // matrix is square so we don't worry about it.
00571         if (this->use_blas_lapack)
00572           this->_lu_decompose_lapack();
00573         else
00574           this->_lu_decompose ();
00575       }
00576     case LU:
00577     case LU_BLAS_LAPACK:
00578       {
00579         // Already decomposed, don't do anything
00580         break;
00581       }
00582     default:
00583       {
00584       libMesh::err << "Error! Can't compute the determinant under "
00585                     << "the current decomposition."
00586                     << std::endl;
00587       libmesh_error();
00588       }
00589     }
00590 
00591   // A variable to keep track of the running product of diagonal terms.
00592   T determinant = 1.;
00593 
00594   // Loop over diagonal terms, computing the product.  In practice,
00595   // be careful because this value could easily become too large to
00596   // fit in a double or float.  To be safe, one should keep track of
00597   // the power (of 10) of the determinant in a separate variable
00598   // and maintain an order 1 value for the determinant itself.
00599   unsigned int n_interchanges = 0;
00600   for (unsigned int i=0; i<this->m(); i++)
00601     {
00602       if (this->_decomposition_type==LU)
00603         if (_pivots[i] != static_cast<int>(i))
00604           n_interchanges++;
00605 
00606       // Lapack pivots are 1-based!
00607       if (this->_decomposition_type==LU_BLAS_LAPACK)
00608         if (_pivots[i] != static_cast<int>(i+1))
00609           n_interchanges++;
00610 
00611       determinant *= (*this)(i,i);
00612     }
00613 
00614   // Compute sign of determinant, depends on number of row interchanges!
00615   // The sign should be (-1)^{n}, where n is the number of interchanges.
00616   Real sign = n_interchanges % 2 == 0 ? 1. : -1.;
00617 
00618   return sign*determinant;
00619 }

template<typename T>

virtual T& libMesh::DenseMatrix< T >::el	(	const unsigned int	i,
		const unsigned int	j
	)			`[inline, virtual]`

Returns:: the (i,j) element of the matrix as a writeable reference.

Implements libMesh::DenseMatrixBase< T >.

Definition at line 98 of file dense_matrix.h.

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

template<typename T>

virtual T libMesh::DenseMatrix< T >::el	(	const unsigned int	i,
		const unsigned int	j
	)			const `[inline, virtual]`

Returns:: the (i,j) element of the matrix as a writeable reference.

Implements libMesh::DenseMatrixBase< T >.

Definition at line 92 of file dense_matrix.h.

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

template<typename T>

void libMesh::DenseMatrix< T >::get_principal_submatrix	(	unsigned int	sub_m,
		DenseMatrix< T > &	dest
	)			const `[inline]`

Put the sub_m x sub_m principal submatrix into dest.

Definition at line 330 of file dense_matrix.C.

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

00331 {
00332   get_principal_submatrix(sub_m, sub_m, dest);
00333 }

template<typename T>

void libMesh::DenseMatrix< T >::get_principal_submatrix	(	unsigned int	sub_m,
		unsigned int	sub_n,
		DenseMatrix< T > &	dest
	)			const `[inline]`

Put the sub_m x sub_n principal submatrix into dest.

Definition at line 315 of file dense_matrix.C.

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

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

00318 {
00319   libmesh_assert( (sub_m <= this->m()) && (sub_n <= this->n()) );
00320 
00321   dest.resize(sub_m, sub_n);
00322   for(unsigned int i=0; i<sub_m; i++)
00323     for(unsigned int j=0; j<sub_n; j++)
00324       dest(i,j) = (*this)(i,j);
00325 }

template<typename T>

void libMesh::DenseMatrix< T >::get_transpose ( DenseMatrix< T > & dest ) const [inline]

Put the tranposed matrix into dest.

Definition at line 338 of file dense_matrix.C.

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

00339 {
00340   dest.resize(this->n(), this->m());
00341 
00342   for (unsigned int i=0; i<dest.m(); i++)
00343     for (unsigned int j=0; j<dest.n(); j++)
00344       dest(i,j) = (*this)(j,i);
00345 }

template<typename T>

const std::vector<T>& libMesh::DenseMatrix< T >::get_values ( ) const [inline]

Return a constant reference to the matrix values.

Definition at line 265 of file dense_matrix.h.

00265 { return _val; }

template<typename T>

std::vector<T>& libMesh::DenseMatrix< T >::get_values ( ) [inline]

Access to the values array. This should be used with caution but can be used to speed up code compilation significantly.

Definition at line 260 of file dense_matrix.h.

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

00260 { return _val; }

template<typename T >

Real libMesh::DenseMatrix< T >::l1_norm ( ) const [inline]

Return the l1-norm of the matrix, that is $|M|_1=max_{all columns j}\sum_{all rows i} |M_ij|$ , (max. sum of columns). This is the natural matrix norm that is compatible to the l1-norm for vectors, i.e. $|Mv|_1\leq |M|_1 |v|_1$ .

Definition at line 793 of file dense_matrix.h.

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

00794 {
00795   libmesh_assert (this->_m);
00796   libmesh_assert (this->_n);
00797 
00798   Real columnsum = 0.;
00799   for (unsigned int i=0; i!=this->_m; i++)
00800     {
00801       columnsum += std::abs((*this)(i,0));
00802     }
00803   Real my_max = columnsum;
00804   for (unsigned int j=1; j!=this->_n; j++)
00805     {
00806       columnsum = 0.;
00807       for (unsigned int i=0; i!=this->_m; i++)
00808         {
00809           columnsum += std::abs((*this)(i,j));
00810         }
00811       my_max = (my_max > columnsum? my_max : columnsum);
00812     }
00813   return my_max;
00814 }

template<typename T>

virtual void libMesh::DenseMatrixBase< T >::left_multiply ( const DenseMatrixBase< T > & M2 ) [pure virtual, inherited]

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

Implemented in libMesh::DenseMatrix< Number >.

template<typename T>

void libMesh::DenseMatrix< T >::left_multiply ( const DenseMatrixBase< T > & M2 ) [inline, virtual]

Left multipliess by the matrix M2.

Definition at line 37 of file dense_matrix.C.

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

00038 {
00039   if (this->use_blas_lapack)
00040     this->_multiply_blas(M2, LEFT_MULTIPLY);
00041   else
00042     {
00043       // (*this) <- M2 * (*this)
00044       // Where:
00045       // (*this) = (m x n),
00046       // M2      = (m x p),
00047       // M3      = (p x n)
00048 
00049       // M3 is a copy of *this before it gets resize()d
00050       DenseMatrix<T> M3(*this);
00051 
00052       // Resize *this so that the result can fit
00053       this->resize (M2.m(), M3.n());
00054 
00055       // Call the multiply function in the base class
00056       this->multiply(*this, M2, M3);
00057     }
00058 }

template<typename T>

void libMesh::DenseMatrix< T >::left_multiply_transpose ( const DenseMatrix< T > & A ) [inline]

Left multiplies by the transpose of the matrix A.

Definition at line 64 of file dense_matrix.C.

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

Referenced by libMesh::DofMap::constrain_element_matrix(), libMesh::DofMap::constrain_element_matrix_and_vector(), and libMesh::DofMap::heterogenously_constrain_element_matrix_and_vector().

00065 {
00066   if (this->use_blas_lapack)
00067     this->_multiply_blas(A, LEFT_MULTIPLY_TRANSPOSE);
00068   else
00069     {
00070       //Check to see if we are doing (A^T)*A
00071       if (this == &A)
00072         {
00073           //libmesh_here();
00074           DenseMatrix<T> B(*this);
00075 
00076           // Simple but inefficient way
00077           // return this->left_multiply_transpose(B);
00078 
00079           // More efficient, but more code way
00080           // If A is mxn, the result will be a square matrix of Size n x n.
00081           const unsigned int n_rows = A.m();
00082           const unsigned int n_cols = A.n();
00083 
00084           // resize() *this and also zero out all entries.
00085           this->resize(n_cols,n_cols);
00086 
00087           // Compute the lower-triangular part
00088           for (unsigned int i=0; i<n_cols; ++i)
00089             for (unsigned int j=0; j<=i; ++j)
00090               for (unsigned int k=0; k<n_rows; ++k) // inner products are over n_rows
00091                 (*this)(i,j) += B(k,i)*B(k,j);
00092 
00093           // Copy lower-triangular part into upper-triangular part
00094           for (unsigned int i=0; i<n_cols; ++i)
00095             for (unsigned int j=i+1; j<n_cols; ++j)
00096               (*this)(i,j) = (*this)(j,i);
00097         }
00098 
00099       else
00100         {
00101           DenseMatrix<T> B(*this);
00102 
00103           this->resize (A.n(), B.n());
00104 
00105           libmesh_assert_equal_to (A.m(), B.m());
00106           libmesh_assert_equal_to (this->m(), A.n());
00107           libmesh_assert_equal_to (this->n(), B.n());
00108 
00109           const unsigned int m_s = A.n();
00110           const unsigned int p_s = A.m();
00111           const unsigned int n_s = this->n();
00112 
00113           // Do it this way because there is a
00114           // decent chance (at least for constraint matrices)
00115           // that A.transpose(i,k) = 0.
00116           for (unsigned int i=0; i<m_s; i++)
00117             for (unsigned int k=0; k<p_s; k++)
00118               if (A.transpose(i,k) != 0.)
00119                 for (unsigned int j=0; j<n_s; j++)
00120                   (*this)(i,j) += A.transpose(i,k)*B(k,j);
00121         }
00122     }
00123 
00124 }

template<typename T >

Real libMesh::DenseMatrix< T >::linfty_norm ( ) const [inline]

Return the linfty-norm of the matrix, that is $|M|_\infty=max_{all rows i}\sum_{all columns j} |M_ij|$ , (max. sum of rows). This is the natural matrix norm that is compatible to the linfty-norm of vectors, i.e. $|Mv|_\infty \leq |M|_\infty |v|_\infty$ .

Definition at line 820 of file dense_matrix.h.

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

00821 {
00822   libmesh_assert (this->_m);
00823   libmesh_assert (this->_n);
00824 
00825   Real rowsum = 0.;
00826   for (unsigned int j=0; j!=this->_n; j++)
00827     {
00828       rowsum += std::abs((*this)(0,j));
00829     }
00830   Real my_max = rowsum;
00831   for (unsigned int i=1; i!=this->_m; i++)
00832     {
00833       rowsum = 0.;
00834       for (unsigned int j=0; j!=this->_n; j++)
00835         {
00836           rowsum += std::abs((*this)(i,j));
00837         }
00838       my_max = (my_max > rowsum? my_max : rowsum);
00839     }
00840   return my_max;
00841 }

template<typename T>

void libMesh::DenseMatrix< T >::lu_solve	(	const DenseVector< T > &	b,
		DenseVector< T > &	x
	)			`[inline]`

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

Referenced by libMesh::WeightedPatchRecoveryErrorEstimator::EstimateError::operator()(), and libMesh::PatchRecoveryErrorEstimator::EstimateError::operator()().

00353 {
00354   // Check to be sure that the matrix is square before attempting
00355   // an LU-solve.  In general, one can compute the LU factorization of
00356   // a non-square matrix, but:
00357   //
00358   // Overdetermined systems (m>n) have a solution only if enough of
00359   // the equations are linearly-dependent.
00360   //
00361   // Underdetermined systems (m<n) typically have infinitely many
00362   // solutions.
00363   //
00364   // We don't want to deal with either of these ambiguous cases here...
00365   libmesh_assert_equal_to (this->m(), this->n());
00366 
00367   switch(this->_decomposition_type)
00368     {
00369     case NONE:
00370       {
00371         if (this->use_blas_lapack)
00372           this->_lu_decompose_lapack();
00373         else
00374           this->_lu_decompose ();
00375         break;
00376       }
00377 
00378     case LU_BLAS_LAPACK:
00379       {
00380         // Already factored, just need to call back_substitute.
00381         if (this->use_blas_lapack)
00382           break;
00383       }
00384 
00385     case LU:
00386       {
00387         // Already factored, just need to call back_substitute.
00388         if ( !(this->use_blas_lapack) )
00389           break;
00390       }
00391 
00392     default:
00393       {
00394         libMesh::err << "Error! This matrix already has a "
00395                       << "different decomposition..."
00396                       << std::endl;
00397         libmesh_error();
00398       }
00399     }
00400 
00401   if (this->use_blas_lapack)
00402      this->_lu_back_substitute_lapack (b, x);
00403   else
00404     this->_lu_back_substitute (b, x);
00405 }

template<typename T>

unsigned int libMesh::DenseMatrixBase< T >::m ( ) const [inline, inherited]

Returns:: the row-dimension of the matrix.

Definition at line 99 of file dense_matrix_base.h.

00099 { return _m; }

template<typename T >

Real libMesh::DenseMatrix< T >::max ( ) const [inline]

Returns:: the maximum element in the matrix. In case of complex numbers, this returns the maximum Real part.

Definition at line 772 of file dense_matrix.h.

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

00773 {
00774   libmesh_assert (this->_m);
00775   libmesh_assert (this->_n);
00776   Real my_max = libmesh_real((*this)(0,0));
00777 
00778   for (unsigned int i=0; i!=this->_m; i++)
00779     {
00780       for (unsigned int j=0; j!=this->_n; j++)
00781         {
00782           Real current = libmesh_real((*this)(i,j));
00783           my_max = (my_max > current? my_max : current);
00784         }
00785     }
00786   return my_max;
00787 }

template<typename T >

Real libMesh::DenseMatrix< T >::min ( ) const [inline]

Returns:: the minimum element in the matrix. In case of complex numbers, this returns the minimum Real part.

Definition at line 751 of file dense_matrix.h.

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

00752 {
00753   libmesh_assert (this->_m);
00754   libmesh_assert (this->_n);
00755   Real my_min = libmesh_real((*this)(0,0));
00756 
00757   for (unsigned int i=0; i!=this->_m; i++)
00758     {
00759       for (unsigned int j=0; j!=this->_n; j++)
00760         {
00761           Real current = libmesh_real((*this)(i,j));
00762           my_min = (my_min < current? my_min : current);
00763         }
00764     }
00765   return my_min;
00766 }

template<typename T>

void libMesh::DenseMatrixBase< T >::multiply	(	DenseMatrixBase< T > &	M1,
		const DenseMatrixBase< T > &	M2,
		const DenseMatrixBase< T > &	M3
	)			`[inline, protected, inherited]`

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

Referenced by libMesh::DenseMatrix< T >::left_multiply(), and libMesh::DenseMatrix< T >::right_multiply().

00034 {
00035   // Assertions to make sure we have been
00036   // passed matrices of the correct dimension.
00037   libmesh_assert_equal_to (M1.m(), M2.m());
00038   libmesh_assert_equal_to (M1.n(), M3.n());
00039   libmesh_assert_equal_to (M2.n(), M3.m());
00040 
00041   const unsigned int m_s = M2.m();
00042   const unsigned int p_s = M2.n();
00043   const unsigned int n_s = M1.n();
00044 
00045   // Do it this way because there is a
00046   // decent chance (at least for constraint matrices)
00047   // that M3(k,j) = 0. when right-multiplying.
00048   for (unsigned int k=0; k<p_s; k++)
00049     for (unsigned int j=0; j<n_s; j++)
00050       if (M3.el(k,j) != 0.)
00051         for (unsigned int i=0; i<m_s; i++)
00052           M1.el(i,j) += M2.el(i,k) * M3.el(k,j);
00053 }

template<typename T>

unsigned int libMesh::DenseMatrixBase< T >::n ( ) const [inline, inherited]

Returns:: the column-dimension of the matrix.

Definition at line 104 of file dense_matrix_base.h.

00104 { return _n; }

template<typename T>

bool libMesh::DenseMatrix< T >::operator!= ( const DenseMatrix< T > & mat ) const [inline]

Tests if mat is not exactly equal to this matrix.

Definition at line 714 of file dense_matrix.h.

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

00715 {
00716   for (unsigned int i=0; i<_val.size(); i++)
00717     if (_val[i] != mat._val[i])
00718       return true;
00719 
00720   return false;
00721 }

template<typename T >

T & libMesh::DenseMatrix< T >::operator()	(	const unsigned int	i,
		const unsigned int	j
	)			`[inline]`

Returns:: the (i,j) element of the matrix as a writeable reference.

Definition at line 654 of file dense_matrix.h.

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

00656 {
00657   libmesh_assert_less (i*j, _val.size());
00658   libmesh_assert_less (i, this->_m);
00659   libmesh_assert_less (j, this->_n);
00660 
00661   //return _val[(i) + (this->_m)*(j)]; // col-major
00662   return _val[(i)*(this->_n) + (j)]; // row-major
00663 }

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

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

00640 {
00641   libmesh_assert_less (i*j, _val.size());
00642   libmesh_assert_less (i, this->_m);
00643   libmesh_assert_less (j, this->_n);
00644 
00645 
00646   //  return _val[(i) + (this->_m)*(j)]; // col-major
00647   return _val[(i)*(this->_n) + (j)]; // row-major
00648 }

template<typename T>

DenseMatrix< T > & libMesh::DenseMatrix< T >::operator*= ( const T factor ) [inline]

Multiplies every element in the matrix by factor.

Definition at line 681 of file dense_matrix.h.

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

00682 {
00683   this->scale(factor);
00684   return *this;
00685 }

template<typename T>

DenseMatrix< T > & libMesh::DenseMatrix< T >::operator+= ( const DenseMatrix< T > & mat ) [inline]

Adds mat to this matrix.

Definition at line 727 of file dense_matrix.h.

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

00728 {
00729   for (unsigned int i=0; i<_val.size(); i++)
00730     _val[i] += mat._val[i];
00731 
00732   return *this;
00733 }

template<typename T>

DenseMatrix< T > & libMesh::DenseMatrix< T >::operator-= ( const DenseMatrix< T > & mat ) [inline]

Subtracts mat from this matrix.

Definition at line 739 of file dense_matrix.h.

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

00740 {
00741   for (unsigned int i=0; i<_val.size(); i++)
00742     _val[i] -= mat._val[i];
00743 
00744   return *this;
00745 }

template<typename T>

DenseMatrix< T > & libMesh::DenseMatrix< T >::operator= ( const DenseMatrix< T > & other_matrix ) [inline]

Assignment operator.

Definition at line 623 of file dense_matrix.h.

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

00624 {
00625   this->_m = other_matrix._m;
00626   this->_n = other_matrix._n;
00627 
00628   _val                = other_matrix._val;
00629   _decomposition_type = other_matrix._decomposition_type;
00630 
00631   return *this;
00632 }

template<typename T>

bool libMesh::DenseMatrix< T >::operator== ( const DenseMatrix< T > & mat ) const [inline]

Tests if mat is exactly equal to this matrix.

Definition at line 701 of file dense_matrix.h.

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

00702 {
00703   for (unsigned int i=0; i<_val.size(); i++)
00704     if (_val[i] != mat._val[i])
00705       return false;
00706 
00707   return true;
00708 }

template<typename T >

void libMesh::DenseMatrixBase< T >::print ( std::ostream & os = libMesh::out ) const [inline, inherited]

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

Definition at line 128 of file dense_matrix_base.C.

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

00129 {
00130   for (unsigned int i=0; i<this->m(); i++)
00131     {
00132       for (unsigned int j=0; j<this->n(); j++)
00133         os << std::setw(8)
00134            << this->el(i,j) << " ";
00135 
00136       os << std::endl;
00137     }
00138 
00139   return;
00140 }

template<typename T >

void libMesh::DenseMatrixBase< T >::print_scientific ( std::ostream & os ) const [inline, inherited]

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

Definition at line 86 of file dense_matrix_base.C.

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

00087 {
00088 #ifndef LIBMESH_BROKEN_IOSTREAM
00089 
00090   // save the initial format flags
00091   std::ios_base::fmtflags os_flags = os.flags();
00092 
00093   // Print the matrix entries.
00094   for (unsigned int i=0; i<this->m(); i++)
00095     {
00096       for (unsigned int j=0; j<this->n(); j++)
00097         os << std::setw(15)
00098            << std::scientific
00099            << std::setprecision(8)
00100            << this->el(i,j) << " ";
00101 
00102       os << std::endl;
00103     }
00104 
00105   // reset the original format flags
00106   os.flags(os_flags);
00107 
00108 #else
00109 
00110   // Print the matrix entries.
00111   for (unsigned int i=0; i<this->m(); i++)
00112     {
00113       for (unsigned int j=0; j<this->n(); j++)
00114         os << std::setprecision(8)
00115            << this->el(i,j)
00116            << " ";
00117 
00118       os << std::endl;
00119     }
00120 
00121 
00122 #endif
00123 }

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

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

Referenced by 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< T >::get_principal_submatrix(), libMesh::DenseMatrix< T >::get_transpose(), libMesh::DenseMatrix< T >::left_multiply(), libMesh::DenseMatrix< T >::left_multiply_transpose(), libMesh::BoundaryProjectSolution::operator()(), libMesh::ProjectFEMSolution::operator()(), libMesh::ProjectSolution::operator()(), libMesh::FEMContext::pre_fe_reinit(), libMesh::DenseMatrix< T >::right_multiply(), libMesh::DenseMatrix< T >::right_multiply_transpose(), and libMesh::HPCoarsenTest::select_refinement().

00592 {
00593   _val.resize(new_m*new_n);
00594 
00595   this->_m = new_m;
00596   this->_n = new_n;
00597 
00598   _decomposition_type = NONE;
00599   this->zero();
00600 }

template<typename T>

virtual void libMesh::DenseMatrixBase< T >::right_multiply ( const DenseMatrixBase< T > & M3 ) [pure virtual, inherited]

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

Implemented in libMesh::DenseMatrix< Number >.

template<typename T>

void libMesh::DenseMatrix< T >::right_multiply ( const DenseMatrixBase< T > & M3 ) [inline, virtual]

Right multiplies by the matrix M3.

Definition at line 132 of file dense_matrix.C.

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

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(), and libMesh::DofMap::heterogenously_constrain_element_matrix_and_vector().

00133 {
00134   if (this->use_blas_lapack)
00135     this->_multiply_blas(M3, RIGHT_MULTIPLY);
00136   else
00137     {
00138       // (*this) <- M3 * (*this)
00139       // Where:
00140       // (*this) = (m x n),
00141       // M2      = (m x p),
00142       // M3      = (p x n)
00143 
00144       // M2 is a copy of *this before it gets resize()d
00145       DenseMatrix<T> M2(*this);
00146 
00147       // Resize *this so that the result can fit
00148       this->resize (M2.m(), M3.n());
00149 
00150       this->multiply(*this, M2, M3);
00151     }
00152 }

template<typename T>

void libMesh::DenseMatrix< T >::right_multiply_transpose ( const DenseMatrix< T > & A ) [inline]

Right multiplies by the transpose of the matrix A

Definition at line 158 of file dense_matrix.C.

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

00159 {
00160   if (this->use_blas_lapack)
00161     this->_multiply_blas(B, RIGHT_MULTIPLY_TRANSPOSE);
00162   else
00163     {
00164       //Check to see if we are doing B*(B^T)
00165       if (this == &B)
00166         {
00167           //libmesh_here();
00168           DenseMatrix<T> A(*this);
00169 
00170           // Simple but inefficient way
00171           // return this->right_multiply_transpose(A);
00172 
00173           // More efficient, more code
00174           // If B is mxn, the result will be a square matrix of Size m x m.
00175           const unsigned int n_rows = B.m();
00176           const unsigned int n_cols = B.n();
00177 
00178           // resize() *this and also zero out all entries.
00179           this->resize(n_rows,n_rows);
00180 
00181           // Compute the lower-triangular part
00182           for (unsigned int i=0; i<n_rows; ++i)
00183             for (unsigned int j=0; j<=i; ++j)
00184               for (unsigned int k=0; k<n_cols; ++k) // inner products are over n_cols
00185                 (*this)(i,j) += A(i,k)*A(j,k);
00186 
00187           // Copy lower-triangular part into upper-triangular part
00188           for (unsigned int i=0; i<n_rows; ++i)
00189             for (unsigned int j=i+1; j<n_rows; ++j)
00190               (*this)(i,j) = (*this)(j,i);
00191         }
00192 
00193       else
00194         {
00195           DenseMatrix<T> A(*this);
00196 
00197           this->resize (A.m(), B.m());
00198 
00199           libmesh_assert_equal_to (A.n(), B.n());
00200           libmesh_assert_equal_to (this->m(), A.m());
00201           libmesh_assert_equal_to (this->n(), B.m());
00202 
00203           const unsigned int m_s = A.m();
00204           const unsigned int p_s = A.n();
00205           const unsigned int n_s = this->n();
00206 
00207           // Do it this way because there is a
00208           // decent chance (at least for constraint matrices)
00209           // that B.transpose(k,j) = 0.
00210           for (unsigned int j=0; j<n_s; j++)
00211             for (unsigned int k=0; k<p_s; k++)
00212               if (B.transpose(k,j) != 0.)
00213                 for (unsigned int i=0; i<m_s; i++)
00214                   (*this)(i,j) += A(i,k)*B.transpose(k,j);
00215         }
00216     }
00217 }

template<typename T>

void libMesh::DenseMatrix< T >::scale ( const T factor ) [inline]

Multiplies every element in the matrix by factor.

Definition at line 671 of file dense_matrix.h.

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

Referenced by libMesh::DenseMatrix< T >::operator*=().

00672 {
00673   for (unsigned int i=0; i<_val.size(); i++)
00674     _val[i] *= factor;
00675 }

template<typename T>

void libMesh::DenseMatrix< T >::svd	(	DenseVector< T > &	sigma,
		DenseMatrix< T > &	U,
		DenseMatrix< T > &	VT
	)			`[inline]`

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

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

00554 {
00555   // We use the LAPACK svd implementation
00556   _svd_lapack(sigma, U, VT);
00557 }

template<typename T>

void libMesh::DenseMatrix< T >::svd ( DenseVector< T > & sigma ) [inline]

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

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

00546 {
00547   // We use the LAPACK svd implementation
00548   _svd_lapack(sigma);
00549 }

template<typename T>

void libMesh::DenseMatrix< T >::swap ( DenseMatrix< T > & other_matrix ) [inline]

STL-like swap method

Definition at line 576 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< T >::_lu_back_substitute(), libMesh::DenseMatrix< T >::_multiply_blas(), libMesh::EulerSolver::element_residual(), libMesh::Euler2Solver::element_residual(), libMesh::EulerSolver::side_residual(), and libMesh::Euler2Solver::side_residual().

00577 {
00578   std::swap(this->_m, other_matrix._m);
00579   std::swap(this->_n, other_matrix._n);
00580   _val.swap(other_matrix._val);
00581   DecompositionType _temp = _decomposition_type;
00582   _decomposition_type = other_matrix._decomposition_type;
00583   other_matrix._decomposition_type = _temp;
00584 }

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

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

00849 {
00850   // Implement in terms of operator()
00851   return (*this)(j,i);
00852 }

template<typename T>

void libMesh::DenseMatrix< T >::vector_mult	(	DenseVector< T > &	dest,
		const DenseVector< T > &	arg
	)			const `[inline]`

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

Definition at line 222 of file dense_matrix.C.

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

Referenced by libMesh::DofMap::heterogenously_constrain_element_matrix_and_vector(), and libMesh::DenseMatrix< T >::vector_mult_add().

00224 {
00225   // Make sure the input sizes are compatible
00226   libmesh_assert_equal_to (this->n(), arg.size());
00227 
00228   // Resize and clear dest.
00229   // Note: DenseVector::resize() also zeros the vector.
00230   dest.resize(this->m());
00231 
00232   // Short-circuit if the matrix is empty
00233   if(this->m() == 0 || this->n() == 0)
00234     return;
00235 
00236   if (this->use_blas_lapack)
00237     this->_matvec_blas(1., 0., dest, arg);
00238   else
00239     {
00240       const unsigned int n_rows = this->m();
00241       const unsigned int n_cols = this->n();
00242 
00243       for(unsigned int i=0; i<n_rows; i++)
00244         for(unsigned int j=0; j<n_cols; j++)
00245           dest(i) += (*this)(i,j)*arg(j);
00246     }
00247 }

template<typename T>

void libMesh::DenseMatrix< T >::vector_mult_add	(	DenseVector< T > &	dest,
		const T	factor,
		const DenseVector< T > &	arg
	)			const `[inline]`

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

Definition at line 291 of file dense_matrix.C.

References libMesh::DenseMatrix< T >::_matvec_blas(), libMesh::DenseVector< T >::add(), libMesh::DenseMatrixBase< T >::m(), libMesh::DenseVector< T >::resize(), libMesh::DenseVector< T >::size(), libMesh::DenseMatrix< T >::use_blas_lapack, and libMesh::DenseMatrix< T >::vector_mult().

Referenced by libMesh::DofMap::build_constraint_matrix_and_vector().

00294 {
00295   // Short-circuit if the matrix is empty
00296   if(this->m() == 0)
00297   {
00298     dest.resize(0);
00299     return;
00300   }
00301 
00302   if (this->use_blas_lapack)
00303     this->_matvec_blas(factor, 1., dest, arg);
00304   else
00305     {
00306       DenseVector<T> temp(arg.size());
00307       this->vector_mult(temp, arg);
00308       dest.add(factor, temp);
00309     }
00310 }

template<typename T>

void libMesh::DenseMatrix< T >::vector_mult_transpose	(	DenseVector< T > &	dest,
		const DenseVector< T > &	arg
	)			const `[inline]`

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

Definition at line 253 of file dense_matrix.C.

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

Referenced by libMesh::DofMap::constrain_element_dyad_matrix(), libMesh::DofMap::constrain_element_matrix_and_vector(), libMesh::DofMap::constrain_element_vector(), and libMesh::DofMap::heterogenously_constrain_element_matrix_and_vector().

00255 {
00256   // Make sure the input sizes are compatible
00257   libmesh_assert_equal_to (this->m(), arg.size());
00258 
00259   // Resize and clear dest.
00260   // Note: DenseVector::resize() also zeros the vector.
00261   dest.resize(this->n());
00262 
00263   // Short-circuit if the matrix is empty
00264   if(this->m() == 0)
00265     return;
00266 
00267   if (this->use_blas_lapack)
00268     {
00269       this->_matvec_blas(1., 0., dest, arg, /*trans=*/true);
00270     }
00271   else
00272     {
00273       const unsigned int n_rows = this->m();
00274       const unsigned int n_cols = this->n();
00275 
00276       // WORKS
00277       // for(unsigned int j=0; j<n_cols; j++)
00278       //   for(unsigned int i=0; i<n_rows; i++)
00279       //     dest(j) += (*this)(i,j)*arg(i);
00280 
00281       // ALSO WORKS, (i,j) just swapped
00282       for(unsigned int i=0; i<n_cols; i++)
00283         for(unsigned int j=0; j<n_rows; j++)
00284           dest(i) += (*this)(j,i)*arg(j);
00285     }
00286 }

template<typename T >

void libMesh::DenseMatrix< T >::zero ( ) [inline, virtual]

Set every element in the matrix to 0.

Implements libMesh::DenseMatrixBase< T >.

Definition at line 606 of file dense_matrix.h.

References libMesh::DenseMatrix< T >::_decomposition_type, libMesh::DenseMatrix< T >::_val, and libMesh::DenseMatrix< T >::NONE.

Referenced by libMesh::HPCoarsenTest::add_projection(), libMesh::FEGenericBase< OutputType >::coarsened_dof_values(), libMesh::EulerSolver::element_residual(), libMesh::Euler2Solver::element_residual(), libMesh::BoundaryProjectSolution::operator()(), libMesh::ProjectFEMSolution::operator()(), libMesh::ProjectSolution::operator()(), libMesh::DenseMatrix< T >::resize(), libMesh::HPCoarsenTest::select_refinement(), libMesh::EulerSolver::side_residual(), and libMesh::Euler2Solver::side_residual().

00607 {
00608   _decomposition_type = NONE;
00609 
00610   // Just doing this ifdef to be completely safe
00611 #ifndef LIBMESH_USE_COMPLEX_NUMBERS
00612   if(_val.size())
00613     std::memset(&_val[0], 0, sizeof(T) * _val.size());
00614 #else
00615   std::fill (_val.begin(), _val.end(), 0.);
00616 #endif
00617 }

Friends And Related Function Documentation

template<typename T>

std::ostream& operator<<	(	std::ostream &	os,
		const DenseMatrixBase< T > &	m
	)			`[friend, inherited]`

Formatted print as above but allows you to do DenseMatrix K; libMesh::out << K << std::endl;

Definition at line 116 of file dense_matrix_base.h.

00117   {
00118     m.print(os);
00119     return os;
00120   }

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

Referenced by libMesh::DenseMatrix< T >::_cholesky_decompose(), libMesh::DenseMatrix< T >::_lu_decompose(), libMesh::DenseMatrix< T >::_lu_decompose_lapack(), libMesh::DenseMatrix< T >::cholesky_solve(), libMesh::DenseMatrix< T >::det(), libMesh::DenseMatrix< T >::lu_solve(), libMesh::DenseMatrix< T >::operator=(), libMesh::DenseMatrix< T >::resize(), libMesh::DenseMatrix< T >::swap(), and libMesh::DenseMatrix< T >::zero().

template<typename T>

unsigned int libMesh::DenseMatrixBase< T >::_m [protected, inherited]

The row dimension.

Definition at line 165 of file dense_matrix_base.h.

Referenced by libMesh::DenseMatrix< T >::_multiply_blas(), libMesh::DenseMatrixBase< T >::condense(), libMesh::DenseMatrix< T >::l1_norm(), libMesh::DenseMatrix< T >::linfty_norm(), libMesh::DenseMatrixBase< Number >::m(), libMesh::DenseMatrix< T >::max(), libMesh::DenseMatrix< T >::min(), libMesh::DenseMatrix< T >::operator()(), libMesh::DenseMatrix< T >::operator=(), libMesh::DenseMatrix< T >::resize(), and libMesh::DenseMatrix< T >::swap().

template<typename T>

unsigned int libMesh::DenseMatrixBase< T >::_n [protected, inherited]

The column dimension.

Definition at line 170 of file dense_matrix_base.h.

Referenced by libMesh::DenseMatrix< T >::_multiply_blas(), libMesh::DenseMatrix< T >::l1_norm(), libMesh::DenseMatrix< T >::linfty_norm(), libMesh::DenseMatrix< T >::max(), libMesh::DenseMatrix< T >::min(), libMesh::DenseMatrixBase< Number >::n(), libMesh::DenseMatrix< T >::operator()(), libMesh::DenseMatrix< T >::operator=(), libMesh::DenseMatrix< T >::resize(), and libMesh::DenseMatrix< T >::swap().

template<typename T>

std::vector<int> libMesh::DenseMatrix< T >::_pivots [private]

This array is 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 470 of file dense_matrix.h.

Referenced by libMesh::DenseMatrix< T >::_lu_back_substitute(), libMesh::DenseMatrix< T >::_lu_back_substitute_lapack(), libMesh::DenseMatrix< T >::_lu_decompose(), libMesh::DenseMatrix< T >::_lu_decompose_lapack(), and libMesh::DenseMatrix< T >::det().

template<typename T>

std::vector<T> libMesh::DenseMatrix< T >::_val [private]

The actual data values, stored as a 1D array.

Definition at line 357 of file dense_matrix.h.

template<typename T>

bool libMesh::DenseMatrix< T >::use_blas_lapack

Computes the inverse of the dense matrix (assuming it is invertible) by first computing the LU decomposition and then performing multiple back substitution steps. Follows the algorithm from Numerical Recipes in C available on the web. This is not the most memory efficient routine since the inverse is not computed "in place" but is instead placed into a the matrix inv passed in by the user. Run-time selectable option to turn on/off blas support. This was primarily used for testing purposes, and could be removed...

Definition at line 350 of file dense_matrix.h.

Referenced by libMesh::DenseMatrix< T >::det(), libMesh::DenseMatrix< T >::left_multiply(), libMesh::DenseMatrix< T >::left_multiply_transpose(), libMesh::DenseMatrix< T >::lu_solve(), libMesh::DenseMatrix< T >::right_multiply(), libMesh::DenseMatrix< T >::right_multiply_transpose(), libMesh::DenseMatrix< T >::vector_mult(), libMesh::DenseMatrix< T >::vector_mult_add(), and libMesh::DenseMatrix< T >::vector_mult_transpose().

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

libMesh::DenseMatrix< T > Class Template Reference

Public Member Functions

Public Attributes

Protected Member Functions

Protected Attributes

Private Types

Private Member Functions

Private Attributes

Friends

Detailed Description

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

Member Enumeration Documentation

Constructor & Destructor Documentation

Member Function Documentation

Friends And Related Function Documentation

Member Data Documentation

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