From 919414b9fe2ad7fdcb0f2b2cbdf6b5322d0f2034 Mon Sep 17 00:00:00 2001 From: Christoph Hertzberg Date: Mon, 3 Dec 2018 16:18:15 +0100 Subject: [PATCH] bug #785: Make Cholesky decomposition work for empty matrices --- Eigen/src/Cholesky/LDLT.h | 3 ++- Eigen/src/Core/ConditionEstimator.h | 2 +- test/cholesky.cpp | 16 +++++++++++----- 3 files changed, 14 insertions(+), 7 deletions(-) diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h index 2dfeac333..6831eab3d 100644 --- a/Eigen/src/Cholesky/LDLT.h +++ b/Eigen/src/Cholesky/LDLT.h @@ -304,7 +304,8 @@ template<> struct ldlt_inplace if (size <= 1) { transpositions.setIdentity(); - if (numext::real(mat.coeff(0,0)) > static_cast(0) ) sign = PositiveSemiDef; + if(size==0) sign = ZeroSign; + else if (numext::real(mat.coeff(0,0)) > static_cast(0) ) sign = PositiveSemiDef; else if (numext::real(mat.coeff(0,0)) < static_cast(0)) sign = NegativeSemiDef; else sign = ZeroSign; return true; diff --git a/Eigen/src/Core/ConditionEstimator.h b/Eigen/src/Core/ConditionEstimator.h index aa7efdc76..8c1a89129 100644 --- a/Eigen/src/Core/ConditionEstimator.h +++ b/Eigen/src/Core/ConditionEstimator.h @@ -160,7 +160,7 @@ rcond_estimate_helper(typename Decomposition::RealScalar matrix_norm, const Deco { typedef typename Decomposition::RealScalar RealScalar; eigen_assert(dec.rows() == dec.cols()); - if (dec.rows() == 0) return RealScalar(1); + if (dec.rows() == 0) return RealScalar(1)/RealScalar(0); if (matrix_norm == RealScalar(0)) return RealScalar(0); if (dec.rows() == 1) return RealScalar(1); const RealScalar inverse_matrix_norm = rcond_invmatrix_L1_norm_estimate(dec); 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::RealScalar matrix_l1_norm(const MatrixType& m) { + if(m.cols()==0) return typename MatrixType::RealScalar(0); MatrixType symm = m.template selfadjointView(); return symm.cwiseAbs().colwise().sum().maxCoeff(); } @@ -96,7 +97,7 @@ template 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 cholup(symmUp); @@ -112,12 +113,12 @@ template void cholesky(const MatrixType& m) rcond = (RealScalar(1) / matrix_l1_norm(symmUp)) / matrix_l1_norm(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 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 ldltup(symmUp); @@ -183,7 +184,7 @@ template void cholesky(const MatrixType& m) rcond = (RealScalar(1) / matrix_l1_norm(symmUp)) / matrix_l1_norm(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()) ); CALL_SUBTEST_4( cholesky_verify_assert() ); CALL_SUBTEST_7( cholesky_verify_assert() );