GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-08-11 11:19:02 +08:00
Code Issues Packages Projects Releases Wiki Activity

LDLT: Report sign consistent with D for indefinite matrices.

Browse Source 
See http://forum.kde.org/viewtopic.php?f=74&t=106942
This commit is contained in:
Jitse Niesen 2012-07-22 21:39:38 +01:00
parent 907f4562ac
commit fd5749f51c
2 changed files with 21 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
Eigen/src/Cholesky/LDLT.h
View File

@ -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)

14
test/cholesky.cpp
View File

@ -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);