Bugzilla – Attachment 644 Details for
Bug 1095
log det for cholmod
Home
|
New
|
Browse
|
Search
|
[?]
|
Reports
|
Requests
|
Help
|
Log In
[x]
|
Forgot Password
Login:
[x]
This bugzilla service is closed. All entries have been migrated to
https://gitlab.com/libeigen/eigen
[patch]
add log det
z (text/plain), 4.48 KB, created by
Joshua Pritikin
on 2016-01-13 16:58:50 UTC
(
hide
)
Description:
add log det
Filename:
MIME Type:
Creator:
Joshua Pritikin
Created:
2016-01-13 16:58:50 UTC
Size:
4.48 KB
patch
obsolete
>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<const Array<Scalar,1,Dynamic>, 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<typename Stream> > 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<typename T> 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<double>()); > CALL_SUBTEST_2(test_cholmod_T<std::complex<double> >()); > }
You cannot view the attachment while viewing its details because your browser does not support IFRAMEs.
View the attachment on a separate page
.
View Attachment As Diff
View Attachment As Raw
Actions:
View
|
Diff
Attachments on
bug 1095
:
636
|
637
|
639
|
640
| 644