DenseMatrix< T > Class Template Reference

#include <dense_matrix.h>

Inheritance diagram for DenseMatrix< T >:

[legend]

List of all members.

Public Member Functions
	DenseMatrix (const unsigned int m=0, const unsigned int n=0)
virtual	~DenseMatrix ()
virtual void	zero ()
T	operator() (const unsigned int i, const unsigned int j) const
T &	operator() (const unsigned int i, const unsigned int j)
virtual T	el (const unsigned int i, const unsigned int j) const
virtual T &	el (const unsigned int i, const unsigned int j)
virtual void	left_multiply (const DenseMatrixBase< T > &M2)
virtual void	right_multiply (const DenseMatrixBase< T > &M3)
void	vector_mult (DenseVector< T > &dest, const DenseVector< T > &arg) const
void	vector_mult_add (DenseVector< T > &dest, const T factor, const DenseVector< T > &arg) const
void	get_principal_submatrix (unsigned int sub_m, unsigned int sub_n, DenseMatrix< T > &dest) const
void	get_principal_submatrix (unsigned int sub_m, DenseMatrix< T > &dest) const
DenseMatrix< T > &	operator= (const DenseMatrix< T > &other_matrix)
void	swap (DenseMatrix< T > &other_matrix)
void	resize (const unsigned int m, const unsigned int n)
void	scale (const T factor)
DenseMatrix< T > &	operator*= (const T factor)
void	add (const T factor, const DenseMatrix< T > &mat)
bool	operator== (const DenseMatrix< T > &mat) const
bool	operator!= (const DenseMatrix< T > &mat) const
DenseMatrix< T > &	operator+= (const DenseMatrix< T > &mat)
DenseMatrix< T > &	operator-= (const DenseMatrix< T > &mat)
Real	min () const
Real	max () const
Real	l1_norm () const
Real	linfty_norm () const
void	left_multiply_transpose (const DenseMatrix< T > &A)
void	right_multiply_transpose (const DenseMatrix< T > &A)
T	transpose (const unsigned int i, const unsigned int j) const
void	get_transpose (DenseMatrix< T > &dest) const
std::vector< T > &	get_values ()
const std::vector< T > &	get_values () const
void	condense (const unsigned int i, const unsigned int j, const T val, DenseVector< T > &rhs)
void	lu_solve (DenseVector< T > &b, DenseVector< T > &x, const bool partial_pivot=false)
template<typename T2 >
void	cholesky_solve (DenseVector< T2 > &b, DenseVector< T2 > &x)
T	det ()
unsigned int	m () const
unsigned int	n () const
void	print (std::ostream &os) 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 (const bool partial_pivot=false)
void	_lu_back_substitute (DenseVector< T > &b, DenseVector< T > &x, const bool partial_pivot=false) const
void	_cholesky_decompose ()
template<typename T2 >
void	_cholesky_back_substitute (DenseVector< T2 > &b, DenseVector< T2 > &x) const
void	_multiply_blas (const DenseMatrixBase< T > &other, _BLAS_Multiply_Flag flag)
void	_lu_decompose_lapack ()
void	_lu_back_substitute_lapack (DenseVector< T > &b, DenseVector< T > &x)
void	_matvec_blas (T alpha, T beta, DenseVector< T > &dest, const DenseVector< T > &arg) 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 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 49 of file dense_matrix.h.

---

Member Enumeration Documentation

template<typename T >

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

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

template<typename T >

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

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

---

Constructor & Destructor Documentation

template<typename T >

DenseMatrix< T >::DenseMatrix	(	const unsigned int	m = 0,
		const unsigned int	n = 0
	)			[inline]

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

Definition at line 482 of file dense_matrix.h.

References DenseMatrix< T >::resize().

  : DenseMatrixBase<T>(m,n),
#if defined(LIBMESH_HAVE_PETSC) && defined(LIBMESH_USE_REAL_NUMBERS)
    use_blas_lapack(true),
#else
    use_blas_lapack(false),
#endif
    _decomposition_type(NONE)
{
  this->resize(m,n);
}

template<typename T >

virtual DenseMatrix< T >::~DenseMatrix

(

)

[inline, virtual]

Copy-constructor. Destructor. Empty.

Definition at line 67 of file dense_matrix.h.

