MOAB  4.9.3pre
Eigen::SimplicialCholeskyBase< Derived > Class Template Reference

A base class for direct sparse Cholesky factorizations. More...

#include <SimplicialCholesky.h>

Inheritance diagram for Eigen::SimplicialCholeskyBase< Derived >:
Collaboration diagram for Eigen::SimplicialCholeskyBase< Derived >:

List of all members.

Classes

struct  keep_diag

Public Types

enum  { UpLo = internal::traits<Derived>::UpLo }
enum  { ColsAtCompileTime = MatrixType::ColsAtCompileTime, MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime }
typedef internal::traits
< Derived >::MatrixType 
MatrixType
typedef internal::traits
< Derived >::OrderingType 
OrderingType
typedef MatrixType::Scalar Scalar
typedef MatrixType::RealScalar RealScalar
typedef MatrixType::StorageIndex StorageIndex
typedef SparseMatrix< Scalar,
ColMajor, StorageIndex
CholMatrixType
typedef CholMatrixType const * ConstCholMatrixPtr
typedef Matrix< Scalar,
Dynamic, 1 > 
VectorType
typedef Matrix< StorageIndex,
Dynamic, 1 > 
VectorI

Public Member Functions

 SimplicialCholeskyBase ()
 SimplicialCholeskyBase (const MatrixType &matrix)
 ~SimplicialCholeskyBase ()
Derived & derived ()
const Derived & derived () const
Index cols () const
Index rows () const
ComputationInfo info () const
 Reports whether previous computation was successful.
const PermutationMatrix
< Dynamic, Dynamic,
StorageIndex > & 
permutationP () const
const PermutationMatrix
< Dynamic, Dynamic,
StorageIndex > & 
permutationPinv () const
Derived & setShift (const RealScalar &offset, const RealScalar &scale=1)
template<typename Stream >
void dumpMemory (Stream &s)
template<typename Rhs , typename Dest >
void _solve_impl (const MatrixBase< Rhs > &b, MatrixBase< Dest > &dest) const
template<typename Rhs , typename Dest >
void _solve_impl (const SparseMatrixBase< Rhs > &b, SparseMatrixBase< Dest > &dest) const

Protected Member Functions

template<bool DoLDLT>
void compute (const MatrixType &matrix)
template<bool DoLDLT>
void factorize (const MatrixType &a)
template<bool DoLDLT>
void factorize_preordered (const CholMatrixType &a)
void analyzePattern (const MatrixType &a, bool doLDLT)
void analyzePattern_preordered (const CholMatrixType &a, bool doLDLT)
void ordering (const MatrixType &a, ConstCholMatrixPtr &pmat, CholMatrixType &ap)

Protected Attributes

ComputationInfo m_info
bool m_factorizationIsOk
bool m_analysisIsOk
CholMatrixType m_matrix
VectorType m_diag
VectorI m_parent
VectorI m_nonZerosPerCol
PermutationMatrix< Dynamic,
Dynamic, StorageIndex
m_P
PermutationMatrix< Dynamic,
Dynamic, StorageIndex
m_Pinv
RealScalar m_shiftOffset
RealScalar m_shiftScale

Private Types

typedef SparseSolverBase< Derived > Base

Detailed Description

template<typename Derived>
class Eigen::SimplicialCholeskyBase< Derived >

A base class for direct sparse Cholesky factorizations.

This is a base class for LL^T and LDL^T Cholesky factorizations of sparse matrices that are selfadjoint and positive definite. These factorizations allow for solving A.X = B where X and B can be either dense or sparse.

In order to reduce the fill-in, a symmetric permutation P is applied prior to the factorization such that the factorized matrix is P A P^-1.

Template Parameters:
Derivedthe type of the derived class, that is the actual factorization type.

Definition at line 55 of file SimplicialCholesky.h.


Member Typedef Documentation

template<typename Derived>
typedef CholMatrixType const* Eigen::SimplicialCholeskyBase< Derived >::ConstCholMatrixPtr

Definition at line 68 of file SimplicialCholesky.h.

template<typename Derived>
typedef internal::traits<Derived>::OrderingType Eigen::SimplicialCholeskyBase< Derived >::OrderingType

Definition at line 62 of file SimplicialCholesky.h.

template<typename Derived>
typedef Matrix<StorageIndex,Dynamic,1> Eigen::SimplicialCholeskyBase< Derived >::VectorI

Definition at line 70 of file SimplicialCholesky.h.


Member Enumeration Documentation

template<typename Derived>
anonymous enum
Enumerator:
UpLo 

Definition at line 63 of file SimplicialCholesky.h.

