Re: [eigen] Instability in LLT and LDLT methods.

[Benoit Jacob <jacob.benoit.1@xxxxxxxxx>]

2009/1/28 Gael Guennebaud <gael.guennebaud@xxxxxxxxx>:
>>>> I think that LAPACK and GSL do a hybrid of partial and full pivoting
>>>> that I have yet to understand.
>>>
>>> That's not what their documentation says:
>>> http://netlib.org/lapack/lug/node38.html
>>> http://www.gnu.org/software/gsl/manual/html_node/LU-Decomposition.html
>>
>> Ah OK. Here I was just repeating stuff I heared over IRC. So I don't
>> know. Perhaps they still do kind of a hybrid but don't mention it in
>> the docs. If they really do plain partial pivoting, then let's try
>> inverting, say, a 2000x2000 matrix, it will fail. For this reason I
>> doubt that they actually do plain partial pivoting.
>
> yeah, GSL clearly only does partial pivoting. About LAPACK I thought I
> read somewhere that the high level block LU algorithm was based on top
> of the sgetc2 auxiliary routines (un-blocked full pivoting), but
> searching over the net again, I found source code using the dgetf2
> routine (un-blocked partial pivoting)... That would make sense because
> I don't know how hybrid full/partial pivoting could work.
>
>
>>> I think a 50% speedup is a worthwhile target. At least for vector solves,
>>> you can compute the residual and redo the computation with complete pivoting
>>> if necessary.
>
> another argument in favor of having partial pivoting is that would
> allow to implement a blocked version with only ~10 lines of code, and
> for large matrix a blocked LU with partial pivoting is much much
> faster than the current one.

OK so I'll really make this change in LU.

I'm really sorry to have generated a lot of noise about partial
pivoting being unusable for large sizes, that was wrong, partial
pivoting is usable as soon as the matrix is invertible.

I'll also correct my false statements on the forum.

Benoit