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 @@ -318,16 +318,74 @@ class CholmodBase : public SparseSolverB * \returns a reference to \c *this. */ Derived& setShift(const RealScalar& offset) { m_shiftOffset[0] = offset; return derived(); } + /** Efficiently obtains the determinant of the factor + * + * \returns the determinant of the factor + */ + RealScalar determinant() const { + return std::exp(logDeterminant()); + } + + /** Efficiently obtains the log determinant of the factor + * + * \returns the log determinant of the factor + */ + RealScalar logDeterminant() const { + eigen_assert((internal::is_same::value || internal::is_same::value) && + "Determinant is not implemented for complex matrices"); + 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. + + int *super = (int*) cf->super; + int *pi = (int*) cf->pi; + int *px = (int*) cf->px; + for (size_t sx=0; sx < cf->nsuper; ++sx) { + int ncols = super[sx + 1] - super[sx]; + int nrows = pi[sx + 1] - pi[sx]; + for (int cx=px[sx]; cx < px[sx] + nrows * ncols; cx += nrows+1) { + logDet += std::log(x[cx]); + } + } + } 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) + int *p = (int*) cf->p; + for (size_t ex=0; ex < cf->n; ++ex) { + logDet += std::log( 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 @@ -7,17 +7,17 @@ // Public License v. 2.0. If a copy of the MPL was not distributed // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. #define EIGEN_NO_DEBUG_SMALL_PRODUCT_BLOCKS #include "sparse_solver.h" #include -template void test_cholmod_T() +template void test1_cholmod_T() { CholmodDecomposition, Lower> g_chol_colmajor_lower; g_chol_colmajor_lower.setMode(CholmodSupernodalLLt); CholmodDecomposition, Upper> g_chol_colmajor_upper; g_chol_colmajor_upper.setMode(CholmodSupernodalLLt); CholmodDecomposition, Lower> g_llt_colmajor_lower; g_llt_colmajor_lower.setMode(CholmodSimplicialLLt); CholmodDecomposition, Upper> g_llt_colmajor_upper; g_llt_colmajor_upper.setMode(CholmodSimplicialLLt); CholmodDecomposition, Lower> g_ldlt_colmajor_lower; g_ldlt_colmajor_lower.setMode(CholmodLDLt); CholmodDecomposition, Upper> g_ldlt_colmajor_upper; g_ldlt_colmajor_upper.setMode(CholmodLDLt); @@ -36,22 +36,34 @@ 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); +} + +template void test2_cholmod_T() +{ + test1_cholmod_T(); + + CholmodSupernodalLLT, Lower> chol_colmajor_lower; + CholmodSupernodalLLT, Upper> chol_colmajor_upper; + CholmodSimplicialLLT, Lower> llt_colmajor_lower; + CholmodSimplicialLLT, Upper> llt_colmajor_upper; + CholmodSimplicialLDLT, Lower> ldlt_colmajor_lower; + CholmodSimplicialLDLT, Upper> 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 >()); + CALL_SUBTEST_1(test2_cholmod_T()); + CALL_SUBTEST_2(test1_cholmod_T >()); }