{ UpLo = internal::traits<Derived>::UpLo };
template<typename Derived>
anonymous enum
Enumerator:
ColsAtCompileTime 
MaxColsAtCompileTime 

Definition at line 72 of file SimplicialCholesky.h.

         {
      ColsAtCompileTime = MatrixType::ColsAtCompileTime,
      MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
    };

Constructor & Destructor Documentation

template<typename Derived>
Eigen::SimplicialCholeskyBase< Derived >::SimplicialCholeskyBase ( ) [inline]

Default constructor

Definition at line 82 of file SimplicialCholesky.h.

template<typename Derived>
Eigen::SimplicialCholeskyBase< Derived >::SimplicialCholeskyBase ( const MatrixType matrix) [inline, explicit]

Definition at line 86 of file SimplicialCholesky.h.

      : m_info(Success), m_shiftOffset(0), m_shiftScale(1)
    {
      derived().compute(matrix);
    }
template<typename Derived>
Eigen::SimplicialCholeskyBase< Derived >::~SimplicialCholeskyBase ( ) [inline]

Definition at line 92 of file SimplicialCholesky.h.

    {
    }

Member Function Documentation

template<typename Derived>
template<typename Rhs , typename Dest >
void Eigen::SimplicialCholeskyBase< Derived >::_solve_impl ( const MatrixBase< Rhs > &  b,
MatrixBase< Dest > &  dest 
) const [inline]

Reimplemented in Eigen::SimplicialCholesky< _MatrixType, _UpLo, _Ordering >.

Definition at line 156 of file SimplicialCholesky.h.

    {
      eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
      eigen_assert(m_matrix.rows()==b.rows());

      if(m_info!=Success)
        return;

      if(m_P.size()>0)
        dest = m_P * b;
      else
        dest = b;

      if(m_matrix.nonZeros()>0) // otherwise L==I
        derived().matrixL().solveInPlace(dest);

      if(m_diag.size()>0)
        dest = m_diag.asDiagonal().inverse() * dest;

      if (m_matrix.nonZeros()>0) // otherwise U==I
        derived().matrixU().solveInPlace(dest);

      if(m_P.size()>0)
        dest = m_Pinv * dest;
    }
template<typename Derived>
template<typename Rhs , typename Dest >
void Eigen::SimplicialCholeskyBase< Derived >::_solve_impl ( const SparseMatrixBase< Rhs > &  b,
SparseMatrixBase< Dest > &  dest 
) const [inline]

default implementation of solving with a sparse rhs

Reimplemented from Eigen::SparseSolverBase< Derived >.

Reimplemented in Eigen::SimplicialCholesky< _MatrixType, _UpLo, _Ordering >.

Definition at line 183 of file SimplicialCholesky.h.

template<typename Derived>
void Eigen::SimplicialCholeskyBase< Derived >::analyzePattern ( const MatrixType a,
bool  doLDLT 
) [inline, protected]

Definition at line 230 of file SimplicialCholesky.h.

    {
      eigen_assert(a.rows()==a.cols());
      Index size = a.cols();
      CholMatrixType tmp(size,size);
      ConstCholMatrixPtr pmat;
      ordering(a, pmat, tmp);
      analyzePattern_preordered(*pmat,doLDLT);
    }
template<typename Derived >
void Eigen::SimplicialCholeskyBase< Derived >::analyzePattern_preordered ( const CholMatrixType a,
bool  doLDLT 
) [protected]

Definition at line 51 of file SimplicialCholesky_impl.h.

{
  const StorageIndex size = StorageIndex(ap.rows());
  m_matrix.resize(size, size);
  m_parent.resize(size);
  m_nonZerosPerCol.resize(size);

  ei_declare_aligned_stack_constructed_variable(StorageIndex, tags, size, 0);

  for(StorageIndex k = 0; k < size; ++k)
  {
    /* L(k,:) pattern: all nodes reachable in etree from nz in A(0:k-1,k) */
    m_parent[k] = -1;             /* parent of k is not yet known */
    tags[k] = k;                  /* mark node k as visited */
    m_nonZerosPerCol[k] = 0;      /* count of nonzeros in column k of L */
    for(typename CholMatrixType::InnerIterator it(ap,k); it; ++it)
    {
      StorageIndex i = it.index();
      if(i < k)
      {
        /* follow path from i to root of etree, stop at flagged node */
        for(; tags[i] != k; i = m_parent[i])
        {
          /* find parent of i if not yet determined */
          if (m_parent[i] == -1)
            m_parent[i] = k;
          m_nonZerosPerCol[i]++;        /* L (k,i) is nonzero */
          tags[i] = k;                  /* mark i as visited */
        }
      }
    }
  }

  /* construct Lp index array from m_nonZerosPerCol column counts */
  StorageIndex* Lp = m_matrix.outerIndexPtr();
  Lp[0] = 0;
  for(StorageIndex k = 0; k < size; ++k)
    Lp[k+1] = Lp[k] + m_nonZerosPerCol[k] + (doLDLT ? 0 : 1);

  m_matrix.resizeNonZeros(Lp[size]);

  m_isInitialized     = true;
  m_info              = Success;
  m_analysisIsOk      = true;
  m_factorizationIsOk = false;
}
template<typename Derived>
Index Eigen::SimplicialCholeskyBase< Derived >::cols ( void  ) const [inline]

