This bugzilla service is closed. All entries have been migrated to https://gitlab.com/libeigen/eigen
View | Details | Raw Unified | Return to bug 1095 | Differences between
and this patch

Collapse All | Expand All

(-)a/Eigen/src/CholmodSupport/CholmodSupport.h (+53 lines)
Lines 318-333 class CholmodBase : public SparseSolverB Link Here
318
      * \returns a reference to \c *this.
318
      * \returns a reference to \c *this.
319
      */
319
      */
320
    Derived& setShift(const RealScalar& offset)
320
    Derived& setShift(const RealScalar& offset)
321
    {
321
    {
322
      m_shiftOffset[0] = offset;
322
      m_shiftOffset[0] = offset;
323
      return derived();
323
      return derived();
324
    }
324
    }
325
    
325
    
326
    /** \returns the determinant of the underlying matrix from the current factorization */
327
    Scalar determinant() const {
328
      using std::exp;
329
      return exp(logDeterminant());
330
    }
331
332
    /** \returns the log determinant of the underlying matrix from the current factorization */
333
    Scalar logDeterminant() const {
334
      using std::log;
335
      using numext::real;
336
      eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
337
338
      // Based on https://github.com/njsmith/scikits-sparse/blob/master/scikits/sparse/cholmod.pyx
339
      cholmod_factor *cf = m_cholmodFactor;
340
      RealScalar logDet = 0;
341
      Scalar *x = (Scalar*) cf->x;
342
      if (cf->is_super) {
343
	// This is a supernodal factorization, which is stored as a bunch
344
	// of dense, lower-triangular, column-major arrays packed into the
345
	// x vector. This is not documented in the CHOLMOD user-guide, or
346
	// anywhere else as far as I can tell; I got the details from
347
	// CVXOPT's C/cholmod.c.
348
349
	StorageIndex *super = (StorageIndex*) cf->super;
350
	StorageIndex *pi = (StorageIndex*) cf->pi;
351
	StorageIndex *px = (StorageIndex*) cf->px;
352
	for (size_t sx=0; sx < cf->nsuper; ++sx) {
353
	  StorageIndex ncols = super[sx + 1] - super[sx];
354
	  StorageIndex nrows = pi[sx + 1] - pi[sx];
355
	  for (StorageIndex cx=px[sx]; cx < px[sx] + nrows * ncols; cx += nrows+1) {
356
	    logDet += log(real(x[cx]));
357
	  }
358
	}
359
      } else {
360
	// This is a simplicial factorization, which is simply stored as a
361
	// sparse CSC matrix in x, p, i. We want the diagonal, which is
362
	// just the first entry in each column; p gives the offsets in x to
363
	// the beginning of each column.
364
	//
365
	// The ->p array actually has n+1 entries, but only the first n
366
	// entries actually point to real columns (the last entry is a
367
	// sentinel)
368
	StorageIndex *p = (StorageIndex*) cf->p;
369
	for (size_t ex=0; ex < cf->n; ++ex) {
370
	  logDet += log(real( x[p[ex]] ));
371
	}
372
      }
373
      if (cf->is_ll) {
374
	logDet *= 2.0;
375
      }
376
      return logDet;
377
    };
378
326
    template<typename Stream>
379
    template<typename Stream>
327
    void dumpMemory(Stream& /*s*/)
380
    void dumpMemory(Stream& /*s*/)
328
    {}
381
    {}
329
    
382
    
330
  protected:
383
  protected:
331
    mutable cholmod_common m_cholmod;
384
    mutable cholmod_common m_cholmod;
332
    cholmod_factor* m_cholmodFactor;
385
    cholmod_factor* m_cholmodFactor;
333
    RealScalar m_shiftOffset[2];
386
    RealScalar m_shiftOffset[2];
(-)a/test/cholmod_support.cpp (-7 / +7 lines)
Lines 36-57 template<typename T> void test_cholmod_T Link Here
36
  check_sparse_spd_solving(g_ldlt_colmajor_upper);
36
  check_sparse_spd_solving(g_ldlt_colmajor_upper);
37
  
37
  
38
  check_sparse_spd_solving(chol_colmajor_lower);
38
  check_sparse_spd_solving(chol_colmajor_lower);
39
  check_sparse_spd_solving(chol_colmajor_upper);
39
  check_sparse_spd_solving(chol_colmajor_upper);
40
  check_sparse_spd_solving(llt_colmajor_lower);
40
  check_sparse_spd_solving(llt_colmajor_lower);
41
  check_sparse_spd_solving(llt_colmajor_upper);
41
  check_sparse_spd_solving(llt_colmajor_upper);
42
  check_sparse_spd_solving(ldlt_colmajor_lower);
42
  check_sparse_spd_solving(ldlt_colmajor_lower);
43
  check_sparse_spd_solving(ldlt_colmajor_upper);
43
  check_sparse_spd_solving(ldlt_colmajor_upper);
44
  
44
45
//  check_sparse_spd_determinant(chol_colmajor_lower);
45
  check_sparse_spd_determinant(chol_colmajor_lower);
46
//  check_sparse_spd_determinant(chol_colmajor_upper);
46
  check_sparse_spd_determinant(chol_colmajor_upper);
47
//  check_sparse_spd_determinant(llt_colmajor_lower);
47
  check_sparse_spd_determinant(llt_colmajor_lower);
48
//  check_sparse_spd_determinant(llt_colmajor_upper);
48
  check_sparse_spd_determinant(llt_colmajor_upper);
49
//  check_sparse_spd_determinant(ldlt_colmajor_lower);
49
  check_sparse_spd_determinant(ldlt_colmajor_lower);
50
//  check_sparse_spd_determinant(ldlt_colmajor_upper);
50
  check_sparse_spd_determinant(ldlt_colmajor_upper);
51
}
51
}
52
52
53
void test_cholmod_support()
53
void test_cholmod_support()
54
{
54
{
55
  CALL_SUBTEST_1(test_cholmod_T<double>());
55
  CALL_SUBTEST_1(test_cholmod_T<double>());
56
  CALL_SUBTEST_2(test_cholmod_T<std::complex<double> >());
56
  CALL_SUBTEST_2(test_cholmod_T<std::complex<double> >());
57
}
57
}

Return to bug 1095