DenseMatrix

Section: C Library Functions (3) Updated: Thu Apr 7 2011 Local index Up

NAME

DenseMatrix -

SYNOPSIS

Inherits DenseMatrixBase< T >.

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

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, 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)

Private Attributes

std::vector< T > _val

DecompositionType _decomposition_type

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]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.

Enumerator:

LEFT_MULTIPLY
RIGHT_MULTIPLY
LEFT_MULTIPLY_TRANSPOSE
RIGHT_MULTIPLY_TRANSPOSE

Definition at line 377 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
NONE

Definition at line 364 of file dense_matrix.h.

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

Constructor & Destructor Documentation

template<typename T > DenseMatrix< T >::DenseMatrix (const unsigned intm = 0, const unsigned intn = 0) [inline]Constructor. Creates a dense matrix of dimension m by n.

Definition at line 430 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(true),
#else
    use_blas(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 [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 561 of file dense_matrix.C.

References DenseVector< T >::resize().

{
  // 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 () [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 510 of file dense_matrix.C.

{
  // 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 boolpartial_pivot = false) 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 321 of file dense_matrix.C.

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

{
  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_decompose (const boolpartial_pivot = false) [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 380 of file dense_matrix.C.

References libMesh::zero.

{
  // 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();

  libmesh_assert (m == this->n());

  // 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 >::_multiply_blas (const DenseMatrixBase< T > &other, _BLAS_Multiply_Flagflag) [private]

Definition at line 606 of file dense_matrix.C.

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

{
  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> void DenseMatrix< T >::add (const Tfactor, const DenseMatrix< T > &mat) [inline]Adds factor times mat to this matrix.

Definition at line 569 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 , typename T3 > boostcopy::enable_if_c< ScalarTraits< T2 >::value, void >::type DenseMatrixBase< T >::add (const T2factor, 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 > template<typename T2 > void DenseMatrix< T >::cholesky_solve (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 473 of file dense_matrix.C.

Referenced by FEBase::coarsened_dof_values(), FEBase::compute_periodic_constraints(), FEBase::compute_proj_constraints(), System::ProjectVector::operator()(), and 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:
      {
        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 DenseMatrix< T >::condense (const unsigned inti, const unsigned intj, const Tval, 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 263 of file dense_matrix.h.

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

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

template<typename T> void DenseMatrixBase< T >::condense (const unsigned inti, const unsigned intj, const Tval, 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 57 of file dense_matrix_base.C.

References DenseVectorBase< T >::el(), 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 > T 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 436 of file dense_matrix.C.

{
  // 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 inti, const unsigned intj) 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> virtual T& DenseMatrix< T >::el (const unsigned inti, const unsigned intj) [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> void DenseMatrix< T >::get_principal_submatrix (unsigned intsub_m, unsigned intsub_n, DenseMatrix< T > &dest) constPut the sub_m x sub_n principal submatrix into dest.

Definition at line 256 of file dense_matrix.C.

References 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 DenseMatrix< T >::get_principal_submatrix (unsigned intsub_m, DenseMatrix< T > &dest) constPut the sub_m x sub_m principal submatrix into dest.

Definition at line 269 of file dense_matrix.C.

{
  get_principal_submatrix(sub_m, sub_m, dest);
}

template<typename T> void DenseMatrix< T >::get_transpose (DenseMatrix< T > &dest) constPut the tranposed matrix into dest.

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

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

{ return _val; }

template<typename T> const std::vector<T>& DenseMatrix< T >::get_values () const [inline]Return a constant reference to the matrix values.

Definition at line 255 of file dense_matrix.h.

{ return _val; }

template<typename T > Real DenseMatrix< T >::l1_norm () const [inline]Return the l1-norm of the matrix, that is $|M|_1=max_{all columns j}um_{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
finition at line 671 of file dense_matrix.h.

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> EXTERN_C_FOR_PETSC_BEGIN EXTERN_C_FOR_PETSC_END void DenseMatrix< T >::left_multiply (const DenseMatrixBase< T > &M2) [virtual]Left multipliess by the matrix M2.

Implements DenseMatrixBase< T >.

Definition at line 44 of file dense_matrix.C.

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

{
  if (this->use_blas)
    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)Left multiplies by the transpose of the matrix A.

Definition at line 71 of file dense_matrix.C.

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

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

{
  if (this->use_blas)
    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 $|M|_infty=max_{all rows i}um_{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
finition at line 698 of file dense_matrix.h.

{
  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 boolpartial_pivot = false)Solve the system Ax=b given the input vector b.

Definition at line 286 of file dense_matrix.C.

Referenced by PatchRecoveryErrorEstimator::EstimateError::operator()().

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

    case LU:
      {
        // 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->_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.

Referenced by DenseMatrix< T >::_multiply_blas(), DenseMatrixBase< T >::add(), EpetraMatrix< T >::add_matrix(), PetscMatrix< T >::add_matrix(), LaspackMatrix< T >::add_matrix(), DofMap::build_constraint_matrix(), 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_transpose(), DenseMatrix< T >::left_multiply(), DenseMatrix< T >::left_multiply_transpose(), DofMap::max_constraint_error(), DenseMatrixBase< T >::multiply(), PatchRecoveryErrorEstimator::EstimateError::operator()(), DenseMatrix< T >::right_multiply(), and DenseMatrix< T >::right_multiply_transpose().

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

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

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

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

{ 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 592 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 inti, const unsigned intj) const [inline]Returns:

: the (i,j) element of the matrix.

Definition at line 516 of file dense_matrix.h.

{
  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 inti, const unsigned intj) [inline]Returns:

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

Definition at line 532 of file dense_matrix.h.

{
  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 Tfactor) [inline]Multiplies every element in the matrix by factor.

Definition at line 559 of file dense_matrix.h.

References MeshTools::Modification::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 605 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 617 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 501 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 579 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 [inherited]Pretty-print the matrix to stdout.

Definition at line 127 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 DenseMatrixBase< T >::print_scientific (std::ostream &os) const [inherited]Prints the matrix entries with more decimal places in scientific notation.

Definition at line 85 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 DenseMatrix< T >::resize (const unsigned intm, const unsigned intn) [inline]Resize the matrix. Will never free memory, but may allocate more. Sets all elements to 0.

Definition at line 474 of file dense_matrix.h.

References libMesh::zero.

Referenced by HPCoarsenTest::add_projection(), 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(), System::ProjectVector::operator()(), FEMContext::reinit(), and HPCoarsenTest::select_refinement().

{
  _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) [virtual]Right multiplies by the matrix M3.

Implements DenseMatrixBase< T >.

Definition at line 139 of file dense_matrix.C.

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

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

{
  if (this->use_blas)
    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)Right multiplies by the transpose of the matrix A

Definition at line 165 of file dense_matrix.C.

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

{
  if (this->use_blas)
    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 Tfactor) [inline]Multiplies every element in the matrix by factor.

Definition at line 549 of file dense_matrix.h.

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

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

Referenced by FEMSystem::assembly(), EulerSolver::element_residual(), Euler2Solver::element_residual(), EulerSolver::side_residual(), and Euler2Solver::side_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 DenseMatrix< T >::transpose (const unsigned inti, const unsigned intj) const [inline]Returns:

: the (i,j) element of the transposed matrix.

Definition at line 725 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) constPerform matrix vector multiplication.

Definition at line 227 of file dense_matrix.C.

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

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

  // Make sure the sizes are compatible
  libmesh_assert(n_cols == arg.size());
  libmesh_assert(n_rows == dest.size());

  dest.zero();
  DenseMatrix<T> A(*this);

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

template<typename T> void DenseMatrix< T >::vector_mult_add (DenseVector< T > &dest, const Tfactor, const DenseVector< T > &arg) constPerform matrix vector multiplication and add scaled result to dest.

Definition at line 246 of file dense_matrix.C.

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

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

Referenced by HPCoarsenTest::add_projection(), FEMSystem::assembly(), FEBase::coarsened_dof_values(), EulerSolver::element_residual(), Euler2Solver::element_residual(), System::ProjectVector::operator()(), HPCoarsenTest::select_refinement(), EulerSolver::side_residual(), and Euler2Solver::side_residual().

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

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

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 DenseMatrixBase< Number >::m(), DenseMatrix< T >::operator=(), 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 DenseMatrixBase< Number >::n(), DenseMatrix< T >::operator=(), and DenseMatrix< T >::swap().

template<typename T> std::vector<T> DenseMatrix< T >::_val [private]The actual data values, stored as a 1D array.

Definition at line 321 of file dense_matrix.h.

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