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 (-1 / +55 lines)
Lines 73-89 cholmod_sparse viewAsCholmod(SparseMatri Link Here
73
73
74
  res.dtype   = 0;
74
  res.dtype   = 0;
75
  res.stype   = -1;
75
  res.stype   = -1;
76
  
76
  
77
  if (internal::is_same<_StorageIndex,int>::value)
77
  if (internal::is_same<_StorageIndex,int>::value)
78
  {
78
  {
79
    res.itype = CHOLMOD_INT;
79
    res.itype = CHOLMOD_INT;
80
  }
80
  }
81
  else if (internal::is_same<_StorageIndex,SuiteSparse_long>::value)
81
  else if (internal::is_same<_StorageIndex,long>::value)
82
  {
82
  {
83
    res.itype = CHOLMOD_LONG;
83
    res.itype = CHOLMOD_LONG;
84
  }
84
  }
85
  else
85
  else
86
  {
86
  {
87
    eigen_assert(false && "Index type not supported yet");
87
    eigen_assert(false && "Index type not supported yet");
88
  }
88
  }
89
89
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
356
	  Map<const Array<Scalar,1,Dynamic>, 0, InnerStride<> >
357
	    s1(x + px[sx], ncols, InnerStride<>(nrows+1));
358
	  logDet += s1.real().log().sum();
359
	}
360
      } else {
361
	// This is a simplicial factorization, which is simply stored as a
362
	// sparse CSC matrix in x, p, i. We want the diagonal, which is
363
	// just the first entry in each column; p gives the offsets in x to
364
	// the beginning of each column.
365
	//
366
	// The ->p array actually has n+1 entries, but only the first n
367
	// entries actually point to real columns (the last entry is a
368
	// sentinel)
369
	StorageIndex *p = (StorageIndex*) cf->p;
370
	for (size_t ex=0; ex < cf->n; ++ex) {
371
	  logDet += log(real( x[p[ex]] ));
372
	}
373
      }
374
      if (cf->is_ll) {
375
	logDet *= 2.0;
376
      }
377
      return logDet;
378
    };
379
326
    template<typename Stream>
380
    template<typename Stream>
327
    void dumpMemory(Stream& /*s*/)
381
    void dumpMemory(Stream& /*s*/)
328
    {}
382
    {}
329
    
383
    
330
  protected:
384
  protected:
331
    mutable cholmod_common m_cholmod;
385
    mutable cholmod_common m_cholmod;
332
    cholmod_factor* m_cholmodFactor;
386
    cholmod_factor* m_cholmodFactor;
333
    RealScalar m_shiftOffset[2];
387
    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