diff options
| -rw-r--r-- | Eigen/src/Cholesky/LDLT.h | 7 | ||||
| -rw-r--r-- | test/cholesky.cpp | 14 |
2 files changed, 21 insertions, 0 deletions
diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h index 68e54b1d4..a73a9c19f 100644 --- a/Eigen/src/Cholesky/LDLT.h +++ b/Eigen/src/Cholesky/LDLT.h @@ -281,6 +281,13 @@ template<> struct ldlt_inplace<Lower> if(sign) *sign = real(mat.diagonal().coeff(index_of_biggest_in_corner)) > 0 ? 1 : -1; } + else if(sign) + { + // LDLT is not guaranteed to work for indefinite matrices, but let's try to get the sign right + int newSign = real(mat.diagonal().coeff(index_of_biggest_in_corner)) > 0; + if(newSign != *sign) + *sign = 0; + } // Finish early if the matrix is not full rank. if(biggest_in_corner < cutoff) diff --git a/test/cholesky.cpp b/test/cholesky.cpp index 14e01c006..49c79f9c8 100644 --- a/test/cholesky.cpp +++ b/test/cholesky.cpp @@ -262,6 +262,19 @@ template<typename MatrixType> void cholesky_bug241(const MatrixType& m) VERIFY_IS_APPROX(matA * vecX, vecB); } +// LDLT is not guaranteed to work for indefinite matrices, but happens to work fine if matrix is diagonal. +// This test checks that LDLT reports correctly that matrix is indefinite. +// See http://forum.kde.org/viewtopic.php?f=74&t=106942 +template<typename MatrixType> void cholesky_indefinite(const MatrixType& m) +{ + eigen_assert(m.rows() == 2 && m.cols() == 2); + MatrixType mat; + mat << 1, 0, 0, -1; + LDLT<MatrixType> ldlt(mat); + VERIFY(!ldlt.isNegative()); + VERIFY(!ldlt.isPositive()); +} + template<typename MatrixType> void cholesky_verify_assert() { MatrixType tmp; @@ -289,6 +302,7 @@ void test_cholesky() CALL_SUBTEST_1( cholesky(Matrix<double,1,1>()) ); CALL_SUBTEST_3( cholesky(Matrix2d()) ); CALL_SUBTEST_3( cholesky_bug241(Matrix2d()) ); + CALL_SUBTEST_3( cholesky_indefinite(Matrix2d()) ); CALL_SUBTEST_4( cholesky(Matrix3f()) ); CALL_SUBTEST_5( cholesky(Matrix4d()) ); s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE); |