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 (+58 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
    /** Efficiently obtains the determinant of the factor
327
      *
328
      * \returns the determinant of the factor
329
      */
330
    RealScalar determinant() const {
331
      return std::exp(logDeterminant());
332
    }
333
334
    /** Efficiently obtains the log determinant of the factor
335
      *
336
      * \returns the log determinant of the factor
337
      */
338
    RealScalar logDeterminant() const {
339
      eigen_assert((internal::is_same<Scalar,float>::value || internal::is_same<Scalar,double>::value) &&
340
		   "Determinant is not implemented for complex matrices");
341
      eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
342
343
      // Based on https://github.com/njsmith/scikits-sparse/blob/master/scikits/sparse/cholmod.pyx
344
      cholmod_factor *cf = m_cholmodFactor;
345
      RealScalar logDet = 0;
346
      Scalar *x = (Scalar*) cf->x;
347
      if (cf->is_super) {
348
	// This is a supernodal factorization, which is stored as a bunch
349
	// of dense, lower-triangular, column-major arrays packed into the
350
	// x vector. This is not documented in the CHOLMOD user-guide, or
351
	// anywhere else as far as I can tell; I got the details from
352
	// CVXOPT's C/cholmod.c.
353
354
	int *super = (int*) cf->super;
355
	int *pi = (int*) cf->pi;
356
	int *px = (int*) cf->px;
357
	for (size_t sx=0; sx < cf->nsuper; ++sx) {
358
	  int ncols = super[sx + 1] - super[sx];
359
	  int nrows = pi[sx + 1] - pi[sx];
360
	  for (int cx=px[sx]; cx < px[sx] + nrows * ncols; cx += nrows+1) {
361
	    logDet += std::log(x[cx]);
362
	  }
363
	}
364
      } else {
365
	// This is a simplicial factorization, which is simply stored as a
366
	// sparse CSC matrix in x, p, i. We want the diagonal, which is
367
	// just the first entry in each column; p gives the offsets in x to
368
	// the beginning of each column.
369
	//
370
	// The ->p array actually has n+1 entries, but only the first n
371
	// entries actually point to real columns (the last entry is a
372
	// sentinel)
373
	int *p = (int*) cf->p;
374
	for (size_t ex=0; ex < cf->n; ++ex) {
375
	  logDet += std::log( x[p[ex]] );
376
	}
377
      }
378
      if (cf->is_ll) {
379
	logDet *= 2.0;
380
      }
381
      return logDet;
382
    };
383
326
    template<typename Stream>
384
    template<typename Stream>
327
    void dumpMemory(Stream& /*s*/)
385
    void dumpMemory(Stream& /*s*/)
328
    {}
386
    {}
329
    
387
    
330
  protected:
388
  protected:
331
    mutable cholmod_common m_cholmod;
389
    mutable cholmod_common m_cholmod;
332
    cholmod_factor* m_cholmodFactor;
390
    cholmod_factor* m_cholmodFactor;
333
    RealScalar m_shiftOffset[2];
391
    RealScalar m_shiftOffset[2];
(-)a/test/cholmod_support.cpp (-10 / +22 lines)
Lines 7-23 Link Here
7
// Public License v. 2.0. If a copy of the MPL was not distributed
7
// Public License v. 2.0. If a copy of the MPL was not distributed
8
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
8
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9
9
10
#define EIGEN_NO_DEBUG_SMALL_PRODUCT_BLOCKS
10
#define EIGEN_NO_DEBUG_SMALL_PRODUCT_BLOCKS
11
#include "sparse_solver.h"
11
#include "sparse_solver.h"
12
12
13
#include <Eigen/CholmodSupport>
13
#include <Eigen/CholmodSupport>
14
14
15
template<typename T> void test_cholmod_T()
15
template<typename T> void test1_cholmod_T()
16
{
16
{
17
  CholmodDecomposition<SparseMatrix<T>, Lower> g_chol_colmajor_lower; g_chol_colmajor_lower.setMode(CholmodSupernodalLLt);
17
  CholmodDecomposition<SparseMatrix<T>, Lower> g_chol_colmajor_lower; g_chol_colmajor_lower.setMode(CholmodSupernodalLLt);
18
  CholmodDecomposition<SparseMatrix<T>, Upper> g_chol_colmajor_upper; g_chol_colmajor_upper.setMode(CholmodSupernodalLLt);
18
  CholmodDecomposition<SparseMatrix<T>, Upper> g_chol_colmajor_upper; g_chol_colmajor_upper.setMode(CholmodSupernodalLLt);
19
  CholmodDecomposition<SparseMatrix<T>, Lower> g_llt_colmajor_lower;  g_llt_colmajor_lower.setMode(CholmodSimplicialLLt);
19
  CholmodDecomposition<SparseMatrix<T>, Lower> g_llt_colmajor_lower;  g_llt_colmajor_lower.setMode(CholmodSimplicialLLt);
20
  CholmodDecomposition<SparseMatrix<T>, Upper> g_llt_colmajor_upper;  g_llt_colmajor_upper.setMode(CholmodSimplicialLLt);
20
  CholmodDecomposition<SparseMatrix<T>, Upper> g_llt_colmajor_upper;  g_llt_colmajor_upper.setMode(CholmodSimplicialLLt);
21
  CholmodDecomposition<SparseMatrix<T>, Lower> g_ldlt_colmajor_lower; g_ldlt_colmajor_lower.setMode(CholmodLDLt);
21
  CholmodDecomposition<SparseMatrix<T>, Lower> g_ldlt_colmajor_lower; g_ldlt_colmajor_lower.setMode(CholmodLDLt);
22
  CholmodDecomposition<SparseMatrix<T>, Upper> g_ldlt_colmajor_upper; g_ldlt_colmajor_upper.setMode(CholmodLDLt);
22
  CholmodDecomposition<SparseMatrix<T>, Upper> g_ldlt_colmajor_upper; g_ldlt_colmajor_upper.setMode(CholmodLDLt);
23
  
23
  
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
46
//  check_sparse_spd_determinant(chol_colmajor_upper);
46
template<typename T> void test2_cholmod_T()
47
//  check_sparse_spd_determinant(llt_colmajor_lower);
47
{
48
//  check_sparse_spd_determinant(llt_colmajor_upper);
48
  test1_cholmod_T<T>();
49
//  check_sparse_spd_determinant(ldlt_colmajor_lower);
49
50
//  check_sparse_spd_determinant(ldlt_colmajor_upper);
50
  CholmodSupernodalLLT<SparseMatrix<T>, Lower> chol_colmajor_lower;
51
  CholmodSupernodalLLT<SparseMatrix<T>, Upper> chol_colmajor_upper;
52
  CholmodSimplicialLLT<SparseMatrix<T>, Lower> llt_colmajor_lower;
53
  CholmodSimplicialLLT<SparseMatrix<T>, Upper> llt_colmajor_upper;
54
  CholmodSimplicialLDLT<SparseMatrix<T>, Lower> ldlt_colmajor_lower;
55
  CholmodSimplicialLDLT<SparseMatrix<T>, Upper> ldlt_colmajor_upper;
56
57
  check_sparse_spd_determinant(chol_colmajor_lower);
58
  check_sparse_spd_determinant(chol_colmajor_upper);
59
  check_sparse_spd_determinant(llt_colmajor_lower);
60
  check_sparse_spd_determinant(llt_colmajor_upper);
61
  check_sparse_spd_determinant(ldlt_colmajor_lower);
62
  check_sparse_spd_determinant(ldlt_colmajor_upper);
51
}
63
}
52
64
53
void test_cholmod_support()
65
void test_cholmod_support()
54
{
66
{
55
  CALL_SUBTEST_1(test_cholmod_T<double>());
67
  CALL_SUBTEST_1(test2_cholmod_T<double>());
56
  CALL_SUBTEST_2(test_cholmod_T<std::complex<double> >());
68
  CALL_SUBTEST_2(test1_cholmod_T<std::complex<double> >());
57
}
69
}

Return to bug 1095