* 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.

5 files changed, 27 insertions, 14 deletions

diff --git a/Eigen/src/Cholesky/Cholesky.h b/Eigen/src/Cholesky/Cholesky.h
index af5dfb430..891a86a79 100644
--- a/Eigen/src/Cholesky/Cholesky.h
+++ b/Eigen/src/Cholesky/Cholesky.h
@@ -93,17 +93,18 @@ void Cholesky<MatrixType>::compute(const MatrixType& a)
   assert(a.rows()==a.cols());
   const int size = a.rows();
   m_matrix.resize(size, size);
+  const RealScalar eps = ei_sqrt(precision<Scalar>());
 
   RealScalar x;
   x = ei_real(a.coeff(0,0));
-  m_isPositiveDefinite = x > precision<Scalar>() && ei_isMuchSmallerThan(ei_imag(a.coeff(0,0)), RealScalar(1));
+  m_isPositiveDefinite = x > eps && ei_isMuchSmallerThan(ei_imag(a.coeff(0,0)), RealScalar(1));
   m_matrix.coeffRef(0,0) = ei_sqrt(x);
   m_matrix.col(0).end(size-1) = a.row(0).end(size-1).adjoint() / ei_real(m_matrix.coeff(0,0));
   for (int j = 1; j < size; ++j)
   {
     Scalar tmp = ei_real(a.coeff(j,j)) - m_matrix.row(j).start(j).norm2();
     x = ei_real(tmp);
-    if (x < precision<Scalar>() || (!ei_isMuchSmallerThan(ei_imag(tmp), RealScalar(1))))
+    if (x < eps || (!ei_isMuchSmallerThan(ei_imag(tmp), RealScalar(1))))
     {
       m_isPositiveDefinite = false;
       return;
diff --git a/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h b/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h
index b00dc0a11..db33b04f9 100644
--- a/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h
+++ b/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h
@@ -94,6 +94,7 @@ void CholeskyWithoutSquareRoot<MatrixType>::compute(const MatrixType& a)
   const int size = a.rows();
   m_matrix.resize(size, size);
   m_isPositiveDefinite = true;
+  const RealScalar eps = ei_sqrt(precision<Scalar>());
 
   // Let's preallocate a temporay vector to evaluate the matrix-vector product into it.
   // Unlike the standard Cholesky decomposition, here we cannot evaluate it to the destination
@@ -111,7 +112,7 @@ void CholeskyWithoutSquareRoot<MatrixType>::compute(const MatrixType& a)
     RealScalar tmp = ei_real(a.coeff(j,j) - (m_matrix.row(j).start(j) * m_matrix.col(j).start(j).conjugate()).coeff(0,0));
     m_matrix.coeffRef(j,j) = tmp;
 
-    if (ei_isMuchSmallerThan(tmp,RealScalar(1)))
+    if (tmp < eps)
     {
       m_isPositiveDefinite = false;
       return;
diff --git a/Eigen/src/Core/Dot.h b/Eigen/src/Core/Dot.h
index eb25185b6..a3e1229f8 100644
--- a/Eigen/src/Core/Dot.h
+++ b/Eigen/src/Core/Dot.h
@@ -229,9 +229,9 @@ struct ei_dot_impl<Derived1, Derived2, LinearVectorization, CompleteUnrolling>
   };
   static Scalar run(const Derived1& v1, const Derived2& v2)
   {
-    Scalar res =  ei_predux(ei_dot_vec_unroller<Derived1, Derived2, 0, VectorizationSize>::run(v1, v2));
+    Scalar res = ei_predux(ei_dot_vec_unroller<Derived1, Derived2, 0, VectorizationSize>::run(v1, v2));
     if (VectorizationSize != Size)
-      res += ei_dot_novec_unroller<Derived1, Derived2, VectorizationSize, Size>::run(v1, v2);
+      res += ei_dot_novec_unroller<Derived1, Derived2, VectorizationSize, Size-VectorizationSize>::run(v1, v2);
     return res;
   }
 };
diff --git a/Eigen/src/Geometry/AngleAxis.h b/Eigen/src/Geometry/AngleAxis.h
index 733f273d7..cd18bfdec 100644
--- a/Eigen/src/Geometry/AngleAxis.h
+++ b/Eigen/src/Geometry/AngleAxis.h
@@ -131,7 +131,7 @@ template<typename Scalar>
 AngleAxis<Scalar>& AngleAxis<Scalar>::operator=(const QuaternionType& q)
 {
   Scalar n2 = q.vec().norm2();
-  if (ei_isMuchSmallerThan(n2,Scalar(1)))
+  if (n2 < precision<Scalar>()*precision<Scalar>())
   {
     m_angle = 0;
     m_axis << 1, 0, 0;
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();
   }
 }