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 ()
 
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)
 
template<typename T2 >
DenseMatrix< T > & operator= (const DenseMatrix< T2 > &other_matrix)
 
void swap (DenseMatrix< T > &other_matrix)
 
void resize (const unsigned int new_m, const unsigned int new_n)
 
void scale (const T factor)
 
void scale_column (const unsigned int col, const T factor)
 
DenseMatrix< T > & operator*= (const T factor)
 
template<typename T2 , typename T3 >
boostcopy::enable_if_c
< ScalarTraits< T2 >::value,
void >::type 
add (const T2 factor, const DenseMatrix< T3 > &mat)
 
bool operator== (const DenseMatrix< T > &mat) const
 
bool operator!= (const DenseMatrix< T > &mat) const
 
DenseMatrix< T > & operator+= (const DenseMatrix< T > &mat)
 
DenseMatrix< T > & operator-= (const DenseMatrix< T > &mat)
 
Real min () const
 
Real max () const
 
Real l1_norm () const
 
Real linfty_norm () const
 
void left_multiply_transpose (const DenseMatrix< T > &A)
 
void right_multiply_transpose (const DenseMatrix< T > &A)
 
transpose (const unsigned int i, const unsigned int j) const
 
void get_transpose (DenseMatrix< T > &dest) const
 
std::vector< T > & get_values ()
 
const std::vector< T > & get_values () const
 
void condense (const unsigned int i, const unsigned int j, const T val, DenseVector< T > &rhs)
 
void lu_solve (const DenseVector< T > &b, DenseVector< T > &x)
 
template<typename T2 >
void cholesky_solve (const DenseVector< T2 > &b, DenseVector< T2 > &x)
 
void svd (DenseVector< T > &sigma)
 
void svd (DenseVector< T > &sigma, DenseMatrix< T > &U, DenseMatrix< T > &VT)
 
det ()
 
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
 

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

429  {
430  LEFT_MULTIPLY = 0,
434  };
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 417 of file dense_matrix.h.