Definition at line 99 of file SimplicialCholesky.h.

{ return m_matrix.cols(); }
template<typename Derived>
template<bool DoLDLT>
void Eigen::SimplicialCholeskyBase< Derived >::compute ( const MatrixType matrix) [inline, protected]

Computes the sparse Cholesky decomposition of matrix

Reimplemented in Eigen::SimplicialCholesky< _MatrixType, _UpLo, _Ordering >, Eigen::SimplicialLDLT< _MatrixType, _UpLo, _Ordering >, and Eigen::SimplicialLLT< _MatrixType, _UpLo, _Ordering >.

Definition at line 194 of file SimplicialCholesky.h.

    {
      eigen_assert(matrix.rows()==matrix.cols());
      Index size = matrix.cols();
      CholMatrixType tmp(size,size);
      ConstCholMatrixPtr pmat;
      ordering(matrix, pmat, tmp);
      analyzePattern_preordered(*pmat, DoLDLT);
      factorize_preordered<DoLDLT>(*pmat);
    }
template<typename Derived>
Derived& Eigen::SimplicialCholeskyBase< Derived >::derived ( ) [inline]

Reimplemented from Eigen::SparseSolverBase< Derived >.

Definition at line 96 of file SimplicialCholesky.h.

{ return *static_cast<Derived*>(this); }
template<typename Derived>
const Derived& Eigen::SimplicialCholeskyBase< Derived >::derived ( ) const [inline]

Reimplemented from Eigen::SparseSolverBase< Derived >.

Definition at line 97 of file SimplicialCholesky.h.

{ return *static_cast<const Derived*>(this); }
template<typename Derived>
template<typename Stream >
void Eigen::SimplicialCholeskyBase< Derived >::dumpMemory ( Stream &  s) [inline]

Definition at line 142 of file SimplicialCholesky.h.

    {
      int total = 0;
      s << "  L:        " << ((total+=(m_matrix.cols()+1) * sizeof(int) + m_matrix.nonZeros()*(sizeof(int)+sizeof(Scalar))) >> 20) << "Mb" << "\n";
      s << "  diag:     " << ((total+=m_diag.size() * sizeof(Scalar)) >> 20) << "Mb" << "\n";
      s << "  tree:     " << ((total+=m_parent.size() * sizeof(int)) >> 20) << "Mb" << "\n";
      s << "  nonzeros: " << ((total+=m_nonZerosPerCol.size() * sizeof(int)) >> 20) << "Mb" << "\n";
      s << "  perm:     " << ((total+=m_P.size() * sizeof(int)) >> 20) << "Mb" << "\n";
      s << "  perm^-1:  " << ((total+=m_Pinv.size() * sizeof(int)) >> 20) << "Mb" << "\n";
      s << "  TOTAL:    " << (total>> 20) << "Mb" << "\n";
    }
template<typename Derived>
template<bool DoLDLT>
void Eigen::SimplicialCholeskyBase< Derived >::factorize ( const MatrixType a) [inline, protected]

Reimplemented in Eigen::SimplicialCholesky< _MatrixType, _UpLo, _Ordering >, Eigen::SimplicialLDLT< _MatrixType, _UpLo, _Ordering >, and Eigen::SimplicialLLT< _MatrixType, _UpLo, _Ordering >.

Definition at line 206 of file SimplicialCholesky.h.

    {
      eigen_assert(a.rows()==a.cols());
      Index size = a.cols();
      CholMatrixType tmp(size,size);
      ConstCholMatrixPtr pmat;
      
      if(m_P.size()==0 && (UpLo&Upper)==Upper)
      {
        // If there is no ordering, try to directly use the input matrix without any copy
        internal::simplicial_cholesky_grab_input<CholMatrixType,MatrixType>::run(a, pmat, tmp);
      }
      else
      {
        tmp.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
        pmat = &tmp;
      }
      
      factorize_preordered<DoLDLT>(*pmat);
    }
