Re: [eigen] Pivoting for LDLT

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

2009/2/3 Keir Mierle <mierle@xxxxxxxxx>:
> I did some more experiments, and I'm convinced. Not only does the example
> you give here fail badly with nopivot, the nopivot version is not
> rank-revealing in the semidefinite case. What's left is the API decision;
> should I make pivoting a template option (enabled by default)? Since I doubt
> many LDLT objects get instantiated this is probably ok.

I've been thinking about this for LU and I think that pivoting or not
is really such a big difference that it's best handled as 2 different
classes.

However in your case I doubt that the non-pivoting option is needed:
if the user knows no pivoting is needed, he might use the LLt, which
is fast and non-pivoting, right? So a non pivoting LDLt would have to
fit in a very small niche market between the LLt and the pivoting
LDLt.

> However, now the matrix L is not meaningful without the permutation matrix.
> This is not an issue for solving, but could be for users of .matrixL(). On
> the other hand, this is currently the case for LU.

It is meaningful, only you have to say that the matrix is no longer
decomposed as LDLt but rather as
(PL)D(PL)t where P is a permutation matrix. (At least that's how i
understand pivoting LDLt, I haven't yet looked at the code).

In the full LU, I ended up removing the matrixL() and matrixU()
methods but for another reason: in the rectangular case, matrixL() was
really really not meaningful at all. Here, I'd say your matrixL() is
very meaningful, it's just that it must not be confused with PL.

That said, i decided in LU to only provide a method matrixLU() and let
the user do his own cooking with part() to get L and U. Maybe this
approach can be considered also for other kinds of decompositions
(again, what made it especially necessary in LU was that a matrixL()
returning a Part was completely meaningless in the rectangular case).

> p.s. the blocked LDLT which can deal with indefinite matrices is a different
> algorithm (D is block-diagonal with 1x1 and 2x2 blocks), and I am not
> targetting it because I only need to solve semidefinite systems stably.

OK, sure.

Cheers,
Benoit