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]; |