Attachment #644 for bug #1095

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





  res.dtype   = 0;
  res.stype   = -1;
  
  if (internal::is_same<_StorageIndex,int>::value)
  {
    res.itype = CHOLMOD_INT;
  }
  else if (internal::is_same<_StorageIndex,long>::value)
  {
    res.itype = CHOLMOD_LONG;
  }
  else
  {
    eigen_assert(false && "Index type not supported yet");
  }


      * \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 */
    Scalar determinant() const {
      using std::exp;
      return exp(logDeterminant());
    }

    /** \returns the log determinant of the underlying matrix from the current factorization */
    Scalar 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];

	  Map<const Array<Scalar,1,Dynamic>, 0, InnerStride<> >
	    s1(x + px[sx], ncols, InnerStride<>(nrows+1));
	  logDet += s1.real().log().sum();
	}
      } 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];




  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> >());
}

Return to bug 1095

Lines 73-89 cholmod_sparse viewAsCholmod(SparseMatri Link Here

(-)a/Eigen/src/CholmodSupport/CholmodSupport.h (-1 / +55 lines)
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];

Lines 36-57 template<typename T> void test_cholmod_T Link Here

(-)a/test/cholmod_support.cpp (-7 / +7 lines)
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	}