Pulled in the latest changes from the Eigen trunk

1 files changed, 22 insertions, 2 deletions

diff --git a/test/cholesky.cpp b/test/cholesky.cpp
index 569318f83..a883192ab 100644
--- a/test/cholesky.cpp
+++ b/test/cholesky.cpp
@@ -68,6 +68,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
   Index cols = m.cols();
 
   typedef typename MatrixType::Scalar Scalar;
+  typedef typename NumTraits<Scalar>::Real RealScalar;
   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> SquareMatrixType;
   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
 
@@ -180,7 +181,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
     if(rows>=3)
     {
       SquareMatrixType A = symm;
-      int c = internal::random<int>(0,rows-2);
+      Index c = internal::random<Index>(0,rows-2);
       A.bottomRightCorner(c,c).setZero();
       // Make sure a solution exists:
       vecX.setRandom();
@@ -195,7 +196,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
     // check non-full rank matrices
     if(rows>=3)
     {
-      int r = internal::random<int>(1,rows-1);
+      Index r = internal::random<Index>(1,rows-1);
       Matrix<Scalar,Dynamic,Dynamic> a = Matrix<Scalar,Dynamic,Dynamic>::Random(rows,r);
       SquareMatrixType A = a * a.adjoint();
       // Make sure a solution exists:
@@ -207,6 +208,25 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
       vecX = ldltlo.solve(vecB);
       VERIFY_IS_APPROX(A * vecX, vecB);
     }
+    
+    // check matrices with a wide spectrum
+    if(rows>=3)
+    {
+      RealScalar s = (std::min)(16,std::numeric_limits<RealScalar>::max_exponent10/8);
+      Matrix<Scalar,Dynamic,Dynamic> a = Matrix<Scalar,Dynamic,Dynamic>::Random(rows,rows);
+      Matrix<RealScalar,Dynamic,1> d =  Matrix<RealScalar,Dynamic,1>::Random(rows);
+      for(Index k=0; k<rows; ++k)
+        d(k) = d(k)*std::pow(RealScalar(10),internal::random<RealScalar>(-s,s));
+      SquareMatrixType A = a * d.asDiagonal() * a.adjoint();
+      // Make sure a solution exists:
+      vecX.setRandom();
+      vecB = A * vecX;
+      vecX.setZero();
+      ldltlo.compute(A);
+      VERIFY_IS_APPROX(A, ldltlo.reconstructedMatrix());
+      vecX = ldltlo.solve(vecB);
+      VERIFY_IS_APPROX(A * vecX, vecB);
+    }
   }
 
   // update/downdate