Re: [eigen] SimplicialCholesky - subclass or patch?

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

The multiple class solution is more consistent with the rest of Eigen
that a template parameter.

gael

On Thu, Jun 30, 2011 at 2:24 PM, Douglas Bates <bates@xxxxxxxxxxxxx> wrote:
> An alternative I was considering is to incorporate LDLT as a template
> parameter. Having an instance variable that can be reset seems unnecessary.
>
> On Jun 30, 2011 1:51 AM, "Gael Guennebaud" <gael.guennebaud@xxxxxxxxx>
> wrote:
>> hm, good point. That would suggest to have two classes, SimplicialLDLt
>> and SimplicialLLt that could be based on a common
>> SimplicialCholeskyBase class.
>>
>> gael
>>
>> On Wed, Jun 29, 2011 at 12:19 AM, Douglas Bates <bates@xxxxxxxxxxxxx>
>> wrote:
>>> On Tue, Jun 28, 2011 at 2:59 PM, Tom Makin <t.makin@xxxxxxxxxx> wrote:
>>>> I would very much support any additions to Eigen that allow direct
>>>> access to
>>>> the L matrix in a sparse Cholesky decomposition. I am currently working
>>>> on a
>>>> surrogate modelling application and have recently switched over to Eigen
>>>> in
>>>> all my code. The main reason I need access to L is to update the
>>>> factorisation without having to re-factorise the whole thing when I add
>>>> a
>>>> couple of extra rows. I am comfortable doing this for the dense case but
>>>> I
>>>> am yet to get to grips with the sparse version because I can't get at L
>>>> directly. I have been able to do a limited amount of sparse updating
>>>> using
>>>> cholmod, but I would like to try and implement it in way that is
>>>> independent
>>>> of the sparse solver.
>>>>
>>>> Douglas: Would you be willing to send me your SimplicalCholesky
>>>> subclass?
>>>
>>> I spoke too soon.  I had not yet tested the code that I wrote and,
>>> naturally, it blew up on me.
>>>
>>> The problem I am encountering is that the return type of a matrixL
>>> method will be determined by the MatrixType and the current value of
>>> m_LDLt. If m_LDLt is true the return type should be
>>> TriangularView<MatrixType, UnitLower>, otherwise
>>> TriangularView<MatrixType, Lower>.  Because m_LDLt can be set or unset
>>> for an instance of the class, I don't think there is a way to
>>> determine at compile time which it should be.  The only way I can
>>> think of to give a definite return type is always to return the Lower
>>> form and, when LDLt is true, add the identity to the strict lower
>>> triangle, which wastes some storage, and will have less efficient
>>> solve methods than are otherwise possible.
>>>
>>> Am I missing something here?
>>>
>>>
>>>> On 28/06/2011 19:00, Douglas Bates wrote:
>>>>>
>>>>> I realize that the Sparse module is not stable in Eigen 3.0.1 and the
>>>>> SimplicialCholesky class is in the unsupported directory.
>>>>> Nevertheless, these are exactly what I need for a central calculation
>>>>> in my research.  (If you want to know the details I include a draft of
>>>>> a paper describing the particular calculations being done.  You can
>>>>> skip to sections 2.3, 2.4 and 2.6 for the relevant steps).
>>>>>
>>>>> Current versions of our code are based on SparseSuite with a rather
>>>>> clunky interface from R.  I would much prefer to use Eigen but, as
>>>>> shown in the enclosed description, I need to access the triangular
>>>>> matrix L from the LLt decomposition and I need a method like
>>>>> logAbsDeterminant for the SimplicialCholesky classes.  I can do this
>>>>> by subclassing the SimplicialCholesky templated class and have done
>>>>> so.  A better long term solution would be to produce and submit
>>>>> patches for SimplicialCholesky.h but before doing so I would like to
>>>>> check that it would be worthwhile.  If that code is not considered
>>>>> viable in the long term then I will stick with subclassing rather than
>>>>> modifying that class.  Could someone from the core group indicate if
>>>>> you would want to receive patches for that file?  I can't guarantee
>>>>> that I will get all the myriad details of the Traits right the first
>>>>> time around but I am willing to learn.
>>>>
>>>>
>>>
>>
>>
>