#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)
template<typename T2 >
void	left_multiply (const DenseMatrixBase< T2 > &M2)
virtual void	right_multiply (const DenseMatrixBase< T > &M2)
template<typename T2 >
void	right_multiply (const DenseMatrixBase< T2 > &M2)
void	vector_mult (DenseVector< T > &dest, const DenseVector< T > &arg) const
template<typename T2 >
void	vector_mult (DenseVector< typename CompareTypes< T, T2 >::supertype > &dest, const DenseVector< T2 > &arg) const
void	vector_mult_transpose (DenseVector< T > &dest, const DenseVector< T > &arg) const
template<typename T2 >
void	vector_mult_transpose (DenseVector< typename CompareTypes< T, T2 >::supertype > &dest, const DenseVector< T2 > &arg) const
void	vector_mult_add (DenseVector< T > &dest, const T factor, const DenseVector< T > &arg) const
template<typename T2 , typename T3 >
void	vector_mult_add (DenseVector< typename CompareTypes< T, typename CompareTypes< T2, T3 >::supertype >::supertype > &dest, const T2 factor, const DenseVector< T3 > &arg) const
void	get_principal_submatrix (unsigned int sub_m, unsigned int sub_n, DenseMatrix< T > &dest) const
void	get_principal_submatrix (unsigned int sub_m, DenseMatrix< T > &dest) const
DenseMatrix< T > &	operator= (const DenseMatrix< T > &other_matrix)
template<typename T2 >
DenseMatrix< T > &	operator= (const DenseMatrix< T2 > &other_matrix)
void	swap (DenseMatrix< T > &other_matrix)
void	resize (const unsigned int new_m, const unsigned int new_n)
void	scale (const T factor)
void	scale_column (const unsigned int col, const T factor)
DenseMatrix< T > &	operator*= (const T factor)
template<typename T2 , typename T3 >
boostcopy::enable_if_c < ScalarTraits< T2 >::value, void >::type	add (const T2 factor, const DenseMatrix< T3 > &mat)
bool	operator== (const DenseMatrix< T > &mat) const
bool	operator!= (const DenseMatrix< T > &mat) const
DenseMatrix< T > &	operator+= (const DenseMatrix< T > &mat)
DenseMatrix< T > &	operator-= (const DenseMatrix< T > &mat)
Real	min () const
Real	max () const
Real	l1_norm () const
Real	linfty_norm () const
void	left_multiply_transpose (const DenseMatrix< T > &A)
template<typename T2 >
void	left_multiply_transpose (const DenseMatrix< T2 > &A)
void	right_multiply_transpose (const DenseMatrix< T > &A)
template<typename T2 >
void	right_multiply_transpose (const DenseMatrix< T2 > &A)
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)
void	evd (DenseVector< T > &lambda_real, DenseVector< T > &lambda_imag)
T	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	_evd_lapack (DenseVector< T > &lambda_real, DenseVector< T > &lambda_imag)
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 50 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 493 of file dense_matrix.h.

                           {
    LEFT_MULTIPLY = 0,
    RIGHT_MULTIPLY,
    LEFT_MULTIPLY_TRANSPOSE,
    RIGHT_MULTIPLY_TRANSPOSE
  };

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

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

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

  : DenseMatrixBase<T>(new_m,new_n),
#if defined(LIBMESH_HAVE_PETSC) && defined(LIBMESH_USE_REAL_NUMBERS) && defined(LIBMESH_DEFAULT_DOUBLE_PRECISION)
    use_blas_lapack(true),
#else
    use_blas_lapack(false),
#endif
    _val(),
    _decomposition_type(NONE),
    _pivots()
{
  this->resize(new_m,new_n);
}

template<typename T>

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

Copy-constructor. Destructor. Empty.

Definition at line 68 of file dense_matrix.h.

{}

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

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

{
  // Shorthand notation for number of rows and columns.
  const unsigned int
    n_rows = this->m(),
    n_cols = this->n();

  // Just to be really sure...
  libmesh_assert_equal_to (n_rows, n_cols);

  // A convenient reference to *this
  const DenseMatrix<T>& A = *this;

  // Now compute the solution to Ax =b using the factorization.
  x.resize(n_rows);

  // Solve for Ly=b
  for (unsigned int i=0; i<n_cols; ++i)
    {
      T2 temp = b(i);

      for (unsigned int k=0; k<i; ++k)
        temp -= A(i,k)*x(k);

      x(i) = temp / A(i,i);
    }

  // Solve for L^T x = y
  for (unsigned int i=0; i<n_cols; ++i)
    {
      const unsigned int ib = (n_cols-1)-i;

      for (unsigned int k=(ib+1); k<n_cols; ++k)
        x(ib) -= A(k,ib) * x(k);

      x(ib) /= A(ib,ib);
    }
}

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

{
  // If we called this function, there better not be any
  // previous decomposition of the matrix.
  libmesh_assert_equal_to (this->_decomposition_type, NONE);

  // Shorthand notation for number of rows and columns.
  const unsigned int
    n_rows = this->m(),
    n_cols = this->n();

  // Just to be really sure...
  libmesh_assert_equal_to (n_rows, n_cols);

  // A convenient reference to *this
  DenseMatrix<T>& A = *this;

  for (unsigned int i=0; i<n_rows; ++i)
    {
      for (unsigned int j=i; j<n_cols; ++j)
        {
          for (unsigned int k=0; k<i; ++k)
            A(i,j) -= A(i,k) * A(j,k);

          if (i == j)
            {
#ifndef LIBMESH_USE_COMPLEX_NUMBERS
              if (A(i,j) <= 0.0)
                libmesh_error_msg("Error! Can only use Cholesky decomposition with symmetric positive definite matrices.");
#endif

              A(i,i) = std::sqrt(A(i,j));
            }
          else
            A(j,i) = A(i,j) / A(i,i);
        }
    }

  // Set the flag for CHOLESKY decomposition
  this->_decomposition_type = CHOLESKY;
}

template<typename T>

template void libMesh::DenseMatrix< T >::_evd_lapack	(	DenseVector< T > &	lambda_real,
		DenseVector< T > &	lambda_imag
	)		`[private]`

Computes the eigenvalues of the matrix using the Lapack routine "DGEEV". [ Implementation in dense_matrix_blas_lapack.C ]

Definition at line 518 of file dense_matrix_blas_lapack.C.

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