{}

---

Member Function Documentation

template<typename T >

template<typename T2 >

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

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

Definition at line 617 of file dense_matrix.C.

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

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

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

  // Just to be really sure...
  libmesh_assert(m==n);

  // A convenient reference to *this
  const DenseMatrix<T>& A = *this;
   
  // Now compute the solution to Ax =b using the factorization.
  x.resize(m);
  
  // Solve for Ly=b
  for (unsigned int i=0; i<n; ++i)
    {
      for (unsigned int k=0; k<i; ++k)
        b(i) -= A(i,k)*x(k);
      
      x(i) = b(i) / A(i,i);
    }
  
  // Solve for L^T x = y
  for (unsigned int i=0; i<n; ++i)
    {
      const unsigned int ib = (n-1)-i;

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

template<typename T >

void DenseMatrix< T >::_cholesky_decompose

(

)

[inline, private]

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

Definition at line 566 of file dense_matrix.C.

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

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

{
  // If we called this function, there better not be any
  // previous decomposition of the matrix.
  libmesh_assert(this->_decomposition_type == NONE);
  
  // Shorthand notation for number of rows and columns.
  const unsigned int
    m = this->m(),
    n = this->n();

  // Just to be really sure...
  libmesh_assert(m==n);

  // A convenient reference to *this
  DenseMatrix<T>& A = *this;
  
  for (unsigned int i=0; i<m; ++i)
    {
      for (unsigned int j=i; j<n; ++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)
                {
                  std::cerr << "Error! Can only use Cholesky decomposition "
                            << "with symmetric positive definite matrices."
                            << std::endl;
                  libmesh_error();
                }
#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 >

void DenseMatrix< T >::_lu_back_substitute	(	DenseVector< T > &	b,
		DenseVector< T > &	x,
		const bool	partial_pivot = false
	)			const [inline, private]

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

Definition at line 368 of file dense_matrix.C.

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

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

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

  libmesh_assert (this->m() == n);
  libmesh_assert (this->m() == b.size());
  
  x.resize (n);

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

  // Transform the RHS vector
  for (unsigned int i=0; i<(n-1); i++)
    {
      // Get the diagonal entry and take its inverse
      const T diag = A(i,i);
      
      libmesh_assert (diag != libMesh::zero);
  
      const T diag_inv = 1./diag;

      // Get the entry b(i) and store it
      const T bi = b(i);
      
      for (unsigned int j=(i+1); j<n; j++)
        b(j) -= bi*A(j,i)*diag_inv;
    }


  // Perform back-substitution
  {
    x(n-1) = b(n-1)/A(n-1,n-1);

    for (unsigned int i=0; i<=(n-1); i++)
      {
        const unsigned int ib = (n-1)-i;

        // Get the diagonal and take its inverse
        const T diag = A(ib,ib);

        libmesh_assert (diag != libMesh::zero);

        const T diag_inv = 1./diag;
        
        for (unsigned int j=(ib+1); j<n; j++)
          {
            b(ib) -= A(ib,j)*x(j);
            x(ib)  = b(ib)*diag_inv;
          }
      }
  }
}

template<typename T >

void DenseMatrix< T >::_lu_back_substitute_lapack	(	DenseVector< T > &	b,
		DenseVector< T > &	x
	)			[inline, private]

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

Definition at line 297 of file dense_matrix_blas_lapack.C.

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

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

{
  // 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)
    {
      std::cout << "INFO="
                << INFO
                << ", Error during Lapack LU solve!" << std::endl;
      libmesh_error();
    }

  // 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 DenseMatrix< T >::_lu_decompose

(

const bool

partial_pivot = false

)

[inline, private]

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

Definition at line 433 of file dense_matrix.C.

References DenseMatrix< T >::_decomposition_type, DenseMatrix< T >::LU, DenseMatrixBase< T >::m(), DenseMatrix< T >::NONE, and libMesh::zero.

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

{
  // If this function was called, there better not be any
  // previous decomposition of the matrix.
  libmesh_assert(this->_decomposition_type == NONE);
  
  // Get the matrix size and make sure it is square
  const unsigned int
    m = this->m();

  // A convenient reference to *this
  DenseMatrix<T>& A = *this;
  
  // Straight, vanilla LU factorization without pivoting
  if (!partial_pivot)
    {
      
      // For each row in the matrix
      for (unsigned int i=0; i<m; i++)
        {
          // Get the diagonal entry and take its inverse
          const T diag = A(i,i);

          libmesh_assert (diag != libMesh::zero);

          const T diag_inv = 1./diag;
          
          // For each row in the submatrix
          for (unsigned int j=i+1; j<m; j++)
            {
              // Get the scale factor for this row
              const T fact = A(j,i)*diag_inv;
              
              // For each column in the subrow scale it
              // by the factor
              for (unsigned int k=i+1; k<m; k++)
                A(j,k) -= fact*A(i,k);        
            }
        }
    }
  
  // Do partial pivoting.
  else
    {
      libmesh_error();
    }
  
  // Set the flag for LU decomposition
  this->_decomposition_type = LU;
}

template<typename T >

void DenseMatrix< T >::_lu_decompose_lapack

(

)

[inline, private]

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

Definition at line 219 of file dense_matrix_blas_lapack.C.

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

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

{
  // If this function was called, there better not be any
  // previous decomposition of the matrix.
  libmesh_assert(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)
    {
      std::cout << "INFO="
                << INFO
                << ", Error during Lapack LU factorization!" << std::endl;
      libmesh_error();
    }
  
  // Set the flag for LU decomposition
  this->_decomposition_type = LU_BLAS_LAPACK;
}

template<typename T >

void DenseMatrix< T >::_matvec_blas	(	T	alpha,
		T	beta,
		DenseVector< T > &	dest,
		const DenseVector< T > &	arg
	)			const [inline, private]

Uses the BLAS GEMV function (through PETSc) to compute

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

where alpha and beta are scalars, A is this matrix, and arg and dest are input vectors of appropriate size.

[ Implementation in dense_matrix_blas_lapack.C ]

Definition at line 392 of file dense_matrix_blas_lapack.C.

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

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

{
  // Ensure that dest and arg are proper size
  // dest  ~ A     * arg
  // (mx1)   (mxn) * (nx1)
  if ((dest.size() != this->m()) || (arg.size() != this->n()))
    {
      std::cout << "Improper input argument sizes!" << std::endl;
      libmesh_error();
    }
  
  // 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
  char TRANS[] = "t";
  
  //   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 >

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

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

Definition at line 39 of file dense_matrix_blas_lapack.C.

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

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

{
  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:
      {
        std::cout << "Unknown flag selected or matrices are ";
        std::cout << "incompatible for multiplication." << std::endl;
        libmesh_error();
      }
    }

  // For this to work, the passed arg. must actually be a DenseMatrix<T>
  const DenseMatrix<T>* const_that = dynamic_cast< const DenseMatrix<T>* >(&other);
  if (!const_that)
    {
      std::cerr << "Unable to cast input matrix to usable type." << std::endl;
      libmesh_error();
    }

  // 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:
      {
        std::cout << "Unknown flag selected." << std::endl;
        libmesh_error();
      }
    }
  
  // Swap my data vector with the result
  this->_val.swap(result);
}

template<typename T >

template<typename T2 , typename T3 >

boostcopy::enable_if_c< ScalarTraits< T2 >::value, void >::type 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 183 of file dense_matrix_base.h.

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

{
  libmesh_assert (this->m() == mat.m());
  libmesh_assert (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 >

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

Adds factor times mat to this matrix.

Definition at line 621 of file dense_matrix.h.

References DenseMatrix< T >::_val.

Referenced by FEMSystem::assembly().

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

template<typename T >

template<typename T2 >

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

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

Definition at line 529 of file dense_matrix.C.

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

Referenced by FEBase::coarsened_dof_values(), FEBase::compute_periodic_constraints(), FEBase::compute_proj_constraints(), and System::ProjectVector::operator()().

{
  // 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:
      {
        std::cerr << "Error! This matrix already has a "
                  << "different decomposition..."
                  << std::endl;
        libmesh_error();
      }
    }

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

template<typename T >

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

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

Definition at line 57 of file dense_matrix_base.C.

References DenseMatrixBase< T >::_m, DenseMatrixBase< T >::el(), DenseVectorBase< T >::el(), DenseMatrixBase< T >::m(), DenseMatrixBase< T >::n(), and DenseVectorBase< T >::size().

{
  libmesh_assert (this->_m == rhs.size());
  libmesh_assert (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 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 265 of file dense_matrix.h.

Referenced by DenseMatrix< Real >::condense().

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

template<typename T >

T DenseMatrix< T >::det

(

)

[inline]

Returns:

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

Definition at line 492 of file dense_matrix.C.

References DenseMatrix< T >::_decomposition_type, DenseMatrix< T >::_lu_decompose(), DenseMatrix< T >::LU, DenseMatrixBase< T >::m(), and DenseMatrix< T >::NONE.

{
  // First LU decompose the matrix (without partial pivoting).
  // Note that the lu_decompose routine will check to see if the
  // matrix is square so we don't worry about it.
  if (this->_decomposition_type == NONE)
    this->_lu_decompose(false);
  else if (this->_decomposition_type != LU)
    {
      std::cerr << "Error! Can't compute the determinant under "
                << "the current decomposition."
                << std::endl;
      libmesh_error();
    }
  
  // 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.
  for (unsigned int i=0; i<this->m(); i++)
    determinant *= (*this)(i,i);

  // Return the determinant
  return determinant;
}

template<typename T >

virtual T& 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 DenseMatrixBase< T >.

Definition at line 96 of file dense_matrix.h.

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

template<typename T >

virtual T 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 DenseMatrixBase< T >.

Definition at line 90 of file dense_matrix.h.

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

template<typename T >

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

Put the sub_m x sub_m principal submatrix into dest.

Definition at line 278 of file dense_matrix.C.

References DenseMatrix< T >::get_principal_submatrix().

{
  get_principal_submatrix(sub_m, sub_m, dest);
}

template<typename T >

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

Put the sub_m x sub_n principal submatrix into dest.

Definition at line 263 of file dense_matrix.C.

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

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

{
  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 DenseMatrix< T >::get_transpose

(

DenseMatrix< T > &

dest

)

const [inline]

Put the tranposed matrix into dest.

Definition at line 286 of file dense_matrix.C.

References DenseMatrixBase< T >::m(), DenseMatrixBase< T >::n(), and 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 >

const std::vector<T>& DenseMatrix< T >::get_values

(

)

const [inline]

Return a constant reference to the matrix values.

Definition at line 257 of file dense_matrix.h.

{ return _val; }

template<typename T >

std::vector<T>& 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 252 of file dense_matrix.h.

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

{ return _val; }

template<typename T >

Real DenseMatrix< T >::l1_norm

(

)

const [inline]

Return the l1-norm of the matrix, that is , (max. sum of columns). This is the natural matrix norm that is compatible to the l1-norm for vectors, i.e. .

Definition at line 723 of file dense_matrix.h.

References DenseMatrixBase< T >::_m, and DenseMatrixBase< T >::_n.

Referenced by 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 DenseMatrix< T >::left_multiply

(

const DenseMatrixBase< T > &

M2

)

[inline, virtual]

Left multipliess by the matrix M2.

Implements DenseMatrixBase< T >.

Definition at line 35 of file dense_matrix.C.

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

{
  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 >

void DenseMatrix< T >::left_multiply_transpose

(

const DenseMatrix< T > &

A

)

[inline]

Left multiplies by the transpose of the matrix A.

Definition at line 62 of file dense_matrix.C.

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

Referenced by DofMap::constrain_element_matrix(), and DofMap::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 m = A.m();
          const unsigned int n = A.n();

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

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

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

      else
        {
          DenseMatrix<T> B(*this);
  
          this->resize (A.n(), B.n());
      
          libmesh_assert (A.m() == B.m());
          libmesh_assert (this->m() == A.n());
          libmesh_assert (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 DenseMatrix< T >::linfty_norm

(

)

const [inline]

Return the linfty-norm of the matrix, that is , (max. sum of rows). This is the natural matrix norm that is compatible to the linfty-norm of vectors, i.e. .

Definition at line 750 of file dense_matrix.h.

References DenseMatrixBase< T >::_m, and DenseMatrixBase< T >::_n.

{
  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 DenseMatrix< T >::lu_solve	(	DenseVector< T > &	b,
		DenseVector< T > &	x,
		const bool	partial_pivot = false
	)			[inline]

Solve the system Ax=b given the input vector b.

Definition at line 299 of file dense_matrix.C.

References DenseMatrix< T >::_decomposition_type, DenseMatrix< T >::_lu_back_substitute(), DenseMatrix< T >::_lu_back_substitute_lapack(), DenseMatrix< T >::_lu_decompose(), DenseMatrix< T >::_lu_decompose_lapack(), DenseMatrix< T >::LU, DenseMatrix< T >::LU_BLAS_LAPACK, DenseMatrixBase< T >::m(), DenseMatrixBase< T >::n(), DenseMatrix< T >::NONE, and DenseMatrix< T >::use_blas_lapack.

Referenced by 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...
  if (this->m() != this->n())
    {
      std::cout << "Error! LU solve only works for square matrices!"
                << std::endl;
      libmesh_error();
    }


  switch(this->_decomposition_type)
    {
    case NONE:
      {
        if (this->use_blas_lapack)
          this->_lu_decompose_lapack();   
        else
          this->_lu_decompose (partial_pivot);
        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:
      {
        std::cerr << "Error! This matrix already has a "
                  << "different decomposition..."
                  << std::endl;
        libmesh_error();
      }
    }

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

template<typename T >

unsigned int DenseMatrixBase< T >::m

(

)

const [inline, inherited]

Returns:

the row-dimension of the matrix.

Definition at line 98 of file dense_matrix_base.h.

References DenseMatrixBase< T >::_m.

Referenced by DenseMatrix< T >::_cholesky_back_substitute(), DenseMatrix< T >::_cholesky_decompose(), DenseMatrix< T >::_lu_back_substitute(), DenseMatrix< T >::_lu_back_substitute_lapack(), DenseMatrix< T >::_lu_decompose(), DenseMatrix< T >::_lu_decompose_lapack(), DenseMatrix< T >::_matvec_blas(), DenseMatrix< T >::_multiply_blas(), DenseMatrixBase< T >::add(), EpetraMatrix< T >::add_matrix(), PetscMatrix< T >::add_matrix(), LaspackMatrix< T >::add_matrix(), DofMap::build_constraint_matrix(), DenseMatrixBase< T >::condense(), DofMap::constrain_element_dyad_matrix(), DofMap::constrain_element_matrix(), DofMap::constrain_element_matrix_and_vector(), DofMap::constrain_element_vector(), DenseMatrix< T >::det(), DofMap::extract_local_vector(), DenseMatrix< T >::get_principal_submatrix(), DenseMatrix< T >::get_transpose(), DenseMatrix< T >::left_multiply(), DenseMatrix< T >::left_multiply_transpose(), DenseMatrix< T >::lu_solve(), DofMap::max_constraint_error(), DenseMatrixBase< T >::multiply(), PatchRecoveryErrorEstimator::EstimateError::operator()(), DenseMatrixBase< T >::print(), DenseMatrixBase< T >::print_scientific(), DenseMatrix< T >::right_multiply(), DenseMatrix< T >::right_multiply_transpose(), and DenseMatrix< T >::vector_mult().

{ return _m; }

template<typename T >

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

References DenseMatrixBase< T >::_m, DenseMatrixBase< T >::_n, 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 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 681 of file dense_matrix.h.

References DenseMatrixBase< T >::_m, DenseMatrixBase< T >::_n, 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 DenseMatrixBase< T >::multiply	(	DenseMatrixBase< T > &	M1,
		const DenseMatrixBase< T > &	M2,
		const DenseMatrixBase< T > &	M3
	)			[inline, protected, inherited]

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

Definition at line 30 of file dense_matrix_base.C.

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

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

{
  // Assertions to make sure we have been
  // passed matrices of the correct dimension.
  libmesh_assert (M1.m() == M2.m());
  libmesh_assert (M1.n() == M3.n());
  libmesh_assert (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 DenseMatrixBase< T >::n

(

)

const [inline, inherited]

Returns:

the column-dimension of the matrix.

Definition at line 103 of file dense_matrix_base.h.

References DenseMatrixBase< T >::_n.

Referenced by DenseMatrix< T >::_cholesky_back_substitute(), DenseMatrix< T >::_cholesky_decompose(), DenseMatrix< T >::_lu_back_substitute(), DenseMatrix< T >::_lu_decompose_lapack(), DenseMatrix< T >::_matvec_blas(), DenseMatrix< T >::_multiply_blas(), DenseMatrixBase< T >::add(), EpetraMatrix< T >::add_matrix(), PetscMatrix< T >::add_matrix(), LaspackMatrix< T >::add_matrix(), DofMap::build_constraint_matrix(), DenseMatrixBase< T >::condense(), DofMap::constrain_element_dyad_matrix(), DofMap::constrain_element_matrix(), DofMap::constrain_element_matrix_and_vector(), DofMap::constrain_element_vector(), DofMap::extract_local_vector(), DenseMatrix< T >::get_principal_submatrix(), DenseMatrix< T >::get_transpose(), DenseMatrix< T >::left_multiply(), DenseMatrix< T >::left_multiply_transpose(), DenseMatrix< T >::lu_solve(), DofMap::max_constraint_error(), DenseMatrixBase< T >::multiply(), PatchRecoveryErrorEstimator::EstimateError::operator()(), DenseMatrixBase< T >::print(), DenseMatrixBase< T >::print_scientific(), DenseMatrix< T >::right_multiply(), DenseMatrix< T >::right_multiply_transpose(), and DenseMatrix< T >::vector_mult().

{ return _n; }

template<typename T >

bool DenseMatrix< T >::operator!=

(

const DenseMatrix< T > &

mat

)

const [inline]

Tests if mat is not exactly equal to this matrix.

Definition at line 644 of file dense_matrix.h.

References 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 & 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 584 of file dense_matrix.h.

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

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

template<typename T >

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

Returns:

the (i,j) element of the matrix.

Definition at line 568 of file dense_matrix.h.

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

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

template<typename T >

DenseMatrix< T > & DenseMatrix< T >::operator*=

(

const T

factor

)

[inline]

Multiplies every element in the matrix by factor.

Definition at line 611 of file dense_matrix.h.

References DenseMatrix< T >::scale().

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

template<typename T >

DenseMatrix< T > & DenseMatrix< T >::operator+=

(

const DenseMatrix< T > &

mat

)

[inline]

Adds mat to this matrix.

Definition at line 657 of file dense_matrix.h.

References DenseMatrix< T >::_val.

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

  return *this;
}

template<typename T >

DenseMatrix< T > & DenseMatrix< T >::operator-=

(

const DenseMatrix< T > &

mat

)

[inline]

Subtracts mat from this matrix.

Definition at line 669 of file dense_matrix.h.

References DenseMatrix< T >::_val.

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

  return *this;
}

template<typename T >

DenseMatrix< T > & DenseMatrix< T >::operator=

(

const DenseMatrix< T > &

other_matrix

)

[inline]

Assignment operator.

Definition at line 553 of file dense_matrix.h.

References DenseMatrix< T >::_decomposition_type, DenseMatrixBase< T >::_m, DenseMatrixBase< T >::_n, and 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 >

bool DenseMatrix< T >::operator==

(

const DenseMatrix< T > &

mat

)

const [inline]

Tests if mat is exactly equal to this matrix.

Definition at line 631 of file dense_matrix.h.

References 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 DenseMatrixBase< T >::print

(

std::ostream &

os

)

const [inline, inherited]

Pretty-print the matrix to stdout.

Definition at line 127 of file dense_matrix_base.C.

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

{  
  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 DenseMatrixBase< T >::print_scientific

(

std::ostream &

os

)

const [inline, inherited]

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

Definition at line 85 of file dense_matrix_base.C.

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

{
#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 DenseMatrix< T >::resize	(	const unsigned int	m,
		const unsigned int	n
	)			[inline]

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

Definition at line 526 of file dense_matrix.h.

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

Referenced by DofMap::build_constraint_matrix(), FEBase::coarsened_dof_values(), FEBase::compute_periodic_constraints(), FEBase::compute_proj_constraints(), DenseMatrix< T >::DenseMatrix(), DenseMatrix< T >::get_principal_submatrix(), DenseMatrix< T >::get_transpose(), DenseMatrix< T >::left_multiply(), DenseMatrix< T >::left_multiply_transpose(), System::ProjectVector::operator()(), DenseMatrix< T >::right_multiply(), and DenseMatrix< T >::right_multiply_transpose().

{
  _val.resize(m*n);

  this->_m = m;
  this->_n = n;

  _decomposition_type = NONE;
  this->zero();
}

template<typename T >

void DenseMatrix< T >::right_multiply

(

const DenseMatrixBase< T > &

M3

)

[inline, virtual]

Right multiplies by the matrix M3.

Implements DenseMatrixBase< T >.

Definition at line 130 of file dense_matrix.C.

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

Referenced by DofMap::build_constraint_matrix(), DofMap::constrain_element_matrix(), and DofMap::constrain_element_matrix_and_vector().

{
  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 >

void DenseMatrix< T >::right_multiply_transpose

(

const DenseMatrix< T > &

A

)

[inline]

Right multiplies by the transpose of the matrix A

Definition at line 156 of file dense_matrix.C.

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

{
  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 m = B.m();
          const unsigned int n = B.n();

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

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

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

      else
        {
          DenseMatrix<T> A(*this);
  
          this->resize (A.m(), B.m());
      
          libmesh_assert (A.n() == B.n());
          libmesh_assert (this->m() == A.m());
          libmesh_assert (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 DenseMatrix< T >::scale

(

const T

factor

)

[inline]

Multiplies every element in the matrix by factor.

Definition at line 601 of file dense_matrix.h.

References DenseMatrix< T >::_val.

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

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

template<typename T >

void DenseMatrix< T >::swap

(

DenseMatrix< T > &

other_matrix

)

[inline]

STL-like swap method

Definition at line 512 of file dense_matrix.h.

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

Referenced by DenseMatrix< T >::_multiply_blas(), and FEMSystem::assembly().

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

Referenced by DofMap::constrain_element_dyad_matrix(), DofMap::constrain_element_matrix_and_vector(), DofMap::constrain_element_vector(), DenseMatrix< T >::left_multiply_transpose(), and DenseMatrix< T >::right_multiply_transpose().

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

template<typename T >

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

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

Definition at line 220 of file dense_matrix.C.

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

Referenced by DenseMatrix< T >::vector_mult_add().

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

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

  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 >

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

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

Definition at line 246 of file dense_matrix.C.

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

{
  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 >

void DenseMatrix< T >::zero

(

)

[inline, virtual]

Set every element in the matrix to 0.

Implements DenseMatrixBase< T >.

Definition at line 542 of file dense_matrix.h.

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

Referenced by FEBase::coarsened_dof_values(), System::ProjectVector::operator()(), and DenseMatrix< T >::resize().

{
  _decomposition_type = NONE;

  std::fill (_val.begin(), _val.end(), 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; std::cout << K << std::endl;

Definition at line 115 of file dense_matrix_base.h.

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

---

Member Data Documentation

template<typename T >

DecompositionType DenseMatrix< T >::_decomposition_type [private]

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

Definition at line 372 of file dense_matrix.h.

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

template<typename T >

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

The row dimension.

Definition at line 164 of file dense_matrix_base.h.

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

template<typename T >

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

The column dimension.

Definition at line 169 of file dense_matrix_base.h.

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

template<typename T >

std::vector<int> 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 411 of file dense_matrix.h.

Referenced by DenseMatrix< T >::_lu_back_substitute_lapack(), and DenseMatrix< T >::_lu_decompose_lapack().

template<typename T >

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

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

Definition at line 323 of file dense_matrix.h.

Referenced by DenseMatrix< T >::_lu_back_substitute_lapack(), DenseMatrix< T >::_lu_decompose_lapack(), DenseMatrix< T >::_multiply_blas(), DenseMatrix< T >::add(), DenseMatrix< Real >::get_values(), DenseMatrix< T >::operator!=(), DenseMatrix< T >::operator()(), DenseMatrix< T >::operator+=(), DenseMatrix< T >::operator-=(), DenseMatrix< T >::operator=(), DenseMatrix< T >::operator==(), DenseMatrix< T >::resize(), DenseMatrix< T >::scale(), DenseMatrix< T >::swap(), and DenseMatrix< T >::zero().

template<typename T >

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

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

---

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

• dense_matrix.h
• dense_matrix.C
• dense_matrix_blas_lapack.C