417 {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 562 of file dense_matrix.h.

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

564  : DenseMatrixBase<T>(new_m,new_n),
565 #if defined(LIBMESH_HAVE_PETSC) && defined(LIBMESH_USE_REAL_NUMBERS) && defined(LIBMESH_DEFAULT_DOUBLE_PRECISION)
566  use_blas_lapack(true),
567 #else
568  use_blas_lapack(false),
569 #endif
570  _val(),
572  _pivots()
573 {
574  this->resize(new_m,new_n);
575 }
template<typename T>
virtual libMesh::DenseMatrix< T >::~DenseMatrix ( )
inlinevirtual

Copy-constructor. Destructor. Empty.

Definition at line 69 of file dense_matrix.h.

69 {}

Member Function Documentation

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

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

Definition at line 716 of file dense_matrix.C.

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

718 {
719  // Shorthand notation for number of rows and columns.
720  const unsigned int
721  n_rows = this->m(),
722  n_cols = this->n();
723 
724  // Just to be really sure...
725  libmesh_assert_equal_to (n_rows, n_cols);
726 
727  // A convenient reference to *this
728  const DenseMatrix<T>& A = *this;
729 
730  // Now compute the solution to Ax =b using the factorization.
731  x.resize(n_rows);
732 
733  // Solve for Ly=b
734  for (unsigned int i=0; i<n_cols; ++i)
735  {
736  T2 temp = b(i);
737 
738  for (unsigned int k=0; k<i; ++k)
739  temp -= A(i,k)*x(k);
740 
741  x(i) = temp / A(i,i);
742  }
743 
744  // Solve for L^T x = y
745  for (unsigned int i=0; i<n_cols; ++i)
746  {
747  const unsigned int ib = (n_cols-1)-i;
748 
749  for (unsigned int k=(ib+1); k<n_cols; ++k)
750  x(ib) -= A(k,ib) * x(k);
751 
752  x(ib) /= A(ib,ib);
753  }
754 }
template<typename T >
void libMesh::DenseMatrix< T >::_cholesky_decompose ( )
private

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

Definition at line 665 of file dense_matrix.C.

References A, and libMesh::err.

666 {
667  // If we called this function, there better not be any
668  // previous decomposition of the matrix.
669  libmesh_assert_equal_to (this->_decomposition_type, NONE);
670 
671  // Shorthand notation for number of rows and columns.
672  const unsigned int
673  n_rows = this->m(),
674  n_cols = this->n();
675 
676  // Just to be really sure...
677  libmesh_assert_equal_to (n_rows, n_cols);
678 
679  // A convenient reference to *this
680  DenseMatrix<T>& A = *this;
681 
682  for (unsigned int i=0; i<n_rows; ++i)
683  {
684  for (unsigned int j=i; j<n_cols; ++j)
685  {
686  for (unsigned int k=0; k<i; ++k)
687  A(i,j) -= A(i,k) * A(j,k);
688 
689  if (i == j)
690  {
691 #ifndef LIBMESH_USE_COMPLEX_NUMBERS
692  if (A(i,j) <= 0.0)
693  {
694  libMesh::err << "Error! Can only use Cholesky decomposition "
695  << "with symmetric positive definite matrices."
696  << std::endl;
697  libmesh_error();
698  }
699 #endif
700 
701  A(i,i) = std::sqrt(A(i,j));
702  }
703  else
704  A(j,i) = A(i,j) / A(i,i);
705  }
706  }
707 
708  // Set the flag for CHOLESKY decomposition
710 }
template<typename T>
void libMesh::DenseMatrix< T >::_lu_back_substitute ( const DenseVector< T > &  b,
DenseVector< T > &  x 
) const
private

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

Definition at line 413 of file dense_matrix.C.

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

415 {
416  const unsigned int
417  n_cols = this->n();
418 
419  libmesh_assert_equal_to (this->m(), n_cols);
420  libmesh_assert_equal_to (this->m(), b.size());
421 
422  x.resize (n_cols);
423 
424  // A convenient reference to *this
425  const DenseMatrix<T>& A = *this;
426 
427  // Temporary vector storage. We use this instead of
428  // modifying the RHS.
429  DenseVector<T> z = b;
430 
431  // Lower-triangular "top to bottom" solve step, taking into account pivots
432  for (unsigned int i=0; i<n_cols; ++i)
433  {
434  // Swap
435  if (_pivots[i] != static_cast<int>(i))
436  std::swap( z(i), z(_pivots[i]) );
437 
438  x(i) = z(i);
439 
440  for (unsigned int j=0; j<i; ++j)
441  x(i) -= A(i,j)*x(j);
442 
443  x(i) /= A(i,i);
444  }
445 
446  // Upper-triangular "bottom to top" solve step
447  const unsigned int last_row = n_cols-1;
448 
449  for (int i=last_row; i>=0; --i)
450  {
451  for (int j=i+1; j<static_cast<int>(n_cols); ++j)
452  x(i) -= A(i,j)*x(j);
453  }
454 }
template<typename T>
template void libMesh::DenseMatrix< T >::_lu_back_substitute_lapack ( const DenseVector< T > &  b,
DenseVector< T > &  x 
)
private

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

Definition at line 535 of file dense_matrix_blas_lapack.C.

References libMesh::DenseVector< T >::get_values(), and libMesh::out.

537 {
538  // The calling sequence for getrs is:
539  // dgetrs(TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
540 
541  // trans (input) char*
542  // 'n' for no tranpose, 't' for transpose
543  char TRANS[] = "t";
544 
545  // N (input) int*
546  // The order of the matrix A. N >= 0.
547  int N = this->m();
548 
549 
550  // NRHS (input) int*
551  // The number of right hand sides, i.e., the number of columns
552  // of the matrix B. NRHS >= 0.
553  int NRHS = 1;
554 
555  // A (input) DOUBLE PRECISION array, dimension (LDA,N)
556  // The factors L and U from the factorization A = P*L*U
557  // as computed by dgetrf.
558  // Here, we pass &(_val[0])
559 
560  // LDA (input) int*
561  // The leading dimension of the array A. LDA >= max(1,N).
562  int LDA = N;
563 
564  // ipiv (input) int array, dimension (N)
565  // The pivot indices from DGETRF; for 1<=i<=N, row i of the
566  // matrix was interchanged with row IPIV(i).
567  // Here, we pass &(_pivots[0]) which was computed in _lu_decompose_lapack
568 
569  // B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
570  // On entry, the right hand side matrix B.
571  // On exit, the solution matrix X.
572  // Here, we pass a copy of the rhs vector's data array in x, so that the
573  // passed right-hand side b is unmodified. I don't see a way around this
574  // copy if we want to maintain an unmodified rhs in LibMesh.
575  x = b;
576  std::vector<T>& x_vec = x.get_values();
577 
578  // We can avoid the copy if we don't care about overwriting the RHS: just
579  // pass b to the Lapack routine and then swap with x before exiting
580  // std::vector<T>& x_vec = b.get_values();
581 
582  // LDB (input) int*
583  // The leading dimension of the array B. LDB >= max(1,N).
584  int LDB = N;
585 
586  // INFO (output) int*
587  // = 0: successful exit
588  // < 0: if INFO = -i, the i-th argument had an illegal value
589  int INFO = 0;
590 
591  // Finally, ready to call the Lapack getrs function
592  LAPACKgetrs_(TRANS, &N, &NRHS, &(_val[0]), &LDA, &(_pivots[0]), &(x_vec[0]), &LDB, &INFO);
593 
594  // Check return value for errors
595  if (INFO != 0)
596  {
597  libMesh::out << "INFO="
598  << INFO
599  << ", Error during Lapack LU solve!" << std::endl;
600  libmesh_error();
601  }
602 
603  // Don't do this if you already made a copy of b above
604  // Swap b and x. The solution will then be in x, and whatever was originally
605  // in x, maybe garbage, maybe nothing, will be in b.
606  // FIXME: Rewrite the LU and Cholesky solves to just take one input, and overwrite
607  // the input. This *should* make user code simpler, as they don't have to create
608  // an extra vector just to pass it in to the solve function!
609  // b.swap(x);
610 }
template<typename T >
void libMesh::DenseMatrix< T >::_lu_decompose ( )
private

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

Definition at line 464 of file dense_matrix.C.

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

465 {
466  // If this function was called, there better not be any
467  // previous decomposition of the matrix.
468  libmesh_assert_equal_to (this->_decomposition_type, NONE);
469 
470  // Get the matrix size and make sure it is square
471  const unsigned int
472  n_rows = this->m();
473 
474  // A convenient reference to *this
475  DenseMatrix<T>& A = *this;
476 
477  _pivots.resize(n_rows);
478 
479  for (unsigned int i=0; i<n_rows; ++i)
480  {
481  // Find the pivot row by searching down the i'th column
482  _pivots[i] = i;
483 
484  // std::abs(complex) must return a Real!
485  Real the_max = std::abs( A(i,i) );
486  for (unsigned int j=i+1; j<n_rows; ++j)
487  {
488  Real candidate_max = std::abs( A(j,i) );
489  if (the_max < candidate_max)
490  {
491  the_max = candidate_max;
492  _pivots[i] = j;
493  }
494  }
495 
496  // libMesh::out << "the_max=" << the_max << " found at row " << _pivots[i] << std::endl;
497 
498  // If the max was found in a different row, interchange rows.
499  // Here we interchange the *entire* row, in Gaussian elimination
500  // you would only interchange the subrows A(i,j) and A(p(i),j), for j>i
501  if (_pivots[i] != static_cast<int>(i))
502  {
503  for (unsigned int j=0; j<n_rows; ++j)
504  std::swap( A(i,j), A(_pivots[i], j) );
505  }
506 
507 
508  // If the max abs entry found is zero, the matrix is singular
509  if (A(i,i) == libMesh::zero)
510  {
511  libMesh::out << "Matrix A is singular!" << std::endl;
512  libmesh_error();
513  }
514 
515  // Scale upper triangle entries of row i by the diagonal entry
516  // Note: don't scale the diagonal entry itself!
517  const T diag_inv = 1. / A(i,i);
518  for (unsigned int j=i+1; j<n_rows; ++j)
519  A(i,j) *= diag_inv;
520 
521  // Update the remaining sub-matrix A[i+1:m][i+1:m]
522  // by subtracting off (the diagonal-scaled)
523  // upper-triangular part of row i, scaled by the
524  // i'th column entry of each row. In terms of
525  // row operations, this is:
526  // for each r > i
527  // SubRow(r) = SubRow(r) - A(r,i)*SubRow(i)
528  //
529  // If we were scaling the i'th column as well, like
530  // in Gaussian elimination, this would 'zero' the
531  // entry in the i'th column.
532  for (unsigned int row=i+1; row<n_rows; ++row)
533  for (unsigned int col=i+1; col<n_rows; ++col)
534  A(row,col) -= A(row,i) * A(i,col);
535 
536  } // end i loop
537 
538  // Set the flag for LU decomposition
539  this->_decomposition_type = LU;
540 }
template<typename T >
template void libMesh::DenseMatrix< T >::_lu_decompose_lapack ( )
private

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

Definition at line 215 of file dense_matrix_blas_lapack.C.

References std::min(), and libMesh::out.

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

Uses the BLAS GEMV function (through PETSc) to compute

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

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

[ Implementation in dense_matrix_blas_lapack.C ]

Definition at line 631 of file dense_matrix_blas_lapack.C.

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

635 {
636  // Ensure that dest and arg sizes are compatible
637  if (!trans)
638  {
639  // dest ~ A * arg
640  // (mx1) (mxn) * (nx1)
641  if ((dest.size() != this->m()) || (arg.size() != this->n()))
642  {
643  libMesh::out << "Improper input argument sizes!" << std::endl;
644  libmesh_error();
645  }
646  }
647 
648  else // trans == true
649  {
650  // Ensure that dest and arg are proper size
651  // dest ~ A^T * arg
652  // (nx1) (nxm) * (mx1)
653  if ((dest.size() != this->n()) || (arg.size() != this->m()))
654  {
655  libMesh::out << "Improper input argument sizes!" << std::endl;
656  libmesh_error();
657  }
658  }
659 
660  // Calling sequence for dgemv:
661  //
662  // dgemv(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
663 
664  // TRANS - CHARACTER*1, 't' for transpose, 'n' for non-transpose multiply
665  // We store everything in row-major order, so pass the transpose flag for
666  // non-transposed matvecs and the 'n' flag for transposed matvecs
667  char TRANS[] = "t";
668  if (trans)
669  TRANS[0] = 'n';
670 
671  // M - INTEGER.
672  // On entry, M specifies the number of rows of the matrix A.
673  // In C/C++, pass the number of *cols* of A
674  int M = this->n();
675 
676  // N - INTEGER.
677  // On entry, N specifies the number of columns of the matrix A.
678  // In C/C++, pass the number of *rows* of A
679  int N = this->m();
680 
681  // ALPHA - DOUBLE PRECISION.
682  // The scalar constant passed to this function
683 
684  // A - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
685  // Before entry, the leading m by n part of the array A must
686  // contain the matrix of coefficients.
687  // The matrix, *this. Note that _matvec_blas is called from
688  // a const function, vector_mult(), and so we have made this function const
689  // as well. Since BLAS knows nothing about const, we have to cast it away
690  // now.
691  DenseMatrix<T>& a_ref = const_cast< DenseMatrix<T>& > ( *this );
692  std::vector<T>& a = a_ref.get_values();
693 
694  // LDA - INTEGER.
695  // On entry, LDA specifies the first dimension of A as declared
696  // in the calling (sub) program. LDA must be at least
697  // max( 1, m ).
698  int LDA = M;
699 
700  // X - DOUBLE PRECISION array of DIMENSION at least
701  // ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
702  // and at least
703  // ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
704  // Before entry, the incremented array X must contain the
705  // vector x.
706  // Here, we must cast away the const-ness of "arg" since BLAS knows
707  // nothing about const
708  DenseVector<T>& x_ref = const_cast< DenseVector<T>& > ( arg );
709  std::vector<T>& x = x_ref.get_values();
710 
711  // INCX - INTEGER.
712  // On entry, INCX specifies the increment for the elements of
713  // X. INCX must not be zero.
714  int INCX = 1;
715 
716  // BETA - DOUBLE PRECISION.
717  // On entry, BETA specifies the scalar beta. When BETA is
718  // supplied as zero then Y need not be set on input.
719  // The second scalar constant passed to this function
720 
721  // Y - DOUBLE PRECISION array of DIMENSION at least
722  // ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
723  // and at least
724  // ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
725  // Before entry with BETA non-zero, the incremented array Y
726  // must contain the vector y. On exit, Y is overwritten by the
727  // updated vector y.
728  // The input vector "dest"
729  std::vector<T>& y = dest.get_values();
730 
731  // INCY - INTEGER.
732  // On entry, INCY specifies the increment for the elements of
733  // Y. INCY must not be zero.
734  int INCY = 1;
735 
736  // Finally, ready to call the BLAS function
737  BLASgemv_(TRANS, &M, &N, &alpha, &(a[0]), &LDA, &(x[0]), &INCX, &beta, &(y[0]), &INCY);
738 }
template<typename T>
template void libMesh::DenseMatrix< T >::_multiply_blas ( const DenseMatrixBase< T > &  other,
_BLAS_Multiply_Flag  flag 
)
private

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

Definition at line 40 of file dense_matrix_blas_lapack.C.

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

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

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

290 {
291  // The calling sequence for dgetrf is:
292  // DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, INFO )
293 
294 
295  // JOBU (input) CHARACTER*1
296  // Specifies options for computing all or part of the matrix U:
297  // = 'A': all M columns of U are returned in array U:
298  // = 'S': the first min(m,n) columns of U (the left singular
299  // vectors) are returned in the array U;
300  // = 'O': the first min(m,n) columns of U (the left singular
301  // vectors) are overwritten on the array A;
302  // = 'N': no columns of U (no left singular vectors) are
303  // computed.
304  char JOBU = 'N';
305 
306  // JOBVT (input) CHARACTER*1
307  // Specifies options for computing all or part of the matrix
308  // V**T:
309  // = 'A': all N rows of V**T are returned in the array VT;
310  // = 'S': the first min(m,n) rows of V**T (the right singular
311  // vectors) are returned in the array VT;
312  // = 'O': the first min(m,n) rows of V**T (the right singular
313  // vectors) are overwritten on the array A;
314  // = 'N': no rows of V**T (no right singular vectors) are
315  // computed.
316  char JOBVT = 'N';
317 
318  std::vector<T> sigma_val;
319  std::vector<T> U_val;
320  std::vector<T> VT_val;
321 
322  _svd_helper(JOBU, JOBVT, sigma_val, U_val, VT_val);
323 
324  // Load the singular values into sigma, ignore U_val and VT_val
325  const unsigned int n_sigma_vals =
326  libmesh_cast_int<unsigned int>(sigma_val.size());
327  sigma.resize(n_sigma_vals);
328  for(unsigned int i=0; i<n_sigma_vals; i++)
329  sigma(i) = sigma_val[i];
330 
331 }
template<typename T>
void libMesh::DenseMatrix< T >::_svd_lapack ( DenseVector< T > &  sigma,
DenseMatrix< T > &  U,
DenseMatrix< T > &  VT 
)
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::DenseMatrixBase< T >::m(), libMesh::DenseMatrixBase< T >::n(), libMesh::DenseVector< T >::resize(), and libMesh::DenseMatrix< T >::resize().

335 {
336  // The calling sequence for dgetrf is:
337  // DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, INFO )
338 
339 
340  // JOBU (input) CHARACTER*1
341  // Specifies options for computing all or part of the matrix U:
342  // = 'A': all M columns of U are returned in array U:
343  // = 'S': the first min(m,n) columns of U (the left singular
344  // vectors) are returned in the array U;
345  // = 'O': the first min(m,n) columns of U (the left singular
346  // vectors) are overwritten on the array A;
347  // = 'N': no columns of U (no left singular vectors) are
348  // computed.
349  char JOBU = 'S';
350 
351  // JOBVT (input) CHARACTER*1
352  // Specifies options for computing all or part of the matrix
353  // V**T:
354  // = 'A': all N rows of V**T are returned in the array VT;
355  // = 'S': the first min(m,n) rows of V**T (the right singular
356  // vectors) are returned in the array VT;
357  // = 'O': the first min(m,n) rows of V**T (the right singular
358  // vectors) are overwritten on the array A;
359  // = 'N': no rows of V**T (no right singular vectors) are
360  // computed.
361  char JOBVT = 'S';
362 
363  std::vector<T> sigma_val;
364  std::vector<T> U_val;
365  std::vector<T> VT_val;
366 
367  _svd_helper(JOBU, JOBVT, sigma_val, U_val, VT_val);
368 
369  // Load the singular values into sigma, ignore U_val and VT_val
370  const unsigned int n_sigma_vals =
371  libmesh_cast_int<unsigned int>(sigma_val.size());
372  sigma.resize(n_sigma_vals);
373  for(unsigned int i=0; i<n_sigma_vals; i++)
374  sigma(i) = sigma_val[i];
375 
376  int M = this->n();
377  int N = this->m();
378  int min_MN = (M < N) ? M : N;
379  U.resize(M,min_MN);
380  for(unsigned int i=0; i<U.m(); i++)
381  for(unsigned int j=0; j<U.n(); j++)
382  {
383  unsigned int index = i + j*U.n(); // Column major storage
384  U(i,j) = U_val[index];
385  }
386 
387  VT.resize(min_MN,N);
388  for(unsigned int i=0; i<VT.m(); i++)
389  for(unsigned int j=0; j<VT.n(); j++)
390  {
391  unsigned int index = i + j*U.n(); // Column major storage
392  VT(i,j) = VT_val[index];
393  }
394 
395 }
template<typename T >
template<typename T2 , typename T3 >
boostcopy::enable_if_c< ScalarTraits< T2 >::value, void >::type libMesh::DenseMatrixBase< T >::add ( const T2  factor,
const DenseMatrixBase< T3 > &  mat 
)
inlineinherited

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

Definition at line 184 of file dense_matrix_base.h.

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

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

Adds factor times mat to this matrix.

Definition at line 737 of file dense_matrix.h.

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

739 {
740  libmesh_assert_equal_to (this->m(), mat.m());
741  libmesh_assert_equal_to (this->n(), mat.n());
742 
743  for (unsigned int i=0; i<this->m(); i++)
744  for (unsigned int j=0; j<this->n(); j++)
745  (*this)(i,j) += factor * mat(i,j);
746 }
template<typename T >
template<typename T2 >
template void libMesh::DenseMatrix< T >::cholesky_solve ( const DenseVector< T2 > &  b,
DenseVector< T2 > &  x 
)

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

Definition at line 628 of file dense_matrix.C.

References libMesh::err.

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

630 {
631  // Check for a previous decomposition
632  switch(this->_decomposition_type)
633  {
634  case NONE:
635  {
636  this->_cholesky_decompose ();
637  break;
638  }
639 
640  case CHOLESKY:
641  {
642  // Already factored, just need to call back_substitute.
643  break;
644  }
645 
646  default:
647  {
648  libMesh::err << "Error! This matrix already has a "
649  << "different decomposition..."
650  << std::endl;
651  libmesh_error();
652  }
653  }
654 
655  // Perform back substitution
656  this->_cholesky_back_substitute (b, x);
657 }
template<typename T >
void libMesh::DenseMatrixBase< T >::condense ( const unsigned int  i,
const unsigned int  j,
const T  val,
DenseVectorBase< T > &  rhs 
)
protectedinherited

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

Definition at line 58 of file dense_matrix_base.C.

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

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

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

295  { DenseMatrixBase<T>::condense (i, j, val, rhs); }
template<typename T >
T libMesh::DenseMatrix< T >::det ( )
Returns
the determinant of the matrix. Note that this means doing an LU decomposition and then computing the product of the diagonal terms. Therefore this is a non-const method.

Definition at line 562 of file dense_matrix.C.

References libMesh::err, and libMesh::Real.

563 {
564  switch(this->_decomposition_type)
565  {
566  case NONE:
567  {
568  // First LU decompose the matrix.
569  // Note that the lu_decompose routine will check to see if the
570  // matrix is square so we don't worry about it.
571  if (this->use_blas_lapack)
572  this->_lu_decompose_lapack();
573  else
574  this->_lu_decompose ();
575  }
576  case LU:
577  case LU_BLAS_LAPACK:
578  {
579  // Already decomposed, don't do anything
580  break;
581  }
582  default:
583  {
584  libMesh::err << "Error! Can't compute the determinant under "
585  << "the current decomposition."
586  << std::endl;
587  libmesh_error();
588  }
589  }
590 
591  // A variable to keep track of the running product of diagonal terms.
592  T determinant = 1.;
593 
594  // Loop over diagonal terms, computing the product. In practice,
595  // be careful because this value could easily become too large to
596  // fit in a double or float. To be safe, one should keep track of
597  // the power (of 10) of the determinant in a separate variable
598  // and maintain an order 1 value for the determinant itself.
599  unsigned int n_interchanges = 0;
600  for (unsigned int i=0; i<this->m(); i++)
601  {
602  if (this->_decomposition_type==LU)
603  if (_pivots[i] != static_cast<int>(i))
604  n_interchanges++;
605 
606  // Lapack pivots are 1-based!
608  if (_pivots[i] != static_cast<int>(i+1))
609  n_interchanges++;
610 
611  determinant *= (*this)(i,i);
612  }
613 
614  // Compute sign of determinant, depends on number of row interchanges!
615  // The sign should be (-1)^{n}, where n is the number of interchanges.
616  Real sign = n_interchanges % 2 == 0 ? 1. : -1.;
617 
618  return sign*determinant;
619 }
template<typename T>
virtual T libMesh::DenseMatrix< T >::el ( const unsigned int  i,
const unsigned int  j 
) const
inlinevirtual
Returns
the (i,j) element of the matrix as a writeable reference.

Implements libMesh::DenseMatrixBase< T >.

Definition at line 92 of file dense_matrix.h.

93  { return (*this)(i,j); }
template<typename T>
virtual T& libMesh::DenseMatrix< T >::el ( const unsigned int  i,
const unsigned int  j 
)
inlinevirtual
Returns
the (i,j) element of the matrix as a writeable reference.

Implements libMesh::DenseMatrixBase< T >.

Definition at line 98 of file dense_matrix.h.

99  { return (*this)(i,j); }
template<typename T>
void libMesh::DenseMatrix< T >::get_principal_submatrix ( unsigned int  sub_m,
unsigned int  sub_n,
DenseMatrix< T > &  dest 
) const

Put the sub_m x sub_n principal submatrix into dest.

Definition at line 315 of file dense_matrix.C.

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

318 {
319  libmesh_assert( (sub_m <= this->m()) && (sub_n <= this->n()) );
320 
321  dest.resize(sub_m, sub_n);
322  for(unsigned int i=0; i<sub_m; i++)
323  for(unsigned int j=0; j<sub_n; j++)
324  dest(i,j) = (*this)(i,j);
325 }
template<typename T>
void libMesh::DenseMatrix< T >::get_principal_submatrix ( unsigned int  sub_m,
DenseMatrix< T > &  dest 
) const

Put the sub_m x sub_m principal submatrix into dest.

Definition at line 330 of file dense_matrix.C.

331 {
332  get_principal_submatrix(sub_m, sub_m, dest);
333 }
template<typename T>
void libMesh::DenseMatrix< T >::get_transpose ( DenseMatrix< T > &  dest) const

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

339 {
340  dest.resize(this->n(), this->m());
341 
342  for (unsigned int i=0; i<dest.m(); i++)
343  for (unsigned int j=0; j<dest.n(); j++)
344  dest(i,j) = (*this)(j,i);
345 }
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 278 of file dense_matrix.h.

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

278 { return _val; }
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 283 of file dense_matrix.h.

283 { 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 844 of file dense_matrix.h.

References std::abs(), libMesh::libmesh_assert(), and libMesh::Real.

845 {
846  libmesh_assert (this->_m);
847  libmesh_assert (this->_n);
848 
849  Real columnsum = 0.;
850  for (unsigned int i=0; i!=this->_m; i++)
851  {
852  columnsum += std::abs((*this)(i,0));
853  }
854  Real my_max = columnsum;
855  for (unsigned int j=1; j!=this->_n; j++)
856  {
857  columnsum = 0.;
858  for (unsigned int i=0; i!=this->_m; i++)
859  {
860  columnsum += std::abs((*this)(i,j));
861  }
862  my_max = (my_max > columnsum? my_max : columnsum);
863  }
864  return my_max;
865 }
template<typename T>
void libMesh::DenseMatrix< T >::left_multiply ( const DenseMatrixBase< T > &  M2)
virtual

Left multipliess by the matrix M2.

Implements libMesh::DenseMatrixBase< T >.

Definition at line 37 of file dense_matrix.C.

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

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

Left multiplies by the transpose of the matrix A.

Definition at line 64 of file dense_matrix.C.

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

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

65 {
66  if (this->use_blas_lapack)
68  else
69  {
70  //Check to see if we are doing (A^T)*A
71  if (this == &A)
72  {
73  //libmesh_here();
74  DenseMatrix<T> B(*this);
75 
76  // Simple but inefficient way
77  // return this->left_multiply_transpose(B);
78 
79  // More efficient, but more code way
80  // If A is mxn, the result will be a square matrix of Size n x n.
81  const unsigned int n_rows = A.m();
82  const unsigned int n_cols = A.n();
83 
84  // resize() *this and also zero out all entries.
85  this->resize(n_cols,n_cols);
86 
87  // Compute the lower-triangular part
88  for (unsigned int i=0; i<n_cols; ++i)
89  for (unsigned int j=0; j<=i; ++j)
90  for (unsigned int k=0; k<n_rows; ++k) // inner products are over n_rows
91  (*this)(i,j) += B(k,i)*B(k,j);
92 
93  // Copy lower-triangular part into upper-triangular part
94  for (unsigned int i=0; i<n_cols; ++i)
95  for (unsigned int j=i+1; j<n_cols; ++j)
96  (*this)(i,j) = (*this)(j,i);
97  }
98 
99  else
100  {
101  DenseMatrix<T> B(*this);
102 
103  this->resize (A.n(), B.n());
104 
105  libmesh_assert_equal_to (A.m(), B.m());
106  libmesh_assert_equal_to (this->m(), A.n());
107  libmesh_assert_equal_to (this->n(), B.n());
108 
109  const unsigned int m_s = A.n();
110  const unsigned int p_s = A.m();
111  const unsigned int n_s = this->n();
112 
113  // Do it this way because there is a
114  // decent chance (at least for constraint matrices)
115  // that A.transpose(i,k) = 0.
116  for (unsigned int i=0; i<m_s; i++)
117  for (unsigned int k=0; k<p_s; k++)
118  if (A.transpose(i,k) != 0.)
119  for (unsigned int j=0; j<n_s; j++)
120  (*this)(i,j) += A.transpose(i,k)*B(k,j);
121  }
122  }
123 
124 }
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 871 of file dense_matrix.h.

References std::abs(), libMesh::libmesh_assert(), and libMesh::Real.

872 {
873  libmesh_assert (this->_m);
874  libmesh_assert (this->_n);
875 
876  Real rowsum = 0.;
877  for (unsigned int j=0; j!=this->_n; j++)
878  {
879  rowsum += std::abs((*this)(0,j));
880  }
881  Real my_max = rowsum;
882  for (unsigned int i=1; i!=this->_m; i++)
883  {
884  rowsum = 0.;
885  for (unsigned int j=0; j!=this->_n; j++)
886  {
887  rowsum += std::abs((*this)(i,j));
888  }
889  my_max = (my_max > rowsum? my_max : rowsum);
890  }
891  return my_max;
892 }
template<typename T>
void libMesh::DenseMatrix< T >::lu_solve ( const DenseVector< T > &  b,
DenseVector< T > &  x 
)

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

Definition at line 351 of file dense_matrix.C.

References libMesh::err.

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

353 {
354  // Check to be sure that the matrix is square before attempting
355  // an LU-solve. In general, one can compute the LU factorization of
356  // a non-square matrix, but:
357  //
358  // Overdetermined systems (m>n) have a solution only if enough of
359  // the equations are linearly-dependent.
360  //
361  // Underdetermined systems (m<n) typically have infinitely many
362  // solutions.
363  //
364  // We don't want to deal with either of these ambiguous cases here...
365  libmesh_assert_equal_to (this->m(), this->n());
366 
367  switch(this->_decomposition_type)
368  {
369  case NONE:
370  {
371  if (this->use_blas_lapack)
372  this->_lu_decompose_lapack();
373  else
374  this->_lu_decompose ();
375  break;
376  }
377 
378  case LU_BLAS_LAPACK:
379  {
380  // Already factored, just need to call back_substitute.
381  if (this->use_blas_lapack)
382  break;
383  }
384 
385  case LU:
386  {
387  // Already factored, just need to call back_substitute.
388  if ( !(this->use_blas_lapack) )
389  break;
390  }
391 
392  default:
393  {
394  libMesh::err << "Error! This matrix already has a "
395  << "different decomposition..."
396  << std::endl;
397  libmesh_error();
398  }
399  }
400 
401  if (this->use_blas_lapack)
402  this->_lu_back_substitute_lapack (b, x);
403  else
404  this->_lu_back_substitute (b, x);
405 }
template<typename T>
unsigned int libMesh::DenseMatrixBase< T >::m ( ) const
inlineinherited
Returns
the row-dimension of the matrix.

Definition at line 99 of file dense_matrix_base.h.

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

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

99 { 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 823 of file dense_matrix.h.

References libMesh::libmesh_assert(), libMesh::libmesh_real(), and libMesh::Real.

824 {
825  libmesh_assert (this->_m);
826  libmesh_assert (this->_n);
827  Real my_max = libmesh_real((*this)(0,0));
828 
829  for (unsigned int i=0; i!=this->_m; i++)
830  {
831  for (unsigned int j=0; j!=this->_n; j++)
832  {
833  Real current = libmesh_real((*this)(i,j));
834  my_max = (my_max > current? my_max : current);
835  }
836  }
837  return my_max;
838 }
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 802 of file dense_matrix.h.

References libMesh::libmesh_assert(), libMesh::libmesh_real(), and libMesh::Real.

803 {
804  libmesh_assert (this->_m);
805  libmesh_assert (this->_n);
806  Real my_min = libmesh_real((*this)(0,0));
807 
808  for (unsigned int i=0; i!=this->_m; i++)
809  {
810  for (unsigned int j=0; j!=this->_n; j++)
811  {
812  Real current = libmesh_real((*this)(i,j));
813  my_min = (my_min < current? my_min : current);
814  }
815  }
816  return my_min;
817 }
template<typename T >
void libMesh::DenseMatrixBase< T >::multiply ( DenseMatrixBase< T > &  M1,
const DenseMatrixBase< T > &  M2,
const DenseMatrixBase< T > &  M3 
)
protectedinherited

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

Definition at line 31 of file dense_matrix_base.C.

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

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

Definition at line 104 of file dense_matrix_base.h.

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

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

104 { 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 765 of file dense_matrix.h.

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

766 {
767  for (unsigned int i=0; i<_val.size(); i++)
768  if (_val[i] != mat._val[i])
769  return true;
770 
771  return false;
772 }
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 672 of file dense_matrix.h.

674 {
675  libmesh_assert_less (i*j, _val.size());
676  libmesh_assert_less (i, this->_m);
677  libmesh_assert_less (j, this->_n);
678 
679 
680  // return _val[(i) + (this->_m)*(j)]; // col-major
681  return _val[(i)*(this->_n) + (j)]; // row-major
682 }
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 688 of file dense_matrix.h.

690 {
691  libmesh_assert_less (i*j, _val.size());
692  libmesh_assert_less (i, this->_m);
693  libmesh_assert_less (j, this->_n);
694 
695  //return _val[(i) + (this->_m)*(j)]; // col-major
696  return _val[(i)*(this->_n) + (j)]; // row-major
697 }
template<typename T>
DenseMatrix< T > & libMesh::DenseMatrix< T >::operator*= ( const T  factor)
inline

Multiplies every element in the matrix by factor.

Definition at line 724 of file dense_matrix.h.

References libMesh::MeshTools::Modification::scale().

725 {
726  this->scale(factor);
727  return *this;
728 }
template<typename T>
DenseMatrix< T > & libMesh::DenseMatrix< T >::operator+= ( const DenseMatrix< T > &  mat)
inline

Adds mat to this matrix.

Definition at line 778 of file dense_matrix.h.

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

779 {
780  for (unsigned int i=0; i<_val.size(); i++)
781  _val[i] += mat._val[i];
782 
783  return *this;
784 }
template<typename T>
DenseMatrix< T > & libMesh::DenseMatrix< T >::operator-= ( const DenseMatrix< T > &  mat)
inline

Subtracts mat from this matrix.

Definition at line 790 of file dense_matrix.h.

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

791 {
792  for (unsigned int i=0; i<_val.size(); i++)
793  _val[i] -= mat._val[i];
794 
795  return *this;
796 }
template<typename T>
DenseMatrix< T > & libMesh::DenseMatrix< T >::operator= ( const DenseMatrix< T > &  other_matrix)
inline

Assignment operator.

Definition at line 657 of file dense_matrix.h.

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

658 {
659  this->_m = other_matrix._m;
660  this->_n = other_matrix._n;
661 
662  _val = other_matrix._val;
663  _decomposition_type = other_matrix._decomposition_type;
664 
665  return *this;
666 }
template<typename T >
template<typename T2 >
DenseMatrix< T > & libMesh::DenseMatrix< T >::operator= ( const DenseMatrix< T2 > &  other_matrix)
inline

Assignment-from-other-matrix-type operator

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

Definition at line 609 of file dense_matrix.h.

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

610 {
611  unsigned int mat_m = mat.m(), mat_n = mat.n();
612  this->resize(mat_m, mat_n);
613  for (unsigned int i=0; i<mat_m; i++)
614  for (unsigned int j=0; j<mat_n; j++)
615  (*this)(i,j) = mat(i,j);
616 
617  return *this;
618 }
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 752 of file dense_matrix.h.

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

753 {
754  for (unsigned int i=0; i<_val.size(); i++)
755  if (_val[i] != mat._val[i])
756  return false;
757 
758  return true;
759 }
template<typename T >
void libMesh::DenseMatrixBase< T >::print ( std::ostream &  os = libMesh::out) const
inherited

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

Definition at line 128 of file dense_matrix_base.C.

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

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

Definition at line 86 of file dense_matrix_base.C.

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

Right multiplies by the matrix M3.

Implements libMesh::DenseMatrixBase< T >.

Definition at line 132 of file dense_matrix.C.

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

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

133 {
134  if (this->use_blas_lapack)
135  this->_multiply_blas(M3, RIGHT_MULTIPLY);
136  else
137  {
138  // (*this) <- M3 * (*this)
139  // Where:
140  // (*this) = (m x n),
141  // M2 = (m x p),
142  // M3 = (p x n)
143 
144  // M2 is a copy of *this before it gets resize()d
145  DenseMatrix<T> M2(*this);
146 
147  // Resize *this so that the result can fit
148  this->resize (M2.m(), M3.n());
149 
150  this->multiply(*this, M2, M3);
151  }
152 }
template<typename T>
void libMesh::DenseMatrix< T >::right_multiply_transpose ( const DenseMatrix< T > &  A)

Right multiplies by the transpose of the matrix A

Definition at line 158 of file dense_matrix.C.

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

159 {
160  if (this->use_blas_lapack)
162  else
163  {
164  //Check to see if we are doing B*(B^T)
165  if (this == &B)
166  {
167  //libmesh_here();
168  DenseMatrix<T> A(*this);
169 
170  // Simple but inefficient way
171  // return this->right_multiply_transpose(A);
172 
173  // More efficient, more code
174  // If B is mxn, the result will be a square matrix of Size m x m.
175  const unsigned int n_rows = B.m();
176  const unsigned int n_cols = B.n();
177 
178  // resize() *this and also zero out all entries.
179  this->resize(n_rows,n_rows);
180 
181  // Compute the lower-triangular part
182  for (unsigned int i=0; i<n_rows; ++i)
183  for (unsigned int j=0; j<=i; ++j)
184  for (unsigned int k=0; k<n_cols; ++k) // inner products are over n_cols
185  (*this)(i,j) += A(i,k)*A(j,k);
186 
187  // Copy lower-triangular part into upper-triangular part
188  for (unsigned int i=0; i<n_rows; ++i)
189  for (unsigned int j=i+1; j<n_rows; ++j)
190  (*this)(i,j) = (*this)(j,i);
191  }
192 
193  else
194  {
195  DenseMatrix<T> A(*this);
196 
197  this->resize (A.m(), B.m());
198 
199  libmesh_assert_equal_to (A.n(), B.n());
200  libmesh_assert_equal_to (this->m(), A.m());
201  libmesh_assert_equal_to (this->n(), B.m());
202 
203  const unsigned int m_s = A.m();
204  const unsigned int p_s = A.n();
205  const unsigned int n_s = this->n();
206 
207  // Do it this way because there is a
208  // decent chance (at least for constraint matrices)
209  // that B.transpose(k,j) = 0.
210  for (unsigned int j=0; j<n_s; j++)
211  for (unsigned int k=0; k<p_s; k++)
212  if (B.transpose(k,j) != 0.)
213  for (unsigned int i=0; i<m_s; i++)
214  (*this)(i,j) += A(i,k)*B.transpose(k,j);
215  }
216  }
217 }
template<typename T>
void libMesh::DenseMatrix< T >::scale ( const T  factor)
inline

Multiplies every element in the matrix by factor.

Definition at line 705 of file dense_matrix.h.

706 {
707  for (unsigned int i=0; i<_val.size(); i++)
708  _val[i] *= factor;
709 }
template<typename T>
void libMesh::DenseMatrix< T >::scale_column ( const unsigned int  col,
const T  factor 
)
inline

Multiplies every element in the column col matrix by factor.

Definition at line 714 of file dense_matrix.h.

715 {
716  for (unsigned int i=0; i<this->m(); i++)
717  (*this)(i, col) *= factor;
718 }
template<typename T>
void libMesh::DenseMatrix< T >::svd ( DenseVector< T > &  sigma)

Compute the Singular Value Decomposition of the matrix. On exit, sigma holds all of the singular values (in descending order).

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

Definition at line 545 of file dense_matrix.C.

546 {
547  // We use the LAPACK svd implementation
548  _svd_lapack(sigma);
549 }
template<typename T>
void libMesh::DenseMatrix< T >::svd ( DenseVector< T > &  sigma,
DenseMatrix< T > &  U,
DenseMatrix< T > &  VT 
)

Compute the "reduced" Singular Value Decomposition of the matrix. On exit, sigma holds all of the singular values (in descending order), U holds the left singular vectors, and VT holds the transpose of the right singular vectors. In the reduced SVD, U has min(m,n) columns and VT has min(m,n) rows. (In the "full" SVD, U and VT would be square.)

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

Definition at line 553 of file dense_matrix.C.

554 {
555  // We use the LAPACK svd implementation
556  _svd_lapack(sigma, U, VT);
557 }
template<typename T>
void libMesh::DenseMatrix< T >::swap ( DenseMatrix< T > &  other_matrix)
inline

STL-like swap method

Definition at line 594 of file dense_matrix.h.

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

Referenced by libMesh::EulerSolver::element_residual(), libMesh::Euler2Solver::element_residual(), libMesh::EulerSolver::side_residual(), and libMesh::Euler2Solver::side_residual().

595 {
596  std::swap(this->_m, other_matrix._m);
597  std::swap(this->_n, other_matrix._n);
598  _val.swap(other_matrix._val);
600  _decomposition_type = other_matrix._decomposition_type;
601  other_matrix._decomposition_type = _temp;
602 }
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 898 of file dense_matrix.h.

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

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

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

Definition at line 222 of file dense_matrix.C.

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

Referenced by libMesh::DofMap::heterogenously_constrain_element_matrix_and_vector(), and libMesh::DofMap::heterogenously_constrain_element_vector().

224 {
225  // Make sure the input sizes are compatible
226  libmesh_assert_equal_to (this->n(), arg.size());
227 
228  // Resize and clear dest.
229  // Note: DenseVector::resize() also zeros the vector.
230  dest.resize(this->m());
231 
232  // Short-circuit if the matrix is empty
233  if(this->m() == 0 || this->n() == 0)
234  return;
235 
236  if (this->use_blas_lapack)
237  this->_matvec_blas(1., 0., dest, arg);
238  else
239  {
240  const unsigned int n_rows = this->m();
241  const unsigned int n_cols = this->n();
242 
243  for(unsigned int i=0; i<n_rows; i++)
244  for(unsigned int j=0; j<n_cols; j++)
245  dest(i) += (*this)(i,j)*arg(j);
246  }
247 }
template<typename T>
void libMesh::DenseMatrix< T >::vector_mult_add ( DenseVector< T > &  dest,
const T  factor,
const DenseVector< T > &  arg 
) const

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

Definition at line 291 of file dense_matrix.C.

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

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

294 {
295  // Short-circuit if the matrix is empty
296  if(this->m() == 0)
297  {
298  dest.resize(0);
299  return;
300  }
301 
302  if (this->use_blas_lapack)
303  this->_matvec_blas(factor, 1., dest, arg);
304  else
305  {
306  DenseVector<T> temp(arg.size());
307  this->vector_mult(temp, arg);
308  dest.add(factor, temp);
309  }
310 }
template<typename T>
void libMesh::DenseMatrix< T >::vector_mult_transpose ( DenseVector< T > &  dest,
const DenseVector< T > &  arg 
) const

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

Definition at line 253 of file dense_matrix.C.

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

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

255 {
256  // Make sure the input sizes are compatible
257  libmesh_assert_equal_to (this->m(), arg.size());
258 
259  // Resize and clear dest.
260  // Note: DenseVector::resize() also zeros the vector.
261  dest.resize(this->n());
262 
263  // Short-circuit if the matrix is empty
264  if(this->m() == 0)
265  return;
266 
267  if (this->use_blas_lapack)
268  {
269  this->_matvec_blas(1., 0., dest, arg, /*trans=*/true);
270  }
271  else
272  {
273  const unsigned int n_rows = this->m();
274  const unsigned int n_cols = this->n();
275 
276  // WORKS
277  // for(unsigned int j=0; j<n_cols; j++)
278  // for(unsigned int i=0; i<n_rows; i++)
279  // dest(j) += (*this)(i,j)*arg(i);
280 
281  // ALSO WORKS, (i,j) just swapped
282  for(unsigned int i=0; i<n_cols; i++)
283  for(unsigned int j=0; j<n_rows; j++)
284  dest(i) += (*this)(j,i)*arg(j);
285  }
286 }
template<typename T >
void libMesh::DenseMatrix< T >::zero ( )
inlinevirtual

Set every element in the matrix to 0.

Implements libMesh::DenseMatrixBase< T >.

Definition at line 640 of file dense_matrix.h.

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

641 {
643 
644  // Just doing this ifdef to be completely safe
645 #ifndef LIBMESH_USE_COMPLEX_NUMBERS
646  if(_val.size())
647  std::memset(&_val[0], 0, sizeof(T) * _val.size());
648 #else
649  std::fill (_val.begin(), _val.end(), 0.);
650 #endif
651 }

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

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

template<typename T>
unsigned int libMesh::DenseMatrixBase< T >::_m
protectedinherited
template<typename T>
unsigned int libMesh::DenseMatrixBase< T >::_n
protectedinherited
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 488 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 368 of file dense_matrix.h.


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

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:59 UTC

Hosted By:
SourceForge.net Logo