Logo Search packages:      
Sourcecode: eigen3 version File versions  Download package

template<typename _MatrixType , int _UpLo>
void SimplicialCholesky< _MatrixType, _UpLo >::analyzePattern ( const MatrixType &  a)
Returns:
the solution x of $ A x = b $ using the current decomposition of A.
See also:
compute() Performs a symbolic decomposition on the sparcity of matrix.

This function is particularly useful when solving for several problems having the same structure.

See also:
factorize()

Definition at line 290 of file SimplicialCholesky.h.

References internal::minimum_degree_ordering(), and SparseMatrix< _Scalar, _Options, _Index >::prune().

Referenced by SimplicialCholesky< _MatrixType, _UpLo >::compute().

{
  eigen_assert(a.rows()==a.cols());
  const Index size = a.rows();
  m_matrix.resize(size, size);
  m_parent.resize(size);
  m_nonZerosPerCol.resize(size);
  
  ei_declare_aligned_stack_constructed_variable(Index, tags, size, 0);
  
  // TODO allows to configure the permutation
  {
    CholMatrixType C;
    C = a.template selfadjointView<UpLo>();
    // remove diagonal entries:
    C.prune(keep_diag());
    internal::minimum_degree_ordering(C, m_P);
  }
  
  if(m_P.size()>0)
    m_Pinv  = m_P.inverse();
  else
    m_Pinv.resize(0);
  
  SparseMatrix<Scalar,ColMajor,Index> ap(size,size);
  ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_Pinv);
  
  for(Index 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)
    {
      Index 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 */
  Index* Lp = m_matrix._outerIndexPtr();
  Lp[0] = 0;
  for(Index k = 0; k < size; ++k)
    Lp[k+1] = Lp[k] + m_nonZerosPerCol[k] + (m_LDLt ? 0 : 1);

  m_matrix.resizeNonZeros(Lp[size]);
  
  m_isInitialized     = true;
  m_info              = Success;
  m_analysisIsOk      = true;
  m_factorizationIsOk = false;
}

Here is the call graph for this function:

Here is the caller graph for this function:


Generated by  Doxygen 1.6.0   Back to index