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