diff --git a/Eigen/src/CholmodSupport/CholmodSupport.h b/Eigen/src/CholmodSupport/CholmodSupport.h --- a/Eigen/src/CholmodSupport/CholmodSupport.h +++ b/Eigen/src/CholmodSupport/CholmodSupport.h @@ -73,17 +73,17 @@ cholmod_sparse viewAsCholmod(SparseMatri res.dtype = 0; res.stype = -1; if (internal::is_same<_StorageIndex,int>::value) { res.itype = CHOLMOD_INT; } - else if (internal::is_same<_StorageIndex,SuiteSparse_long>::value) + else if (internal::is_same<_StorageIndex,long>::value) { res.itype = CHOLMOD_LONG; } else { eigen_assert(false && "Index type not supported yet"); } @@ -318,16 +318,70 @@ class CholmodBase : public SparseSolverB * \returns a reference to \c *this. */ Derived& setShift(const RealScalar& offset) { m_shiftOffset[0] = offset; return derived(); } + /** \returns the determinant of the underlying matrix from the current factorization */ + Scalar determinant() const { + using std::exp; + return exp(logDeterminant()); + } + + /** \returns the log determinant of the underlying matrix from the current factorization */ + Scalar logDeterminant() const { + using std::log; + using numext::real; + eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()"); + + // Based on https://github.com/njsmith/scikits-sparse/blob/master/scikits/sparse/cholmod.pyx + cholmod_factor *cf = m_cholmodFactor; + RealScalar logDet = 0; + Scalar *x = (Scalar*) cf->x; + if (cf->is_super) { + // This is a supernodal factorization, which is stored as a bunch + // of dense, lower-triangular, column-major arrays packed into the + // x vector. This is not documented in the CHOLMOD user-guide, or + // anywhere else as far as I can tell; I got the details from + // CVXOPT's C/cholmod.c. + + StorageIndex *super = (StorageIndex*) cf->super; + StorageIndex *pi = (StorageIndex*) cf->pi; + StorageIndex *px = (StorageIndex*) cf->px; + for (size_t sx=0; sx < cf->nsuper; ++sx) { + StorageIndex ncols = super[sx + 1] - super[sx]; + StorageIndex nrows = pi[sx + 1] - pi[sx]; + + Map, 0, InnerStride<> > + s1(x + px[sx], ncols, InnerStride<>(nrows+1)); + logDet += s1.real().log().sum(); + } + } else { + // This is a simplicial factorization, which is simply stored as a + // sparse CSC matrix in x, p, i. We want the diagonal, which is + // just the first entry in each column; p gives the offsets in x to + // the beginning of each column. + // + // The ->p array actually has n+1 entries, but only the first n + // entries actually point to real columns (the last entry is a + // sentinel) + StorageIndex *p = (StorageIndex*) cf->p; + for (size_t ex=0; ex < cf->n; ++ex) { + logDet += log(real( x[p[ex]] )); + } + } + if (cf->is_ll) { + logDet *= 2.0; + } + return logDet; + }; + template void dumpMemory(Stream& /*s*/) {} protected: mutable cholmod_common m_cholmod; cholmod_factor* m_cholmodFactor; RealScalar m_shiftOffset[2]; diff --git a/test/cholmod_support.cpp b/test/cholmod_support.cpp --- a/test/cholmod_support.cpp +++ b/test/cholmod_support.cpp @@ -36,22 +36,22 @@ template void test_cholmod_T check_sparse_spd_solving(g_ldlt_colmajor_upper); check_sparse_spd_solving(chol_colmajor_lower); check_sparse_spd_solving(chol_colmajor_upper); check_sparse_spd_solving(llt_colmajor_lower); check_sparse_spd_solving(llt_colmajor_upper); check_sparse_spd_solving(ldlt_colmajor_lower); check_sparse_spd_solving(ldlt_colmajor_upper); - -// check_sparse_spd_determinant(chol_colmajor_lower); -// check_sparse_spd_determinant(chol_colmajor_upper); -// check_sparse_spd_determinant(llt_colmajor_lower); -// check_sparse_spd_determinant(llt_colmajor_upper); -// check_sparse_spd_determinant(ldlt_colmajor_lower); -// check_sparse_spd_determinant(ldlt_colmajor_upper); + + check_sparse_spd_determinant(chol_colmajor_lower); + check_sparse_spd_determinant(chol_colmajor_upper); + check_sparse_spd_determinant(llt_colmajor_lower); + check_sparse_spd_determinant(llt_colmajor_upper); + check_sparse_spd_determinant(ldlt_colmajor_lower); + check_sparse_spd_determinant(ldlt_colmajor_upper); } void test_cholmod_support() { CALL_SUBTEST_1(test_cholmod_T()); CALL_SUBTEST_2(test_cholmod_T >()); }