Re: [eigen] LDLt and LLt fixes

[Bill Greene <w.h.greene@xxxxxxxxx>]

Benoit,

   Wow, thanks for the quick work!

   I just downloaded the latest code from svn and tried a quick test. Surprisingly, it didn't work. After some
debugging and investigation, I think this may be an unusual case. My other test cases work correctly--
exactly as I expect. Here is the code for my failed case

void bathe2DOF()
{
  const int n = 2;
  Eigen::Matrix2d K, M;
  K << 6., -2., -2., 4.;
  M << 2., 0., 0., 1.;
  const double omega2 = 3.0;
  Eigen::Matrix2d A = K - omega2*M;
  Eigen::LDLT<Eigen::Matrix2d> factA(A);
  std::cout << "Is invertible: " << factA.isInvertible() << std::endl;
  Eigen::Vector2d b, x;
  b << 0.0, 10.0;
  factA.solve(b, &x);
  prtMat(x, "x");
}

The correct value for x, by the way, is [-5.0, 0]. 

The lapack routines dsytrf and dsytrs compute this solution. I promised to provide more details of the lapack
solution-- and I intend to-- but I'm somewhat puzzled by the factored matrix returned by dsytrf. Even though
the solution from dsytrs is correct, the factored matrix is not what I expected-- or perhaps its just that the
documentation is so cryptic I don't know what I'm seeing. At any rate, I'll provide more lapack details as
soon as I understand them.

Bill

On Mon, Mar 30, 2009 at 5:55 PM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx> wrote:

Hi,

the thread with Bill inspired me... so now (trunk revision 947097),
LDLt supports the negative semidefinite case.

There's a couple of things to discuss.

First, the floating-point comparisons in both LLt and LDLt had issues. LLt did

  if(x < eps)

and LDLt did

  if(Djj <= 0)

and replacing these by proper fuzzy comparisons (experimenting with a
few different ones) resulted in a big improvement in the accuracy of
some methods, like LDLt::rank() which used to be untested and had
typically 3% imprecision, and now is exact in 95% of cases (However
for an exact rank computation, it's still not nearly as reliable as
the full-pivoting LU which can be explained by the different between
bi-directional pivoting and the diagonal pivoting that's all one can
do with Cholesky).

The other issue is that LLt and LDLt have claims in their dox that
they use only a triangular half of the matrix and don't need the other
half to be initialized, well that's not currently the case, one can
check that by uncommenting a line I added to the cholesky.cpp
unit-test. Help fixing this welcome...

Cheers,
Benoit