* bug fixes in: Dot, generalized eigen problem, singular matrix detetection in Cholesky

* fix all numerical instabilies in the unit tests, now all tests can be run 2000 times
  with almost zero failures.

1 files changed, 19 insertions, 8 deletions

diff --git a/Eigen/src/QR/SelfAdjointEigenSolver.h b/Eigen/src/QR/SelfAdjointEigenSolver.h
index f8bd7bfad..765af7d21 100644
--- a/Eigen/src/QR/SelfAdjointEigenSolver.h
+++ b/Eigen/src/QR/SelfAdjointEigenSolver.h
@@ -225,22 +225,33 @@ void SelfAdjointEigenSolver<MatrixType>::
 compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors)
 {
   ei_assert(matA.cols()==matA.rows() && matB.rows()==matA.rows() && matB.cols()==matB.rows());
-
-  // Compute the cholesky decomposition of matB = U'U
+  
+  // Compute the cholesky decomposition of matB = L L'
   Cholesky<MatrixType> cholB(matB);
 
-  // compute C = inv(U') A inv(U)
-  MatrixType matC = cholB.matrixL().solveTriangular(matA);
-  // FIXME since we currently do not support A * inv(U),
-  // let's do (inv(U') A')' :
-  matC = (cholB.matrixL().solveTriangular(matC.adjoint())).adjoint();
+  // compute C = inv(L) A inv(L')
+  MatrixType matC = matA;
+  cholB.matrixL().solveTriangularInPlace(matC);
+  // FIXME since we currently do not support A * inv(L'), let's do (inv(L) A')' :
+  matC = matC.adjoint().eval();
+  cholB.matrixL().template marked<Lower>().solveTriangularInPlace(matC);
+  matC = matC.adjoint().eval();
+  // this version works too:
+//   matC = matC.transpose();
+//   cholB.matrixL().conjugate().template marked<Lower>().solveTriangularInPlace(matC);
+//   matC = matC.transpose();
+  // FIXME: this should work: (currently it only does for small matrices)
+//   Transpose<MatrixType> trMatC(matC);
+//   cholB.matrixL().conjugate().eval().template marked<Lower>().solveTriangularInPlace(trMatC);
 
   compute(matC, computeEigenvectors);
 
   if (computeEigenvectors)
   {
     // transform back the eigen vectors: evecs = inv(U) * evecs
-    m_eivec = cholB.matrixL().adjoint().template marked<Upper>().solveTriangular(m_eivec);
+    cholB.matrixL().adjoint().template marked<Upper>().solveTriangularInPlace(m_eivec);
+    for (int i=0; i<m_eivec.cols(); ++i)
+      m_eivec.col(i) = m_eivec.col(i).normalized();
   }
 }