{
  // The calling sequence for dgeev is:
  // DGEEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, WR, WI, VL, LDVL, VR,
  //         LDVR, ILO, IHI, SCALE, ABNRM, RCONDE, RCONDV, WORK, LWORK, IWORK, INFO )


  //  BALANC  (input) CHARACTER*1
  //          Indicates how the input matrix should be diagonally scaled
  //          and/or permuted to improve the conditioning of its
  //          eigenvalues.
  //          = 'N': Do not diagonally scale or permute;
  char BALANC = 'N';

  //  JOBVL   (input) CHARACTER*1
  //          = 'N': left eigenvectors of A are not computed;
  //          = 'V': left eigenvectors of A are computed.
  char JOBVL = 'N';

  //  JOBVR   (input) CHARACTER*1
  //          = 'N': right eigenvectors of A are not computed;
  //          = 'V': right eigenvectors of A are computed.
  char JOBVR = 'N';

  //  SENSE   (input) CHARACTER*1
  //          Determines which reciprocal condition numbers are computed.
  //          = 'N': None are computed;
  //          = 'E': Computed for eigenvalues only;
  //          = 'V': Computed for right eigenvectors only;
  //          = 'B': Computed for eigenvalues and right eigenvectors.
  char SENSE = 'N';

  //    N       (input) int*
  //            The number of rows/cols of the matrix A.  N >= 0.
  libmesh_assert( this->m() == this->n() );
  int N = this->m();

  //  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
  //          On entry, the N-by-N matrix A.
  //          On exit, A has been overwritten.
  // Here, we pass &(_val[0]).

  //    LDA     (input) int*
  //            The leading dimension of the array A.  LDA >= max(1,N).
  int LDA = N;

  //  WR      (output) DOUBLE PRECISION array, dimension (N)
  //  WI      (output) DOUBLE PRECISION array, dimension (N)
  //          WR and WI contain the real and imaginary parts,
  //          respectively, of the computed eigenvalues.  Complex
  //          conjugate pairs of eigenvalues appear consecutively
  //          with the eigenvalue having the positive imaginary part
  //          first.
  lambda_real.resize(N);
  lambda_imag.resize(N);

  //  VL      (output) DOUBLE PRECISION array, dimension (LDVL,N)
  //          If JOBVL = 'V', the left eigenvectors u(j) are stored one
  //          after another in the columns of VL, in the same order
  //          as their eigenvalues.
  //          If JOBVL = 'N', VL is not referenced.
  //          If the j-th eigenvalue is real, then u(j) = VL(:,j),
  //          the j-th column of VL.
  //          If the j-th and (j+1)-st eigenvalues form a complex
  //          conjugate pair, then u(j) = VL(:,j) + i*VL(:,j+1) and
  //          u(j+1) = VL(:,j) - i*VL(:,j+1).
  // Just set to NULL here.

  //  LDVL    (input) INTEGER
  //          The leading dimension of the array VL.  LDVL >= 1; if
  //          JOBVL = 'V', LDVL >= N.
  int LDVL = 1;

  //  VR      (output) DOUBLE PRECISION array, dimension (LDVR,N)
  //          If JOBVR = 'V', the right eigenvectors v(j) are stored one
  //          after another in the columns of VR, in the same order
  //          as their eigenvalues.
  //          If JOBVR = 'N', VR is not referenced.
  //          If the j-th eigenvalue is real, then v(j) = VR(:,j),
  //          the j-th column of VR.
  //          If the j-th and (j+1)-st eigenvalues form a complex
  //          conjugate pair, then v(j) = VR(:,j) + i*VR(:,j+1) and
  //          v(j+1) = VR(:,j) - i*VR(:,j+1).
  // Just set to NULL here.

  //  LDVR    (input) INTEGER
  //          The leading dimension of the array VR.  LDVR >= 1; if
  //          JOBVR = 'V', LDVR >= N.
  int LDVR = 1;

  // Outputs (unused)
  int ILO = 0;
  int IHI = 0;
  std::vector<T> SCALE(N);
  T ABNRM;
  std::vector<T> RCONDE(N);
  std::vector<T> RCONDV(N);

  //  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
  //          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  //
  //  LWORK   (input) INTEGER
  //          The dimension of the array WORK.
  int LWORK = 3*N;
  std::vector<T> WORK( LWORK );

  //  IWORK   (workspace) INTEGER array, dimension (2*N-2)
  //          If SENSE = 'N' or 'E', not referenced.
  // Just set to NULL


  //  INFO    (output) INTEGER
  //          = 0:  successful exit
  //          < 0:  if INFO = -i, the i-th argument had an illegal value.
  //          > 0:  if INFO = i, the QR algorithm failed to compute all the
  //                eigenvalues, and no eigenvectors or condition numbers
  //                have been computed; elements 1:ILO-1 and i+1:N of WR
  //                and WI contain eigenvalues which have converged.
  int INFO = 0;

  // Get references to raw data
  std::vector<T>& lambda_real_val = lambda_real.get_values();
  std::vector<T>& lambda_imag_val = lambda_imag.get_values();

  // Ready to call the actual factorization routine through SLEPc's interface
  LAPACKgeevx_( &BALANC, &JOBVL, &JOBVR, &SENSE, &N, &(_val[0]), &LDA, &lambda_real_val[0],
                &lambda_imag_val[0], NULL, &LDVL, NULL, &LDVR, &ILO, &IHI, &SCALE[0], &ABNRM,
                &RCONDE[0], &RCONDV[0], &WORK[0], &LWORK, NULL, &INFO );

  // Check return value for errors
  if (INFO != 0)
    libmesh_error_msg("INFO=" << INFO << ", Error during Lapack eigenvalue calculation!");
}

template<typename T>

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

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

Definition at line 646 of file dense_matrix.C.

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

{
  const unsigned int
    n_cols = this->n();

  libmesh_assert_equal_to (this->m(), n_cols);
  libmesh_assert_equal_to (this->m(), b.size());

  x.resize (n_cols);

  // A convenient reference to *this
  const DenseMatrix<T>& A = *this;

  // Temporary vector storage.  We use this instead of
  // modifying the RHS.
  DenseVector<T> z = b;

  // Lower-triangular "top to bottom" solve step, taking into account pivots
  for (unsigned int i=0; i<n_cols; ++i)
    {
      // Swap
      if (_pivots[i] != static_cast<int>(i))
        std::swap( z(i), z(_pivots[i]) );

      x(i) = z(i);

      for (unsigned int j=0; j<i; ++j)
        x(i) -= A(i,j)*x(j);

      x(i) /= A(i,i);
    }

  // Upper-triangular "bottom to top" solve step
  const unsigned int last_row = n_cols-1;

  for (int i=last_row; i>=0; --i)
    {
      for (int j=i+1; j<static_cast<int>(n_cols); ++j)
        x(i) -= A(i,j)*x(j);
    }
}

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 669 of file dense_matrix_blas_lapack.C.

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

{
  // The calling sequence for getrs is:
  // dgetrs(TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO)

  //    trans   (input) char*
  //            'n' for no tranpose, 't' for transpose
  char TRANS[] = "t";

  //    N       (input) int*
  //            The order of the matrix A.  N >= 0.
  int N = this->m();


  //    NRHS    (input) int*
  //            The number of right hand sides, i.e., the number of columns
  //            of the matrix B.  NRHS >= 0.
  int NRHS = 1;

  //    A       (input) DOUBLE PRECISION array, dimension (LDA,N)
  //            The factors L and U from the factorization A = P*L*U
  //            as computed by dgetrf.
  // Here, we pass &(_val[0])

  //    LDA     (input) int*
  //            The leading dimension of the array A.  LDA >= max(1,N).
  int LDA = N;

  //    ipiv    (input) int array, dimension (N)
  //            The pivot indices from DGETRF; for 1<=i<=N, row i of the
  //            matrix was interchanged with row IPIV(i).
  // Here, we pass &(_pivots[0]) which was computed in _lu_decompose_lapack

  //    B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
  //            On entry, the right hand side matrix B.
  //            On exit, the solution matrix X.
  // Here, we pass a copy of the rhs vector's data array in x, so that the
  // passed right-hand side b is unmodified.  I don't see a way around this
  // copy if we want to maintain an unmodified rhs in LibMesh.
  x = b;
  std::vector<T>& x_vec = x.get_values();

  // We can avoid the copy if we don't care about overwriting the RHS: just
  // pass b to the Lapack routine and then swap with x before exiting
  // std::vector<T>& x_vec = b.get_values();

  //    LDB     (input) int*
  //            The leading dimension of the array B.  LDB >= max(1,N).
  int LDB = N;

  //    INFO    (output) int*
  //            = 0:  successful exit
  //            < 0:  if INFO = -i, the i-th argument had an illegal value
  int INFO = 0;

  // Finally, ready to call the Lapack getrs function
  LAPACKgetrs_(TRANS, &N, &NRHS, &(_val[0]), &LDA, &(_pivots[0]), &(x_vec[0]), &LDB, &INFO);

  // Check return value for errors
  if (INFO != 0)
    libmesh_error_msg("INFO=" << INFO << ", Error during Lapack LU solve!");

  // Don't do this if you already made a copy of b above
  // Swap b and x.  The solution will then be in x, and whatever was originally
  // in x, maybe garbage, maybe nothing, will be in b.
  // FIXME: Rewrite the LU and Cholesky solves to just take one input, and overwrite
  // the input.  This *should* make user code simpler, as they don't have to create
  // an extra vector just to pass it in to the solve function!
  // b.swap(x);
}

template<typename T >

void libMesh::DenseMatrix< T >::_lu_decompose ( ) [private]

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

Definition at line 697 of file dense_matrix.C.

References std::abs(), libMesh::Real, libMesh::DenseMatrix< T >::resize(), swap(), and libMesh::zero.

