From 919414b9fe2ad7fdcb0f2b2cbdf6b5322d0f2034 Mon Sep 17 00:00:00 2001
From: Christoph Hertzberg <chtz@informatik.uni-bremen.de>
Date: Mon, 3 Dec 2018 16:18:15 +0100
Subject: bug #785: Make Cholesky decomposition work for empty matrices

---
 test/cholesky.cpp | 16 +++++++++++-----
 1 file changed, 11 insertions(+), 5 deletions(-)

(limited to 'test/cholesky.cpp')

diff --git a/test/cholesky.cpp b/test/cholesky.cpp
index b871351e0..e1e8b7bf7 100644
--- a/test/cholesky.cpp
+++ b/test/cholesky.cpp
@@ -19,6 +19,7 @@
 
 template<typename MatrixType, int UpLo>
 typename MatrixType::RealScalar matrix_l1_norm(const MatrixType& m) {
+  if(m.cols()==0) return typename MatrixType::RealScalar(0);
   MatrixType symm = m.template selfadjointView<UpLo>();
   return symm.cwiseAbs().colwise().sum().maxCoeff();
 }
@@ -96,7 +97,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
     RealScalar rcond_est = chollo.rcond();
     // Verify that the estimated condition number is within a factor of 10 of the
     // truth.
-    VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10);
+    VERIFY(rcond_est >= rcond / 10 && rcond_est <= rcond * 10);
 
     // test the upper mode
     LLT<SquareMatrixType,Upper> cholup(symmUp);
@@ -112,12 +113,12 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
     rcond = (RealScalar(1) / matrix_l1_norm<MatrixType, Upper>(symmUp)) /
                              matrix_l1_norm<MatrixType, Upper>(symmUp_inverse);
     rcond_est = cholup.rcond();
-    VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10);
+    VERIFY(rcond_est >= rcond / 10 && rcond_est <= rcond * 10);
 
 
     MatrixType neg = -symmLo;
     chollo.compute(neg);
-    VERIFY(chollo.info()==NumericalIssue);
+    VERIFY(neg.size()==0 || chollo.info()==NumericalIssue);
 
     VERIFY_IS_APPROX(MatrixType(chollo.matrixL().transpose().conjugate()), MatrixType(chollo.matrixU()));
     VERIFY_IS_APPROX(MatrixType(chollo.matrixU().transpose().conjugate()), MatrixType(chollo.matrixL()));
@@ -166,7 +167,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
     RealScalar rcond_est = ldltlo.rcond();
     // Verify that the estimated condition number is within a factor of 10 of the
     // truth.
-    VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10);
+    VERIFY(rcond_est >= rcond / 10 && rcond_est <= rcond * 10);
 
 
     LDLT<SquareMatrixType,Upper> ldltup(symmUp);
@@ -183,7 +184,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
     rcond = (RealScalar(1) / matrix_l1_norm<MatrixType, Upper>(symmUp)) /
                              matrix_l1_norm<MatrixType, Upper>(symmUp_inverse);
     rcond_est = ldltup.rcond();
-    VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10);
+    VERIFY(rcond_est >= rcond / 10 && rcond_est <= rcond * 10);
 
     VERIFY_IS_APPROX(MatrixType(ldltlo.matrixL().transpose().conjugate()), MatrixType(ldltlo.matrixU()));
     VERIFY_IS_APPROX(MatrixType(ldltlo.matrixU().transpose().conjugate()), MatrixType(ldltlo.matrixL()));
@@ -507,6 +508,11 @@ EIGEN_DECLARE_TEST(cholesky)
     CALL_SUBTEST_6( cholesky_cplx(MatrixXcd(s,s)) );
     TEST_SET_BUT_UNUSED_VARIABLE(s)
   }
+  // empty matrix, regression test for Bug 785:
+  CALL_SUBTEST_2( cholesky(MatrixXd(0,0)) );
+
+  // This does not work yet:
+  // CALL_SUBTEST_2( cholesky(Matrix<double,0,0>()) );
 
   CALL_SUBTEST_4( cholesky_verify_assert<Matrix3f>() );
   CALL_SUBTEST_7( cholesky_verify_assert<Matrix3d>() );
-- 
cgit v1.2.3