Bugzilla – Attachment 637 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]
proposed change
z.patch (text/plain), 5.73 KB, created by
Joshua Pritikin
on 2015-12-27 20:26:32 UTC
(
hide
)
Description:
proposed change
Filename:
MIME Type:
Creator:
Joshua Pritikin
Created:
2015-12-27 20:26:32 UTC
Size:
5.73 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 >@@ -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<Scalar,float>::value || internal::is_same<Scalar,double>::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<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 >@@ -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 <Eigen/CholmodSupport> > >-template<typename T> void test_cholmod_T() >+template<typename T> void test1_cholmod_T() > { > CholmodDecomposition<SparseMatrix<T>, Lower> g_chol_colmajor_lower; g_chol_colmajor_lower.setMode(CholmodSupernodalLLt); > CholmodDecomposition<SparseMatrix<T>, Upper> g_chol_colmajor_upper; g_chol_colmajor_upper.setMode(CholmodSupernodalLLt); > CholmodDecomposition<SparseMatrix<T>, Lower> g_llt_colmajor_lower; g_llt_colmajor_lower.setMode(CholmodSimplicialLLt); > CholmodDecomposition<SparseMatrix<T>, Upper> g_llt_colmajor_upper; g_llt_colmajor_upper.setMode(CholmodSimplicialLLt); > CholmodDecomposition<SparseMatrix<T>, Lower> g_ldlt_colmajor_lower; g_ldlt_colmajor_lower.setMode(CholmodLDLt); > CholmodDecomposition<SparseMatrix<T>, Upper> g_ldlt_colmajor_upper; g_ldlt_colmajor_upper.setMode(CholmodLDLt); > >@@ -36,22 +36,34 @@ 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); >+} >+ >+template<typename T> void test2_cholmod_T() >+{ >+ test1_cholmod_T<T>(); >+ >+ CholmodSupernodalLLT<SparseMatrix<T>, Lower> chol_colmajor_lower; >+ CholmodSupernodalLLT<SparseMatrix<T>, Upper> chol_colmajor_upper; >+ CholmodSimplicialLLT<SparseMatrix<T>, Lower> llt_colmajor_lower; >+ CholmodSimplicialLLT<SparseMatrix<T>, Upper> llt_colmajor_upper; >+ CholmodSimplicialLDLT<SparseMatrix<T>, Lower> ldlt_colmajor_lower; >+ CholmodSimplicialLDLT<SparseMatrix<T>, 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<double>()); >- CALL_SUBTEST_2(test_cholmod_T<std::complex<double> >()); >+ CALL_SUBTEST_1(test2_cholmod_T<double>()); >+ CALL_SUBTEST_2(test1_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