{
  // 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_error_msg("Matrix A is singular!");

      // 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 ( ) [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 211 of file dense_matrix_blas_lapack.C.

References std::min().

{
  // If this function was called, there better not be any
  // previous decomposition of the matrix.
  libmesh_assert_equal_to (this->_decomposition_type, NONE);

  // The calling sequence for dgetrf is:
  // dgetrf(M, N, A, lda, ipiv, info)

  //    M       (input) int*
  //            The number of rows of the matrix A.  M >= 0.
  // In C/C++, pass the number of *cols* of A
  int M = this->n();

  //    N       (input) int*
  //            The number of columns of the matrix A.  N >= 0.
  // In C/C++, pass the number of *rows* of A
  int N = this->m();

  //    A (input/output) double precision array, dimension (LDA,N)
  //      On entry, the M-by-N matrix to be factored.
  //      On exit, the factors L and U from the factorization
  //      A = P*L*U; the unit diagonal elements of L are not stored.
  // Here, we pass &(_val[0]).

  //    LDA     (input) int*
  //            The leading dimension of the array A.  LDA >= max(1,M).
  int LDA = M;

  //    ipiv    (output) integer array, dimension (min(m,n))
  //            The pivot indices; for 1 <= i <= min(m,n), row i of the
  //            matrix was interchanged with row IPIV(i).
  // Here, we pass &(_pivots[0]), a private class member used to store pivots
  this->_pivots.resize( std::min(M,N) );

  //    info    (output) int*
  //            = 0:  successful exit
  //            < 0:  if INFO = -i, the i-th argument had an illegal value
  //            > 0:  if INFO = i, U(i,i) is exactly zero. The factorization
  //                  has been completed, but the factor U is exactly
  //                  singular, and division by zero will occur if it is used
  //                  to solve a system of equations.
  int INFO = 0;

  // Ready to call the actual factorization routine through PETSc's interface
  LAPACKgetrf_(&M, &N, &(this->_val[0]), &LDA, &(_pivots[0]), &INFO);

  // Check return value for errors
  if (INFO != 0)
    libmesh_error_msg("INFO=" << INFO << ", Error during Lapack LU factorization!");

  // Set the flag for LU decomposition
  this->_decomposition_type = LU_BLAS_LAPACK;
}

template<typename T>

template void libMesh::DenseMatrix< T >::_matvec_blas	(	T	alpha,
		T	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 759 of file dense_matrix_blas_lapack.C.

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

{
  // Ensure that dest and arg sizes are compatible
  if (!trans)
    {
      // dest  ~ A     * arg
      // (mx1)   (mxn) * (nx1)
      if ((dest.size() != this->m()) || (arg.size() != this->n()))
        libmesh_error_msg("Improper input argument sizes!");
    }

  else // trans == true
    {
      // Ensure that dest and arg are proper size
      // dest  ~ A^T   * arg
      // (nx1)   (nxm) * (mx1)
      if ((dest.size() != this->n()) || (arg.size() != this->m()))
        libmesh_error_msg("Improper input argument sizes!");
    }

  // Calling sequence for dgemv:
  //
  // dgemv(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)

  //   TRANS  - CHARACTER*1, 't' for transpose, 'n' for non-transpose multiply
  // We store everything in row-major order, so pass the transpose flag for
  // non-transposed matvecs and the 'n' flag for transposed matvecs
  char TRANS[] = "t";
  if (trans)
    TRANS[0] = 'n';

  //   M      - INTEGER.
  //            On entry, M specifies the number of rows of the matrix A.
  // In C/C++, pass the number of *cols* of A
  int M = this->n();

  //   N      - INTEGER.
  //            On entry, N specifies the number of columns of the matrix A.
  // In C/C++, pass the number of *rows* of A
  int N = this->m();

  //   ALPHA  - DOUBLE PRECISION.
  // The scalar constant passed to this function

  //   A      - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
  //            Before entry, the leading m by n part of the array A must
  //            contain the matrix of coefficients.
  // The matrix, *this.  Note that _matvec_blas is called from
  // a const function, vector_mult(), and so we have made this function const
  // as well.  Since BLAS knows nothing about const, we have to cast it away
  // now.
  DenseMatrix<T>& a_ref = const_cast< DenseMatrix<T>& > ( *this );
  std::vector<T>& a = a_ref.get_values();

  //   LDA    - INTEGER.
  //            On entry, LDA specifies the first dimension of A as declared
  //            in the calling (sub) program. LDA must be at least
  //            max( 1, m ).
  int LDA = M;

  //   X      - DOUBLE PRECISION array of DIMENSION at least
  //            ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
  //            and at least
  //            ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
  //            Before entry, the incremented array X must contain the
  //            vector x.
  // Here, we must cast away the const-ness of "arg" since BLAS knows
  // nothing about const
  DenseVector<T>& x_ref = const_cast< DenseVector<T>& > ( arg );
  std::vector<T>& x = x_ref.get_values();

  //   INCX   - INTEGER.
  //            On entry, INCX specifies the increment for the elements of
  //            X. INCX must not be zero.
  int INCX = 1;

  //   BETA   - DOUBLE PRECISION.
  //            On entry, BETA specifies the scalar beta. When BETA is
  //            supplied as zero then Y need not be set on input.
  // The second scalar constant passed to this function

  //   Y      - DOUBLE PRECISION array of DIMENSION at least
  //            ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
  //            and at least
  //            ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
  //            Before entry with BETA non-zero, the incremented array Y
  //            must contain the vector y. On exit, Y is overwritten by the
  //            updated vector y.
  // The input vector "dest"
  std::vector<T>& y = dest.get_values();

  //   INCY   - INTEGER.
  //            On entry, INCY specifies the increment for the elements of
  //            Y. INCY must not be zero.
  int INCY = 1;

  // Finally, ready to call the BLAS function
  BLASgemv_(TRANS, &M, &N, &alpha, &(a[0]), &LDA, &(x[0]), &INCX, &beta, &(y[0]), &INCY);
}

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 44 of file dense_matrix_blas_lapack.C.

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

{
  int result_size = 0;

  // For each case, determine the size of the final result make sure
  // that the inner dimensions match
  switch (flag)
    {
    case LEFT_MULTIPLY:
      {
        result_size = other.m() * this->n();
        if (other.n() == this->m())
          break;
      }
    case RIGHT_MULTIPLY:
      {
        result_size = other.n() * this->m();
        if (other.m() == this->n())
          break;
      }
    case LEFT_MULTIPLY_TRANSPOSE:
      {
        result_size = other.n() * this->n();
        if (other.m() == this->m())
          break;
      }
    case RIGHT_MULTIPLY_TRANSPOSE:
      {
        result_size = other.m() * this->m();
        if (other.n() == this->n())
          break;
      }
    default:
      libmesh_error_msg("Unknown flag selected or matrices are incompatible for multiplication.");
    }

  // For this to work, the passed arg. must actually be a DenseMatrix<T>
  const DenseMatrix<T>* const_that = cast_ptr< const DenseMatrix<T>* >(&other);

  // Also, although 'that' is logically const in this BLAS routine,
  // the PETSc BLAS interface does not specify that any of the inputs are
  // const.  To use it, I must cast away const-ness.
  DenseMatrix<T>* that = const_cast< DenseMatrix<T>* > (const_that);

  // Initialize A, B pointers for LEFT_MULTIPLY* cases
  DenseMatrix<T>
    *A = this,
    *B = that;

  // For RIGHT_MULTIPLY* cases, swap the meaning of A and B.
  // Here is a full table of combinations we can pass to BLASgemm, and what the answer is when finished:
  // pass A B   -> (Fortran) -> A^T B^T -> (C++) -> (A^T B^T)^T -> (identity) -> B A   "lt multiply"
  // pass B A   -> (Fortran) -> B^T A^T -> (C++) -> (B^T A^T)^T -> (identity) -> A B   "rt multiply"
  // pass A B^T -> (Fortran) -> A^T B   -> (C++) -> (A^T B)^T   -> (identity) -> B^T A "lt multiply t"
  // pass B^T A -> (Fortran) -> B A^T   -> (C++) -> (B A^T)^T   -> (identity) -> A B^T "rt multiply t"
  if (flag==RIGHT_MULTIPLY || flag==RIGHT_MULTIPLY_TRANSPOSE)
    std::swap(A,B);

  // transa, transb values to pass to blas
  char
    transa[] = "n",
    transb[] = "n";

  // Integer values to pass to BLAS:
  //
  // M
  // In Fortran, the number of rows of op(A),
  // In the BLAS documentation, typically known as 'M'.
  //
  // In C/C++, we set:
  // M = n_cols(A) if (transa='n')
  //     n_rows(A) if (transa='t')
  int M = static_cast<int>( A->n() );

  // N
  // In Fortran, the number of cols of op(B), and also the number of cols of C.
  // In the BLAS documentation, typically known as 'N'.
  //
  // In C/C++, we set:
  // N = n_rows(B) if (transb='n')
  //     n_cols(B) if (transb='t')
  int N = static_cast<int>( B->m() );

  // K
  // In Fortran, the number of cols of op(A), and also
  // the number of rows of op(B). In the BLAS documentation,
  // typically known as 'K'.
  //
  // In C/C++, we set:
  // K = n_rows(A) if (transa='n')
  //     n_cols(A) if (transa='t')
  int K = static_cast<int>( A->m() );

  // LDA (leading dimension of A). In our cases,
  // LDA is always the number of columns of A.
  int LDA = static_cast<int>( A->n() );

  // LDB (leading dimension of B).  In our cases,
  // LDB is always the number of columns of B.
  int LDB = static_cast<int>( B->n() );

  if (flag == LEFT_MULTIPLY_TRANSPOSE)
    {
      transb[0] = 't';
      N = static_cast<int>( B->n() );
    }

  else if (flag == RIGHT_MULTIPLY_TRANSPOSE)
    {
      transa[0] = 't';
      std::swap(M,K);
    }

  // LDC (leading dimension of C).  LDC is the
  // number of columns in the solution matrix.
  int LDC = M;

  // Scalar values to pass to BLAS
  //
  // scalar multiplying the whole product AB
  T alpha = 1.;

  // scalar multiplying C, which is the original matrix.
  T beta  = 0.;

  // Storage for the result
  std::vector<T> result (result_size);

  // Finally ready to call the BLAS
  BLASgemm_(transa, transb, &M, &N, &K, &alpha, &(A->_val[0]), &LDA, &(B->_val[0]), &LDB, &beta, &result[0], &LDC);

  // Update the relevant dimension for this matrix.
  switch (flag)
    {
    case LEFT_MULTIPLY:            { this->_m = other.m(); break; }
    case RIGHT_MULTIPLY:           { this->_n = other.n(); break; }
    case LEFT_MULTIPLY_TRANSPOSE:  { this->_m = other.n(); break; }
    case RIGHT_MULTIPLY_TRANSPOSE: { this->_n = other.m(); break; }
    default:
      libmesh_error_msg("Unknown flag selected.");
    }

  // Swap my data vector with the result
  this->_val.swap(result);
}

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 390 of file dense_matrix_blas_lapack.C.

{

  //    M       (input) int*
  //            The number of rows of the matrix A.  M >= 0.
  // In C/C++, pass the number of *cols* of A
  int M = this->n();

  //    N       (input) int*
  //            The number of columns of the matrix A.  N >= 0.
  // In C/C++, pass the number of *rows* of A
  int N = this->m();

  int min_MN = (M < N) ? M : N;
  int max_MN = (M > N) ? M : N;

  //  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
  //          On entry, the M-by-N matrix A.
  //          On exit,
  //          if JOBU = 'O',  A is overwritten with the first min(m,n)
  //                          columns of U (the left singular vectors,
  //                          stored columnwise);
  //          if JOBVT = 'O', A is overwritten with the first min(m,n)
  //                          rows of V**T (the right singular vectors,
  //                          stored rowwise);
  //          if JOBU .ne. 'O' and JOBVT .ne. 'O', the contents of A
  //                          are destroyed.
  // Here, we pass &(_val[0]).

  //    LDA     (input) int*
  //            The leading dimension of the array A.  LDA >= max(1,M).
  int LDA = M;

  //  S       (output) DOUBLE PRECISION array, dimension (min(M,N))
  //          The singular values of A, sorted so that S(i) >= S(i+1).
  sigma_val.resize( min_MN );

  //  LDU     (input) INTEGER
  //          The leading dimension of the array U.  LDU >= 1; if
  //          JOBU = 'S' or 'A', LDU >= M.
  int LDU = M;

  //  U       (output) DOUBLE PRECISION array, dimension (LDU,UCOL)
  //          (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'.
  //          If JOBU = 'A', U contains the M-by-M orthogonal matrix U;
  //          if JOBU = 'S', U contains the first min(m,n) columns of U
  //          (the left singular vectors, stored columnwise);
  //          if JOBU = 'N' or 'O', U is not referenced.
  U_val.resize( LDU*M );

  //  LDVT    (input) INTEGER
  //          The leading dimension of the array VT.  LDVT >= 1; if
  //          JOBVT = 'A', LDVT >= N; if JOBVT = 'S', LDVT >= min(M,N).
  int LDVT = N;

  //  VT      (output) DOUBLE PRECISION array, dimension (LDVT,N)
  //          If JOBVT = 'A', VT contains the N-by-N orthogonal matrix
  //          V**T;
  //          if JOBVT = 'S', VT contains the first min(m,n) rows of
  //          V**T (the right singular vectors, stored rowwise);
  //          if JOBVT = 'N' or 'O', VT is not referenced.
  VT_val.resize( LDVT*N );

  //  LWORK   (input) INTEGER
  //          The dimension of the array WORK.
  //          LWORK >= MAX(1,3*MIN(M,N)+MAX(M,N),5*MIN(M,N)).
  //          For good performance, LWORK should generally be larger.
  //
  //          If LWORK = -1, then a workspace query is assumed; the routine
  //          only calculates the optimal size of the WORK array, returns
  //          this value as the first entry of the WORK array, and no error
  //          message related to LWORK is issued by XERBLA.
  int larger = (3*min_MN+max_MN > 5*min_MN) ? 3*min_MN+max_MN : 5*min_MN;
  int LWORK  = (larger > 1) ? larger : 1;


  //  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
  //          On exit, if INFO = 0, WORK(1) returns the optimal LWORK;
  //          if INFO > 0, WORK(2:MIN(M,N)) contains the unconverged
  //          superdiagonal elements of an upper bidiagonal matrix B
  //          whose diagonal is in S (not necessarily sorted). B
  //          satisfies A = U * B * VT, so it has the same singular values
  //          as A, and singular vectors related by U and VT.
  std::vector<T> WORK( LWORK );

  //  INFO    (output) INTEGER
  //          = 0:  successful exit.
  //          < 0:  if INFO = -i, the i-th argument had an illegal value.
  //          > 0:  if DBDSQR did not converge, INFO specifies how many
  //                superdiagonals of an intermediate bidiagonal form B
  //                did not converge to zero. See the description of WORK
  //                above for details.
  int INFO = 0;

  // Ready to call the actual factorization routine through PETSc's interface
  LAPACKgesvd_(&JOBU, &JOBVT, &M, &N, &(_val[0]), &LDA, &(sigma_val[0]), &(U_val[0]),
               &LDU, &(VT_val[0]), &LDVT, &(WORK[0]), &LWORK, &INFO);

  // Check return value for errors
  if (INFO != 0)
    libmesh_error_msg("INFO=" << INFO << ", Error during Lapack SVD calculation!");
}

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 279 of file dense_matrix_blas_lapack.C.

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

{
  // The calling sequence for dgetrf is:
  // DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, INFO )


  //  JOBU    (input) CHARACTER*1
  //          Specifies options for computing all or part of the matrix U:
  //          = 'A':  all M columns of U are returned in array U:
  //          = 'S':  the first min(m,n) columns of U (the left singular
  //                  vectors) are returned in the array U;
  //          = 'O':  the first min(m,n) columns of U (the left singular
  //                  vectors) are overwritten on the array A;
  //          = 'N':  no columns of U (no left singular vectors) are
  //                  computed.
  char JOBU = 'N';

  //  JOBVT   (input) CHARACTER*1
  //          Specifies options for computing all or part of the matrix
  //          V**T:
  //          = 'A':  all N rows of V**T are returned in the array VT;
  //          = 'S':  the first min(m,n) rows of V**T (the right singular
  //                  vectors) are returned in the array VT;
  //          = 'O':  the first min(m,n) rows of V**T (the right singular
  //                  vectors) are overwritten on the array A;
  //          = 'N':  no rows of V**T (no right singular vectors) are
  //                  computed.
  char JOBVT = 'N';

  std::vector<T> sigma_val;
  std::vector<T> U_val;
  std::vector<T> VT_val;

  _svd_helper(JOBU, JOBVT, sigma_val, U_val, VT_val);

  // Load the singular values into sigma, ignore U_val and VT_val
  const unsigned int n_sigma_vals =
    cast_int<unsigned int>(sigma_val.size());
  sigma.resize(n_sigma_vals);
  for(unsigned int i=0; i<n_sigma_vals; i++)
    sigma(i) = sigma_val[i];

}

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

{
  // The calling sequence for dgetrf is:
  // DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, INFO )


  //  JOBU    (input) CHARACTER*1
  //          Specifies options for computing all or part of the matrix U:
  //          = 'A':  all M columns of U are returned in array U:
  //          = 'S':  the first min(m,n) columns of U (the left singular
  //                  vectors) are returned in the array U;
  //          = 'O':  the first min(m,n) columns of U (the left singular
  //                  vectors) are overwritten on the array A;
  //          = 'N':  no columns of U (no left singular vectors) are
  //                  computed.
  char JOBU = 'S';

  //  JOBVT   (input) CHARACTER*1
  //          Specifies options for computing all or part of the matrix
  //          V**T:
  //          = 'A':  all N rows of V**T are returned in the array VT;
  //          = 'S':  the first min(m,n) rows of V**T (the right singular
  //                  vectors) are returned in the array VT;
  //          = 'O':  the first min(m,n) rows of V**T (the right singular
  //                  vectors) are overwritten on the array A;
  //          = 'N':  no rows of V**T (no right singular vectors) are
  //                  computed.
  char JOBVT = 'S';

  std::vector<T> sigma_val;
  std::vector<T> U_val;
  std::vector<T> VT_val;

  _svd_helper(JOBU, JOBVT, sigma_val, U_val, VT_val);

  // Load the singular values into sigma, ignore U_val and VT_val
  const unsigned int n_sigma_vals =
    cast_int<unsigned int>(sigma_val.size());
  sigma.resize(n_sigma_vals);
  for(unsigned int i=0; i<n_sigma_vals; i++)
    sigma(i) = sigma_val[i];

  int M = this->n();
  int N = this->m();
  int min_MN = (M < N) ? M : N;
  U.resize(M,min_MN);
  for(unsigned int i=0; i<U.m(); i++)
    for(unsigned int j=0; j<U.n(); j++)
      {
        unsigned int index = i + j*U.n();  // Column major storage
        U(i,j) = U_val[index];
      }

  VT.resize(min_MN,N);
  for(unsigned int i=0; i<VT.m(); i++)
    for(unsigned int j=0; j<VT.n(); j++)
      {
        unsigned int index = i + j*U.n(); // Column major storage
        VT(i,j) = VT_val[index];
      }

}

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

{
  libmesh_assert_equal_to (this->m(), mat.m());
  libmesh_assert_equal_to (this->n(), mat.n());

  for (unsigned int j=0; j<this->n(); j++)
    for (unsigned int i=0; i<this->m(); i++)
      this->el(i,j) += factor*mat.el(i,j);
}

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

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

Referenced by libMesh::FEMSystem::assembly().

{
  libmesh_assert_equal_to (this->m(), mat.m());
  libmesh_assert_equal_to (this->n(), mat.n());

  for (unsigned int i=0; i<this->m(); i++)
    for (unsigned int j=0; j<this->n(); j++)
      (*this)(i,j) += factor * mat(i,j);
}

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

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

{
  // Check for a previous decomposition
  switch(this->_decomposition_type)
    {
    case NONE:
      {
        this->_cholesky_decompose ();
        break;
      }

    case CHOLESKY:
      {
        // Already factored, just need to call back_substitute.
        break;
      }

    default:
      libmesh_error_msg("Error! This matrix already has a different decomposition...");
    }

  // Perform back substitution
  this->_cholesky_back_substitute (b, x);
}

template<typename T>

void libMesh::DenseMatrixBase< T >::condense	(	const unsigned int	i,
		const unsigned int	j,
		const T	val,
		DenseVectorBase< T > &	rhs
	)		`[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::DenseVectorBase< T >::el(), and libMesh::DenseVectorBase< T >::size().

{
  libmesh_assert_equal_to (this->_m, rhs.size());
  libmesh_assert_equal_to (iv, jv);


  // move the known value into the RHS
  // and zero the column
  for (unsigned int i=0; i<this->m(); i++)
    {
      rhs.el(i) -= this->el(i,jv)*val;
      this->el(i,jv) = 0.;
    }

  // zero the row
  for (unsigned int j=0; j<this->n(); j++)
    this->el(iv,j) = 0.;

  this->el(iv,jv) = 1.;
  rhs.el(iv) = val;

}

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

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

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

References libMesh::Real.

{
  switch(this->_decomposition_type)
    {
    case NONE:
      {
        // First LU decompose the matrix.
        // Note that the lu_decompose routine will check to see if the
        // matrix is square so we don't worry about it.
        if (this->use_blas_lapack)
          this->_lu_decompose_lapack();
        else
          this->_lu_decompose ();
      }
    case LU:
    case LU_BLAS_LAPACK:
      {
        // Already decomposed, don't do anything
        break;
      }
    default:
      libmesh_error_msg("Error! Can't compute the determinant under the current decomposition.");
    }

  // A variable to keep track of the running product of diagonal terms.
  T determinant = 1.;

  // Loop over diagonal terms, computing the product.  In practice,
  // be careful because this value could easily become too large to
  // fit in a double or float.  To be safe, one should keep track of
  // the power (of 10) of the determinant in a separate variable
  // and maintain an order 1 value for the determinant itself.
  unsigned int n_interchanges = 0;
  for (unsigned int i=0; i<this->m(); i++)
    {
      if (this->_decomposition_type==LU)
        if (_pivots[i] != static_cast<int>(i))
          n_interchanges++;

      // Lapack pivots are 1-based!
      if (this->_decomposition_type==LU_BLAS_LAPACK)
        if (_pivots[i] != static_cast<int>(i+1))
          n_interchanges++;

      determinant *= (*this)(i,i);
    }

  // Compute sign of determinant, depends on number of row interchanges!
  // The sign should be (-1)^{n}, where n is the number of interchanges.
  Real sign = n_interchanges % 2 == 0 ? 1. : -1.;

  return sign*determinant;
}

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

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

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

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

template<typename T>

void libMesh::DenseMatrix< T >::evd	(	DenseVector< T > &	lambda_real,
		DenseVector< T > &	lambda_imag
	)

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

The implementation requires the LAPACKgeevx_ which is wrapped by SLEPc, and throws an error if called when SLEPc is not available.

Definition at line 792 of file dense_matrix.C.

{
  // We use the LAPACK eigenvalue problem implementation
  _evd_lapack(lambda_real, lambda_imag);
}

template<typename T>

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

Put the sub_m x sub_n principal submatrix into dest.

Definition at line 553 of file dense_matrix.C.

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

{
  libmesh_assert( (sub_m <= this->m()) && (sub_n <= this->n()) );

  dest.resize(sub_m, sub_n);
  for(unsigned int i=0; i<sub_m; i++)
    for(unsigned int j=0; j<sub_n; j++)
      dest(i,j) = (*this)(i,j);
}

template<typename T>

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

Put the sub_m x sub_m principal submatrix into dest.

Definition at line 568 of file dense_matrix.C.

{
  get_principal_submatrix(sub_m, sub_m, dest);
}

template<typename T>

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

Put the tranposed matrix into dest.

Definition at line 576 of file dense_matrix.C.

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

{
  dest.resize(this->n(), this->m());

  for (unsigned int i=0; i<dest.m(); i++)
    for (unsigned int j=0; j<dest.n(); j++)
      dest(i,j) = (*this)(j,i);
}

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

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

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

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

Referenced by libMesh::FEMSystem::assembly().

{
  libmesh_assert (this->_m);
  libmesh_assert (this->_n);

  Real columnsum = 0.;
  for (unsigned int i=0; i!=this->_m; i++)
    {
      columnsum += std::abs((*this)(i,0));
    }
  Real my_max = columnsum;
  for (unsigned int j=1; j!=this->_n; j++)
    {
      columnsum = 0.;
      for (unsigned int i=0; i!=this->_m; i++)
        {
          columnsum += std::abs((*this)(i,j));
        }
      my_max = (my_max > columnsum? my_max : columnsum);
    }
  return my_max;
}

template<typename T>

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

Left multipliess by the matrix M2.

Implements libMesh::DenseMatrixBase< T >.

Definition at line 37 of file dense_matrix.C.

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

{
  if (this->use_blas_lapack)
    this->_multiply_blas(M2, LEFT_MULTIPLY);
  else
    {
      // (*this) <- M2 * (*this)
      // Where:
      // (*this) = (m x n),
      // M2      = (m x p),
      // M3      = (p x n)

      // M3 is a copy of *this before it gets resize()d
      DenseMatrix<T> M3(*this);

      // Resize *this so that the result can fit
      this->resize (M2.m(), M3.n());

      // Call the multiply function in the base class
      this->multiply(*this, M2, M3);
    }
}

template<typename T >

template<typename T2 >

void libMesh::DenseMatrix< T >::left_multiply ( const DenseMatrixBase< T2 > & M2 )

Left multipliess by the matrix M2 of different type

Definition at line 64 of file dense_matrix.C.

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

{
  // (*this) <- M2 * (*this)
  // Where:
  // (*this) = (m x n),
  // M2      = (m x p),
  // M3      = (p x n)

  // M3 is a copy of *this before it gets resize()d
  DenseMatrix<T> M3(*this);

  // Resize *this so that the result can fit
  this->resize (M2.m(), M3.n());

  // Call the multiply function in the base class
  this->multiply(*this, M2, M3);
}

template<typename T>

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

Left multiplies by the transpose of the matrix A.

Definition at line 85 of file dense_matrix.C.

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

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

{
  if (this->use_blas_lapack)
    this->_multiply_blas(A, LEFT_MULTIPLY_TRANSPOSE);
  else
    {
      //Check to see if we are doing (A^T)*A
      if (this == &A)
        {
          //libmesh_here();
          DenseMatrix<T> B(*this);

          // Simple but inefficient way
          // return this->left_multiply_transpose(B);

          // More efficient, but more code way
          // If A is mxn, the result will be a square matrix of Size n x n.
          const unsigned int n_rows = A.m();
          const unsigned int n_cols = A.n();

          // resize() *this and also zero out all entries.
          this->resize(n_cols,n_cols);

          // Compute the lower-triangular part
          for (unsigned int i=0; i<n_cols; ++i)
            for (unsigned int j=0; j<=i; ++j)
              for (unsigned int k=0; k<n_rows; ++k) // inner products are over n_rows
                (*this)(i,j) += B(k,i)*B(k,j);

          // Copy lower-triangular part into upper-triangular part
          for (unsigned int i=0; i<n_cols; ++i)
            for (unsigned int j=i+1; j<n_cols; ++j)
              (*this)(i,j) = (*this)(j,i);
        }

      else
        {
          DenseMatrix<T> B(*this);

          this->resize (A.n(), B.n());

          libmesh_assert_equal_to (A.m(), B.m());
          libmesh_assert_equal_to (this->m(), A.n());
          libmesh_assert_equal_to (this->n(), B.n());

          const unsigned int m_s = A.n();
          const unsigned int p_s = A.m();
          const unsigned int n_s = this->n();

          // Do it this way because there is a
          // decent chance (at least for constraint matrices)
          // that A.transpose(i,k) = 0.
          for (unsigned int i=0; i<m_s; i++)
            for (unsigned int k=0; k<p_s; k++)
              if (A.transpose(i,k) != 0.)
                for (unsigned int j=0; j<n_s; j++)
                  (*this)(i,j) += A.transpose(i,k)*B(k,j);
        }
    }

}

template<typename T >

template<typename T2 >

void libMesh::DenseMatrix< T >::left_multiply_transpose ( const DenseMatrix< T2 > & A )

Left multiplies by the transpose of the matrix A of different type

Definition at line 151 of file dense_matrix.C.

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

{
  //Check to see if we are doing (A^T)*A
  if (this == &A)
    {
      //libmesh_here();
      DenseMatrix<T> B(*this);

      // Simple but inefficient way
      // return this->left_multiply_transpose(B);

      // More efficient, but more code way
      // If A is mxn, the result will be a square matrix of Size n x n.
      const unsigned int n_rows = A.m();
      const unsigned int n_cols = A.n();

      // resize() *this and also zero out all entries.
      this->resize(n_cols,n_cols);

      // Compute the lower-triangular part
      for (unsigned int i=0; i<n_cols; ++i)
        for (unsigned int j=0; j<=i; ++j)
          for (unsigned int k=0; k<n_rows; ++k) // inner products are over n_rows
            (*this)(i,j) += B(k,i)*B(k,j);

      // Copy lower-triangular part into upper-triangular part
      for (unsigned int i=0; i<n_cols; ++i)
        for (unsigned int j=i+1; j<n_cols; ++j)
          (*this)(i,j) = (*this)(j,i);
    }

  else
    {
      DenseMatrix<T> B(*this);

      this->resize (A.n(), B.n());

      libmesh_assert_equal_to (A.m(), B.m());
      libmesh_assert_equal_to (this->m(), A.n());
      libmesh_assert_equal_to (this->n(), B.n());

      const unsigned int m_s = A.n();
      const unsigned int p_s = A.m();
      const unsigned int n_s = this->n();

      // Do it this way because there is a
      // decent chance (at least for constraint matrices)
      // that A.transpose(i,k) = 0.
      for (unsigned int i=0; i<m_s; i++)
        for (unsigned int k=0; k<p_s; k++)
          if (A.transpose(i,k) != 0.)
            for (unsigned int j=0; j<n_s; j++)
              (*this)(i,j) += A.transpose(i,k)*B(k,j);
    }
}

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

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

{
  libmesh_assert (this->_m);
  libmesh_assert (this->_n);

  Real rowsum = 0.;
  for (unsigned int j=0; j!=this->_n; j++)
    {
      rowsum += std::abs((*this)(0,j));
    }
  Real my_max = rowsum;
  for (unsigned int i=1; i!=this->_m; i++)
    {
      rowsum = 0.;
      for (unsigned int j=0; j!=this->_n; j++)
        {
          rowsum += std::abs((*this)(i,j));
        }
      my_max = (my_max > rowsum? my_max : rowsum);
    }
  return my_max;
}

template<typename T>

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

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

Definition at line 589 of file dense_matrix.C.

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

{
  // Check to be sure that the matrix is square before attempting
  // an LU-solve.  In general, one can compute the LU factorization of
  // a non-square matrix, but:
  //
  // Overdetermined systems (m>n) have a solution only if enough of
  // the equations are linearly-dependent.
  //
  // Underdetermined systems (m<n) typically have infinitely many
  // solutions.
  //
  // We don't want to deal with either of these ambiguous cases here...
  libmesh_assert_equal_to (this->m(), this->n());

  switch(this->_decomposition_type)
    {
    case NONE:
      {
        if (this->use_blas_lapack)
          this->_lu_decompose_lapack();
        else
          this->_lu_decompose ();
        break;
      }

    case LU_BLAS_LAPACK:
      {
        // Already factored, just need to call back_substitute.
        if (this->use_blas_lapack)
          break;
      }

    case LU:
      {
        // Already factored, just need to call back_substitute.
        if ( !(this->use_blas_lapack) )
          break;
      }

    default:
      libmesh_error_msg("Error! This matrix already has a different decomposition...");
    }

  if (this->use_blas_lapack)
    this->_lu_back_substitute_lapack (b, x);
  else
    this->_lu_back_substitute (b, x);
}

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.

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

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

{
  libmesh_assert (this->_m);
  libmesh_assert (this->_n);
  Real my_max = libmesh_real((*this)(0,0));

  for (unsigned int i=0; i!=this->_m; i++)
    {
      for (unsigned int j=0; j!=this->_n; j++)
        {
          Real current = libmesh_real((*this)(i,j));
          my_max = (my_max > current? my_max : current);
        }
    }
  return my_max;
}

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

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

{
  libmesh_assert (this->_m);
  libmesh_assert (this->_n);
  Real my_min = libmesh_real((*this)(0,0));

  for (unsigned int i=0; i!=this->_m; i++)
    {
      for (unsigned int j=0; j!=this->_n; j++)
        {
          Real current = libmesh_real((*this)(i,j));
          my_min = (my_min < current? my_min : current);
        }
    }
  return my_min;
}

template<typename T>

void libMesh::DenseMatrixBase< T >::multiply	(	DenseMatrixBase< T > &	M1,
		const DenseMatrixBase< T > &	M2,
		const DenseMatrixBase< T > &	M3
	)		`[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().

{
  // Assertions to make sure we have been
  // passed matrices of the correct dimension.
  libmesh_assert_equal_to (M1.m(), M2.m());
  libmesh_assert_equal_to (M1.n(), M3.n());
  libmesh_assert_equal_to (M2.n(), M3.m());

  const unsigned int m_s = M2.m();
  const unsigned int p_s = M2.n();
  const unsigned int n_s = M1.n();

  // Do it this way because there is a
  // decent chance (at least for constraint matrices)
  // that M3(k,j) = 0. when right-multiplying.
  for (unsigned int k=0; k<p_s; k++)
    for (unsigned int j=0; j<n_s; j++)
      if (M3.el(k,j) != 0.)
        for (unsigned int i=0; i<m_s; i++)
          M1.el(i,j) += M2.el(i,k) * M3.el(k,j);
}

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.

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

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

{
  for (unsigned int i=0; i<_val.size(); i++)
    if (_val[i] != mat._val[i])
      return true;

  return false;
}

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

{
  libmesh_assert_less (i*j, _val.size());
  libmesh_assert_less (i, this->_m);
  libmesh_assert_less (j, this->_n);


  //  return _val[(i) + (this->_m)*(j)]; // col-major
  return _val[(i)*(this->_n) + (j)]; // row-major
}

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

{
  libmesh_assert_less (i*j, _val.size());
  libmesh_assert_less (i, this->_m);
  libmesh_assert_less (j, this->_n);

  //return _val[(i) + (this->_m)*(j)]; // col-major
  return _val[(i)*(this->_n) + (j)]; // row-major
}

template<typename T>

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

Multiplies every element in the matrix by factor.

Definition at line 789 of file dense_matrix.h.

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

{
  this->scale(factor);
  return *this;
}

template<typename T>

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

Adds mat to this matrix.

Definition at line 843 of file dense_matrix.h.

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

{
  for (unsigned int i=0; i<_val.size(); i++)
    _val[i] += mat._val[i];

  return *this;
}

template<typename T>

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

Subtracts mat from this matrix.

Definition at line 855 of file dense_matrix.h.

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

{
  for (unsigned int i=0; i<_val.size(); i++)
    _val[i] -= mat._val[i];

  return *this;
}

template<typename T>

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

Assignment operator.

Definition at line 722 of file dense_matrix.h.

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

{
  this->_m = other_matrix._m;
  this->_n = other_matrix._n;

  _val                = other_matrix._val;
  _decomposition_type = other_matrix._decomposition_type;

  return *this;
}

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

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

{
  unsigned int mat_m = mat.m(), mat_n = mat.n();
  this->resize(mat_m, mat_n);
  for (unsigned int i=0; i<mat_m; i++)
    for (unsigned int j=0; j<mat_n; j++)
      (*this)(i,j) = mat(i,j);

  return *this;
}

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

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

{
  for (unsigned int i=0; i<_val.size(); i++)
    if (_val[i] != mat._val[i])
      return false;

  return true;
}

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.

{
  for (unsigned int i=0; i<this->m(); i++)
    {
      for (unsigned int j=0; j<this->n(); j++)
        os << std::setw(8)
           << this->el(i,j) << " ";

      os << std::endl;
    }

  return;
}

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.

{
#ifndef LIBMESH_BROKEN_IOSTREAM

  // save the initial format flags
  std::ios_base::fmtflags os_flags = os.flags();

  // Print the matrix entries.
  for (unsigned int i=0; i<this->m(); i++)
    {
      for (unsigned int j=0; j<this->n(); j++)
        os << std::setw(15)
           << std::scientific
           << std::setprecision(8)
           << this->el(i,j) << " ";

      os << std::endl;
    }

  // reset the original format flags
  os.flags(os_flags);

#else

  // Print the matrix entries.
  for (unsigned int i=0; i<this->m(); i++)
    {
      for (unsigned int j=0; j<this->n(); j++)
        os << std::setprecision(8)
           << this->el(i,j)
           << " ";

      os << std::endl;
    }


#endif
}

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

References libMesh::zero.

Referenced by libMesh::DenseMatrix< T >::_lu_decompose(), 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::FESubdivision::init_subdivision_matrix(), libMesh::DGFEMContext::neighbor_side_fe_reinit(), libMesh::ProjectSolution::operator()(), libMesh::ProjectFEMSolution::operator()(), libMesh::BoundaryProjectSolution::operator()(), libMesh::FEMContext::pre_fe_reinit(), and libMesh::HPCoarsenTest::select_refinement().

{
  _val.resize(new_m*new_n);

  this->_m = new_m;
  this->_n = new_n;

  // zero and set decomposition_type to NONE
  this->zero();
}

template<typename T>

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

Right multiplies by the matrix M2.

Implements libMesh::DenseMatrixBase< T >.

Definition at line 210 of file dense_matrix.C.

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

Referenced by libMesh::DofMap::build_constraint_matrix(), libMesh::DofMap::build_constraint_matrix_and_vector(), libMesh::DofMap::constrain_element_matrix(), libMesh::DofMap::constrain_element_matrix_and_vector(), libMesh::DofMap::heterogenously_constrain_element_matrix_and_vector(), and libMesh::FESubdivision::init_shape_functions().

{
  if (this->use_blas_lapack)
    this->_multiply_blas(M3, RIGHT_MULTIPLY);
  else
    {
      // (*this) <- M3 * (*this)
      // Where:
      // (*this) = (m x n),
      // M2      = (m x p),
      // M3      = (p x n)

      // M2 is a copy of *this before it gets resize()d
      DenseMatrix<T> M2(*this);

      // Resize *this so that the result can fit
      this->resize (M2.m(), M3.n());

      this->multiply(*this, M2, M3);
    }
}

template<typename T >

template<typename T2 >

void libMesh::DenseMatrix< T >::right_multiply ( const DenseMatrixBase< T2 > & M2 )

Right multiplies by the matrix M2 of different type

Definition at line 236 of file dense_matrix.C.

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

{
  // (*this) <- M3 * (*this)
  // Where:
  // (*this) = (m x n),
  // M2      = (m x p),
  // M3      = (p x n)

  // M2 is a copy of *this before it gets resize()d
  DenseMatrix<T> M2(*this);

  // Resize *this so that the result can fit
  this->resize (M2.m(), M3.n());

  this->multiply(*this, M2, M3);
}

template<typename T>

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

Right multiplies by the transpose of the matrix A

Definition at line 257 of file dense_matrix.C.

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

{
  if (this->use_blas_lapack)
    this->_multiply_blas(B, RIGHT_MULTIPLY_TRANSPOSE);
  else
    {
      //Check to see if we are doing B*(B^T)
      if (this == &B)
        {
          //libmesh_here();
          DenseMatrix<T> A(*this);

          // Simple but inefficient way
          // return this->right_multiply_transpose(A);

          // More efficient, more code
          // If B is mxn, the result will be a square matrix of Size m x m.
          const unsigned int n_rows = B.m();
          const unsigned int n_cols = B.n();

          // resize() *this and also zero out all entries.
          this->resize(n_rows,n_rows);

          // Compute the lower-triangular part
          for (unsigned int i=0; i<n_rows; ++i)
            for (unsigned int j=0; j<=i; ++j)
              for (unsigned int k=0; k<n_cols; ++k) // inner products are over n_cols
                (*this)(i,j) += A(i,k)*A(j,k);

          // Copy lower-triangular part into upper-triangular part
          for (unsigned int i=0; i<n_rows; ++i)
            for (unsigned int j=i+1; j<n_rows; ++j)
              (*this)(i,j) = (*this)(j,i);
        }

      else
        {
          DenseMatrix<T> A(*this);

          this->resize (A.m(), B.m());

          libmesh_assert_equal_to (A.n(), B.n());
          libmesh_assert_equal_to (this->m(), A.m());
          libmesh_assert_equal_to (this->n(), B.m());

          const unsigned int m_s = A.m();
          const unsigned int p_s = A.n();
          const unsigned int n_s = this->n();

          // Do it this way because there is a
          // decent chance (at least for constraint matrices)
          // that B.transpose(k,j) = 0.
          for (unsigned int j=0; j<n_s; j++)
            for (unsigned int k=0; k<p_s; k++)
              if (B.transpose(k,j) != 0.)
                for (unsigned int i=0; i<m_s; i++)
                  (*this)(i,j) += A(i,k)*B.transpose(k,j);
        }
    }
}

template<typename T >

template<typename T2 >

void libMesh::DenseMatrix< T >::right_multiply_transpose ( const DenseMatrix< T2 > & A )

Right multiplies by the transpose of the matrix A of different type.

Definition at line 322 of file dense_matrix.C.

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

{
  //Check to see if we are doing B*(B^T)
  if (this == &B)
    {
      //libmesh_here();
      DenseMatrix<T> A(*this);

      // Simple but inefficient way
      // return this->right_multiply_transpose(A);

      // More efficient, more code
      // If B is mxn, the result will be a square matrix of Size m x m.
      const unsigned int n_rows = B.m();
      const unsigned int n_cols = B.n();

      // resize() *this and also zero out all entries.
      this->resize(n_rows,n_rows);

      // Compute the lower-triangular part
      for (unsigned int i=0; i<n_rows; ++i)
        for (unsigned int j=0; j<=i; ++j)
          for (unsigned int k=0; k<n_cols; ++k) // inner products are over n_cols
            (*this)(i,j) += A(i,k)*A(j,k);

      // Copy lower-triangular part into upper-triangular part
      for (unsigned int i=0; i<n_rows; ++i)
        for (unsigned int j=i+1; j<n_rows; ++j)
          (*this)(i,j) = (*this)(j,i);
    }

  else
    {
      DenseMatrix<T> A(*this);

      this->resize (A.m(), B.m());

      libmesh_assert_equal_to (A.n(), B.n());
      libmesh_assert_equal_to (this->m(), A.m());
      libmesh_assert_equal_to (this->n(), B.m());

      const unsigned int m_s = A.m();
      const unsigned int p_s = A.n();
      const unsigned int n_s = this->n();

      // Do it this way because there is a
      // decent chance (at least for constraint matrices)
      // that B.transpose(k,j) = 0.
      for (unsigned int j=0; j<n_s; j++)
        for (unsigned int k=0; k<p_s; k++)
          if (B.transpose(k,j) != 0.)
            for (unsigned int i=0; i<m_s; i++)
              (*this)(i,j) += A(i,k)*B.transpose(k,j);
    }
}

template<typename T>

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

Multiplies every element in the matrix by factor.

Definition at line 770 of file dense_matrix.h.

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

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

{
  for (unsigned int i=0; i<this->m(); i++)
    (*this)(i, col) *= factor;
}

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

{
  // We use the LAPACK svd implementation
  _svd_lapack(sigma);
}

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

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

template<typename T>

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

STL-like swap method

Definition at line 665 of file dense_matrix.h.

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

Referenced by libMesh::EulerSolver::_general_residual(), and libMesh::Euler2Solver::_general_residual().

{
  std::swap(this->_m, other_matrix._m);
  std::swap(this->_n, other_matrix._n);
  _val.swap(other_matrix._val);
  DecompositionType _temp = _decomposition_type;
  _decomposition_type = other_matrix._decomposition_type;
  other_matrix._decomposition_type = _temp;
}

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

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

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

template<typename T>

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

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

Definition at line 382 of file dense_matrix.C.

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

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

{
  // Make sure the input sizes are compatible
  libmesh_assert_equal_to (this->n(), arg.size());

  // Resize and clear dest.
  // Note: DenseVector::resize() also zeros the vector.
  dest.resize(this->m());

  // Short-circuit if the matrix is empty
  if(this->m() == 0 || this->n() == 0)
    return;

  if (this->use_blas_lapack)
    this->_matvec_blas(1., 0., dest, arg);
  else
    {
      const unsigned int n_rows = this->m();
      const unsigned int n_cols = this->n();

      for(unsigned int i=0; i<n_rows; i++)
        for(unsigned int j=0; j<n_cols; j++)
          dest(i) += (*this)(i,j)*arg(j);
    }
}

template<typename T>

template<typename T2 >

void libMesh::DenseMatrix< T >::vector_mult	(	DenseVector< typename CompareTypes< T, T2 >::supertype > &	dest,
		const DenseVector< T2 > &	arg
	)		const

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

Definition at line 413 of file dense_matrix.C.

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

{
  // Make sure the input sizes are compatible
  libmesh_assert_equal_to (this->n(), arg.size());

  // Resize and clear dest.
  // Note: DenseVector::resize() also zeros the vector.
  dest.resize(this->m());

  // Short-circuit if the matrix is empty
  if(this->m() == 0 || this->n() == 0)
    return;

  const unsigned int n_rows = this->m();
  const unsigned int n_cols = this->n();

  for(unsigned int i=0; i<n_rows; i++)
    for(unsigned int j=0; j<n_cols; j++)
      dest(i) += (*this)(i,j)*arg(j);
}

template<typename T>

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

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

Definition at line 508 of file dense_matrix.C.

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

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

{
  // Short-circuit if the matrix is empty
  if(this->m() == 0)
    {
      dest.resize(0);
      return;
    }

  if (this->use_blas_lapack)
    this->_matvec_blas(factor, 1., dest, arg);
  else
    {
      DenseVector<T> temp(arg.size());
      this->vector_mult(temp, arg);
      dest.add(factor, temp);
    }
}

template<typename T>

template<typename T2 , typename T3 >

void libMesh::DenseMatrix< T >::vector_mult_add	(	DenseVector< typename CompareTypes< T, typename CompareTypes< T2, T3 >::supertype >::supertype > &	dest,
		const T2	factor,
		const DenseVector< T3 > &	arg
	)		const

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

Definition at line 533 of file dense_matrix.C.

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

{
  // Short-circuit if the matrix is empty
  if (this->m() == 0)
    {
      dest.resize(0);
      return;
    }

  DenseVector<typename CompareTypes<T,T3>::supertype>
    temp(arg.size());
  this->vector_mult(temp, arg);
  dest.add(factor, temp);
}

template<typename T>

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

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

Definition at line 438 of file dense_matrix.C.

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

Referenced by libMesh::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().

{
  // Make sure the input sizes are compatible
  libmesh_assert_equal_to (this->m(), arg.size());

  // Resize and clear dest.
  // Note: DenseVector::resize() also zeros the vector.
  dest.resize(this->n());

  // Short-circuit if the matrix is empty
  if(this->m() == 0)
    return;

  if (this->use_blas_lapack)
    {
      this->_matvec_blas(1., 0., dest, arg, /*trans=*/true);
    }
  else
    {
      const unsigned int n_rows = this->m();
      const unsigned int n_cols = this->n();

      // WORKS
      // for(unsigned int j=0; j<n_cols; j++)
      //   for(unsigned int i=0; i<n_rows; i++)
      //     dest(j) += (*this)(i,j)*arg(i);

      // ALSO WORKS, (i,j) just swapped
      for(unsigned int i=0; i<n_cols; i++)
        for(unsigned int j=0; j<n_rows; j++)
          dest(i) += (*this)(j,i)*arg(j);
    }
}

template<typename T>

template<typename T2 >

void libMesh::DenseMatrix< T >::vector_mult_transpose	(	DenseVector< typename CompareTypes< T, T2 >::supertype > &	dest,
		const DenseVector< T2 > &	arg
	)		const

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

Definition at line 477 of file dense_matrix.C.

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

{
  // Make sure the input sizes are compatible
  libmesh_assert_equal_to (this->m(), arg.size());

  // Resize and clear dest.
  // Note: DenseVector::resize() also zeros the vector.
  dest.resize(this->n());

  // Short-circuit if the matrix is empty
  if(this->m() == 0)
    return;

  const unsigned int n_rows = this->m();
  const unsigned int n_cols = this->n();

  // WORKS
  // for(unsigned int j=0; j<n_cols; j++)
  //   for(unsigned int i=0; i<n_rows; i++)
  //     dest(j) += (*this)(i,j)*arg(i);

  // ALSO WORKS, (i,j) just swapped
  for(unsigned int i=0; i<n_cols; i++)
    for(unsigned int j=0; j<n_rows; j++)
      dest(i) += (*this)(j,i)*arg(j);
}

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

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

{
  _decomposition_type = NONE;

  std::fill (_val.begin(), _val.end(), static_cast<T>(0));
}

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.

  {
    m.print(os);
    return os;
  }

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 487 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 [protected, inherited]

The row dimension.

Definition at line 165 of file dense_matrix_base.h.

Referenced by libMesh::DenseMatrixBase< Number >::m(), libMesh::DenseMatrix< T >::operator=(), 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::DenseMatrixBase< Number >::n(), libMesh::DenseMatrix< T >::operator=(), 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 559 of file dense_matrix.h.

template<typename T>

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

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

Definition at line 439 of file dense_matrix.h.

Referenced by libMesh::DenseMatrix< T >::_multiply_blas(), libMesh::DenseMatrix< Number >::get_values(), libMesh::DenseMatrix< T >::operator!=(), libMesh::DenseMatrix< T >::operator+=(), libMesh::DenseMatrix< T >::operator-=(), libMesh::DenseMatrix< T >::operator=(), libMesh::DenseMatrix< T >::operator==(), and libMesh::DenseMatrix< T >::swap().

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

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

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 >