[eigen] Instability in LLT and LDLT methods.

[Keir Mierle <mierle@xxxxxxxxx>]

Here is a testcase that fails with LLT and LDLT but works fine with
all of LU, SVD, and QR solving. Depends on my previous patch for QR
solver (or comment out the qr().solve line).

Keir

#include <Eigen/Core>
#include <Eigen/Cholesky>
#include <Eigen/LU>
#include <Eigen/QR>
#include <Eigen/SVD>

using namespace Eigen;
using namespace std;

template<typename Decomposition, typename TMat, typename TVec>
void Solve(const Decomposition &decomposition, const TMat A, const TVec &b) {
  TVec x;
  if (decomposition.solve(b, &x)) {
    cout << "Solved.\n";
    cout << "x = " << x.transpose() << endl;
    cout << "||Ax - b|| = " << (A*x - b).norm() << endl;
  } else {
    cout << "Failed to solve.\n";
  }
  cout << "\n";
}

int main() {
  Matrix<double, 3, 3> A;
  Matrix<double, 3, 1> b;

  A << 7.685e-06,  0,          -2.316e-09,
       0,          7.686e-06,  -3.276e-09,
       -2.316e-09, -3.276e-09,  2.892e-06;
  b << 1.855e-09, 2.623e-09, -2.689e-12;

  cout << "A = \n" << A << endl;
  cout << "b = " << b.transpose() << endl;
  cout << "\n";

  cout << "ldlt: "; Solve(A.ldlt(), A, b);
  cout << "llt:  "; Solve(A.llt(),  A, b);
  cout << "lu:   "; Solve(A.lu(),   A, b);
  cout << "svd:  "; Solve(A.svd(),  A, b);
  cout << "qr:   "; Solve(A.qr(),   A, b);
}