template<typename Derived >
template<bool DoLDLT>
void Eigen::SimplicialCholeskyBase< Derived >::factorize_preordered ( const CholMatrixType a) [protected]

Definition at line 101 of file SimplicialCholesky_impl.h.

{
  using std::sqrt;

  eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
  eigen_assert(ap.rows()==ap.cols());
  eigen_assert(m_parent.size()==ap.rows());
  eigen_assert(m_nonZerosPerCol.size()==ap.rows());

  const StorageIndex size = StorageIndex(ap.rows());
  const StorageIndex* Lp = m_matrix.outerIndexPtr();
  StorageIndex* Li = m_matrix.innerIndexPtr();
  Scalar* Lx = m_matrix.valuePtr();

  ei_declare_aligned_stack_constructed_variable(Scalar, y, size, 0);
  ei_declare_aligned_stack_constructed_variable(StorageIndex,  pattern, size, 0);
  ei_declare_aligned_stack_constructed_variable(StorageIndex,  tags, size, 0);

  bool ok = true;
  m_diag.resize(DoLDLT ? size : 0);

  for(StorageIndex k = 0; k < size; ++k)
  {
    // compute nonzero pattern of kth row of L, in topological order
    y[k] = 0.0;                     // Y(0:k) is now all zero
    StorageIndex top = size;               // stack for pattern is empty
    tags[k] = k;                    // mark node k as visited
    m_nonZerosPerCol[k] = 0;        // count of nonzeros in column k of L
    for(typename CholMatrixType::InnerIterator it(ap,k); it; ++it)
    {
      StorageIndex i = it.index();
      if(i <= k)
      {
        y[i] += numext::conj(it.value());            /* scatter A(i,k) into Y (sum duplicates) */
        Index len;
        for(len = 0; tags[i] != k; i = m_parent[i])
        {
          pattern[len++] = i;     /* L(k,i) is nonzero */
          tags[i] = k;            /* mark i as visited */
        }
        while(len > 0)
          pattern[--top] = pattern[--len];
      }
    }

    /* compute numerical values kth row of L (a sparse triangular solve) */

    RealScalar d = numext::real(y[k]) * m_shiftScale + m_shiftOffset;    // get D(k,k), apply the shift function, and clear Y(k)
    y[k] = 0.0;
    for(; top < size; ++top)
    {
      Index i = pattern[top];       /* pattern[top:n-1] is pattern of L(:,k) */
      Scalar yi = y[i];             /* get and clear Y(i) */
      y[i] = 0.0;

      /* the nonzero entry L(k,i) */
      Scalar l_ki;
      if(DoLDLT)
        l_ki = yi / m_diag[i];
      else
        yi = l_ki = yi / Lx[Lp[i]];

      Index p2 = Lp[i] + m_nonZerosPerCol[i];
      Index p;
      for(p = Lp[i] + (DoLDLT ? 0 : 1); p < p2; ++p)
        y[Li[p]] -= numext::conj(Lx[p]) * yi;
      d -= numext::real(l_ki * numext::conj(yi));
      Li[p] = k;                          /* store L(k,i) in column form of L */
      Lx[p] = l_ki;
      ++m_nonZerosPerCol[i];              /* increment count of nonzeros in col i */
    }
    if(DoLDLT)
    {
      m_diag[k] = d;
      if(d == RealScalar(0))
      {
        ok = false;                         /* failure, D(k,k) is zero */
        break;
      }
    }
    else
    {
      Index p = Lp[k] + m_nonZerosPerCol[k]++;
      Li[p] = k ;                /* store L(k,k) = sqrt (d) in column k */
      if(d <= RealScalar(0)) {
        ok = false;              /* failure, matrix is not positive definite */
        break;
      }
      Lx[p] = sqrt(d) ;
    }
  }

  m_info = ok ? Success : NumericalIssue;
  m_factorizationIsOk = true;
}
template<typename Derived>
ComputationInfo Eigen::SimplicialCholeskyBase< Derived >::info ( ) const [inline]

Reports whether previous computation was successful.

Returns:
Success if computation was succesful, NumericalIssue if the matrix.appears to be negative.

Definition at line 107 of file SimplicialCholesky.h.

    {
      eigen_assert(m_isInitialized && "Decomposition is not initialized.");
      return m_info;
    }
template<typename Derived >
void Eigen::SimplicialCholeskyBase< Derived >::ordering ( const MatrixType a,
ConstCholMatrixPtr pmat,
CholMatrixType ap 
) [protected]

Definition at line 650 of file SimplicialCholesky.h.

