Remove early termination in LDLT: the zero on the diagonal of the input matrix does not mean the matrix is not full rank. Typical examples are matrices coming from LS with linear equality constraints.

1 files changed, 32 insertions, 0 deletions

diff --git a/test/cholesky.cpp b/test/cholesky.cpp
index b980dc572..d4d90e467 100644
--- a/test/cholesky.cpp
+++ b/test/cholesky.cpp
@@ -179,6 +179,38 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
     // restore
     if(sign == -1)
       symm = -symm;
+    
+    // check matrices coming from linear constraints with Lagrange multipliers
+    if(rows>=3)
+    {
+      SquareMatrixType A = symm;
+      int c = internal::random<int>(0,rows-2);
+      A.bottomRightCorner(c,c).setZero();
+      // 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);
+    }
+    
+    // check non-full rank matrices
+    if(rows>=3)
+    {
+      int r = internal::random<int>(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:
+      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