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@xxxxxxxxxx> 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.
>>>
>>>
>>
>
>