{
  eigen_assert(a.rows()==a.cols());
  const Index size = a.rows();
  pmat = &ap;
  // Note that ordering methods compute the inverse permutation
  if(!internal::is_same<OrderingType,NaturalOrdering<Index> >::value)
  {
    {
      CholMatrixType C;
      C = a.template selfadjointView<UpLo>();
      
      OrderingType ordering;
      ordering(C,m_Pinv);
    }

    if(m_Pinv.size()>0) m_P = m_Pinv.inverse();
    else                m_P.resize(0);
    
    ap.resize(size,size);
    ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
  }
  else
  {
    m_Pinv.resize(0);
    m_P.resize(0);
    if(int(UpLo)==int(Lower) || MatrixType::IsRowMajor)
    {
      // we have to transpose the lower part to to the upper one
      ap.resize(size,size);
      ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>();
    }
    else
      internal::simplicial_cholesky_grab_input<CholMatrixType,MatrixType>::run(a, pmat, ap);
  }  
}
template<typename Derived>
const PermutationMatrix<Dynamic,Dynamic,StorageIndex>& Eigen::SimplicialCholeskyBase< Derived >::permutationP ( ) const [inline]
Returns:
the permutation P
See also:
permutationPinv()

Definition at line 115 of file SimplicialCholesky.h.

    { return m_P; }
template<typename Derived>
const PermutationMatrix<Dynamic,Dynamic,StorageIndex>& Eigen::SimplicialCholeskyBase< Derived >::permutationPinv ( ) const [inline]
Returns:
the inverse P^-1 of the permutation P
See also:
permutationP()

Definition at line 120 of file SimplicialCholesky.h.

    { return m_Pinv; }
template<typename Derived>
Index Eigen::SimplicialCholeskyBase< Derived >::rows ( void  ) const [inline]

Definition at line 100 of file SimplicialCholesky.h.

{ return m_matrix.rows(); }
template<typename Derived>
Derived& Eigen::SimplicialCholeskyBase< Derived >::setShift ( const RealScalar offset,
const RealScalar scale = 1 
) [inline]

Sets the shift parameters that will be used to adjust the diagonal coefficients during the numerical factorization.

During the numerical factorization, the diagonal coefficients are transformed by the following linear model:
d_ii = offset + scale * d_ii

The default is the identity transformation with offset=0, and scale=1.

Returns:
a reference to *this.

Definition at line 132 of file SimplicialCholesky.h.

    {
      m_shiftOffset = offset;
      m_shiftScale = scale;
      return derived();
    }

Member Data Documentation

template<typename Derived>
bool Eigen::SimplicialCholeskyBase< Derived >::m_analysisIsOk [protected]

Definition at line 253 of file SimplicialCholesky.h.

template<typename Derived>
VectorType Eigen::SimplicialCholeskyBase< Derived >::m_diag [protected]

Definition at line 256 of file SimplicialCholesky.h.

template<typename Derived>
bool Eigen::SimplicialCholeskyBase< Derived >::m_factorizationIsOk [protected]

Definition at line 252 of file SimplicialCholesky.h.

template<typename Derived>
ComputationInfo Eigen::SimplicialCholeskyBase< Derived >::m_info [mutable, protected]

Definition at line 251 of file SimplicialCholesky.h.

template<typename Derived>
CholMatrixType Eigen::SimplicialCholeskyBase< Derived >::m_matrix [protected]

Definition at line 255 of file SimplicialCholesky.h.

template<typename Derived>
VectorI Eigen::SimplicialCholeskyBase< Derived >::m_nonZerosPerCol [protected]

Definition at line 258 of file SimplicialCholesky.h.

template<typename Derived>
PermutationMatrix<Dynamic,Dynamic,StorageIndex> Eigen::SimplicialCholeskyBase< Derived >::m_P [protected]

Definition at line 259 of file SimplicialCholesky.h.

template<typename Derived>
VectorI Eigen::SimplicialCholeskyBase< Derived >::m_parent [protected]

Definition at line 257 of file SimplicialCholesky.h.

template<typename Derived>
PermutationMatrix<Dynamic,Dynamic,StorageIndex> Eigen::SimplicialCholeskyBase< Derived >::m_Pinv [protected]

Definition at line 260 of file SimplicialCholesky.h.

template<typename Derived>
RealScalar Eigen::SimplicialCholeskyBase< Derived >::m_shiftOffset [protected]

Definition at line 262 of file SimplicialCholesky.h.

template<typename Derived>
RealScalar Eigen::SimplicialCholeskyBase< Derived >::m_shiftScale [protected]

Definition at line 263 of file SimplicialCholesky.h.


The documentation for this class was generated from the following files:
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines