Summary: | log det for cholmod | ||||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
Product: | Eigen | Reporter: | Joshua Pritikin <jpritikin> | ||||||||||||
Component: | Sparse | Assignee: | Nobody <eigen.nobody> | ||||||||||||
Status: | RESOLVED FIXED | ||||||||||||||
Severity: | Feature Request | CC: | chtz, gael.guennebaud | ||||||||||||
Priority: | Normal | ||||||||||||||
Version: | 3.3 (current stable) | ||||||||||||||
Hardware: | All | ||||||||||||||
OS: | All | ||||||||||||||
Whiteboard: | |||||||||||||||
Attachments: |
|
Description
Joshua Pritikin
2015-10-21 08:55:36 UTC
Sure, patch welcome :) Created attachment 636 [details]
proposed change
How does this look?
Created attachment 637 [details]
proposed change
How does this look?
Thanks for the patch, some comments: - The documentation is a bit misleading, and I'd suggest something like: /** \returns the determinant of the underlying matrix from the current factorization */ - I guess that int* should be replaced by StorageIndex* to support 64 bits integers. - complex scalar type can easily be supported by observing that the diagonal entries are real. Use numext::real(val) to extract the real part (works on reals too). - the nested for loop can probably be replaced by something like: logDet += Array<Scalar,Dynamic,Dynamic>::Map(x+px[sx], nrows, ncols).diagonal().real().log().sum(); - Calls like std::bla(x) should be replaced by: using std::bla; ... bla(x) to support custom scalar types. The int pointers point into the cholmod_factor type. This is a C structure provided by cholmod, not Eigen. I suspect it only implements integers. I will address your other suggestions. I think that they can be long int if cholmod_sparse::itype == CHOLMOD_LONG, which is why those pointer members are void* and not int*. Created attachment 639 [details] implement complex and address feedback > the nested for loop can probably be replaced by something like: If you look closely, we don't add up the diagonal of a dense matrix but entries in the first row. Maybe it is possible to use a Map'd Eigen type with TriangularView, but I didn't research it. Created attachment 640 [details]
unremove tests
right, I'have been confused by the +1, then it's even simpler: Map<const Matrix<Scalar,1,Dynamic>, 0, InnerStride<> >(x+px[sx], ncols, InnerStride<>(nrows+1)) Created attachment 644 [details]
add log det
Use fancy Eigen::Map instead of boring loop
thanks for the patch. -- GitLab Migration Automatic Message -- This bug has been migrated to gitlab.com's GitLab instance and has been closed from further activity. You can subscribe and participate further through the new bug through this link to our GitLab instance: https://gitlab.com/libeigen/eigen/issues/1095. |