Attachment 639 Details for Bug 1095 – implement complex and address feedback

This bugzilla service is closed. All entries have been migrated to https://gitlab.com/libeigen/eigen

[patch] implement complex and address feedback

z.patch (text/plain), 3.61 KB, created by Joshua Pritikin on 2015-12-29 14:19:18 UTC

(hide)

Description:

Filename:

MIME Type:

Creator: Joshua Pritikin

Created: 2015-12-29 14:19:18 UTC

Size: 3.61 KB

patch

obsolete

>diff --git a/Eigen/src/CholmodSupport/CholmodSupport.h b/Eigen/src/CholmodSupport/CholmodSupport.h
>@@ -318,16 +318,69 @@ 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 */
>+    RealScalar determinant() const {
>+      using std::exp;
>+      return exp(logDeterminant());
>+    }
>+
>+    /** \returns the log determinant of the underlying matrix from the current factorization */
>+    RealScalar 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];
>+	  for (StorageIndex cx=px[sx]; cx < px[sx] + nrows * ncols; cx += nrows+1) {
>+	    logDet += log(real(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)
>+	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,15 @@ 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);
> }
> 
> void test_cholmod_support()
> {
>   CALL_SUBTEST_1(test_cholmod_T<double>());
>   CALL_SUBTEST_2(test_cholmod_T<std::complex<double> >());
> }

Actions: View | Diff

Attachments on bug 1095: 636 | 637 | 639 | 640 | 644