aboutsummaryrefslogtreecommitdiffhomepage
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Eigen/src/Cholesky/LDLT.h52
-rw-r--r--test/cholesky.cpp24
2 files changed, 52 insertions, 24 deletions
diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h
index 4faedd257..f34d26465 100644
--- a/Eigen/src/Cholesky/LDLT.h
+++ b/Eigen/src/Cholesky/LDLT.h
@@ -16,7 +16,10 @@
namespace Eigen {
namespace internal {
-template<typename MatrixType, int UpLo> struct LDLT_Traits;
+ template<typename MatrixType, int UpLo> struct LDLT_Traits;
+
+ // PositiveSemiDef means positive semi-definite and non-zero; same for NegativeSemiDef
+ enum SignMatrix { PositiveSemiDef, NegativeSemiDef, ZeroSign, Indefinite };
}
/** \ingroup Cholesky_Module
@@ -69,7 +72,12 @@ template<typename _MatrixType, int _UpLo> class LDLT
* The default constructor is useful in cases in which the user intends to
* perform decompositions via LDLT::compute(const MatrixType&).
*/
- LDLT() : m_matrix(), m_transpositions(), m_isInitialized(false) {}
+ LDLT()
+ : m_matrix(),
+ m_transpositions(),
+ m_sign(internal::ZeroSign),
+ m_isInitialized(false)
+ {}
/** \brief Default Constructor with memory preallocation
*
@@ -81,6 +89,7 @@ template<typename _MatrixType, int _UpLo> class LDLT
: m_matrix(size, size),
m_transpositions(size),
m_temporary(size),
+ m_sign(internal::ZeroSign),
m_isInitialized(false)
{}
@@ -93,6 +102,7 @@ template<typename _MatrixType, int _UpLo> class LDLT
: m_matrix(matrix.rows(), matrix.cols()),
m_transpositions(matrix.rows()),
m_temporary(matrix.rows()),
+ m_sign(internal::ZeroSign),
m_isInitialized(false)
{
compute(matrix);
@@ -139,7 +149,7 @@ template<typename _MatrixType, int _UpLo> class LDLT
inline bool isPositive() const
{
eigen_assert(m_isInitialized && "LDLT is not initialized.");
- return m_sign == 1;
+ return m_sign == internal::PositiveSemiDef || m_sign == internal::ZeroSign;
}
#ifdef EIGEN2_SUPPORT
@@ -153,7 +163,7 @@ template<typename _MatrixType, int _UpLo> class LDLT
inline bool isNegative(void) const
{
eigen_assert(m_isInitialized && "LDLT is not initialized.");
- return m_sign == -1;
+ return m_sign == internal::NegativeSemiDef || m_sign == internal::ZeroSign;
}
/** \returns a solution x of \f$ A x = b \f$ using the current decomposition of A.
@@ -235,7 +245,7 @@ template<typename _MatrixType, int _UpLo> class LDLT
MatrixType m_matrix;
TranspositionType m_transpositions;
TmpMatrixType m_temporary;
- int m_sign;
+ internal::SignMatrix m_sign;
bool m_isInitialized;
};
@@ -246,7 +256,7 @@ template<int UpLo> struct ldlt_inplace;
template<> struct ldlt_inplace<Lower>
{
template<typename MatrixType, typename TranspositionType, typename Workspace>
- static bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, int* sign=0)
+ static bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, SignMatrix& sign)
{
using std::abs;
typedef typename MatrixType::Scalar Scalar;
@@ -258,8 +268,9 @@ template<> struct ldlt_inplace<Lower>
if (size <= 1)
{
transpositions.setIdentity();
- if(sign)
- *sign = numext::real(mat.coeff(0,0))>0 ? 1:-1;
+ if (numext::real(mat.coeff(0,0)) > 0) sign = PositiveSemiDef;
+ else if (numext::real(mat.coeff(0,0)) < 0) sign = NegativeSemiDef;
+ else sign = ZeroSign;
return true;
}
@@ -284,7 +295,6 @@ template<> struct ldlt_inplace<Lower>
if(biggest_in_corner < cutoff)
{
for(Index i = k; i < size; i++) transpositions.coeffRef(i) = i;
- if(sign) *sign = 0;
break;
}
@@ -325,15 +335,15 @@ template<> struct ldlt_inplace<Lower>
}
if((rs>0) && (abs(mat.coeffRef(k,k)) > cutoff))
A21 /= mat.coeffRef(k,k);
-
- if(sign)
- {
- // LDLT is not guaranteed to work for indefinite matrices, but let's try to get the sign right
- int newSign = numext::real(mat.diagonal().coeff(index_of_biggest_in_corner)) > 0;
- if(k == 0)
- *sign = newSign;
- else if(*sign != newSign)
- *sign = 0;
+
+ RealScalar realAkk = numext::real(mat.coeffRef(k,k));
+ if (sign == PositiveSemiDef) {
+ if (realAkk < 0) sign = Indefinite;
+ } else if (sign == NegativeSemiDef) {
+ if (realAkk > 0) sign = Indefinite;
+ } else if (sign == ZeroSign) {
+ if (realAkk > 0) sign = PositiveSemiDef;
+ else if (realAkk < 0) sign = NegativeSemiDef;
}
}
@@ -399,7 +409,7 @@ template<> struct ldlt_inplace<Lower>
template<> struct ldlt_inplace<Upper>
{
template<typename MatrixType, typename TranspositionType, typename Workspace>
- static EIGEN_STRONG_INLINE bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, int* sign=0)
+ static EIGEN_STRONG_INLINE bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, SignMatrix& sign)
{
Transpose<MatrixType> matt(mat);
return ldlt_inplace<Lower>::unblocked(matt, transpositions, temp, sign);
@@ -445,7 +455,7 @@ LDLT<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::compute(const MatrixType& a)
m_isInitialized = false;
m_temporary.resize(size);
- internal::ldlt_inplace<UpLo>::unblocked(m_matrix, m_transpositions, m_temporary, &m_sign);
+ internal::ldlt_inplace<UpLo>::unblocked(m_matrix, m_transpositions, m_temporary, m_sign);
m_isInitialized = true;
return *this;
@@ -473,7 +483,7 @@ LDLT<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::rankUpdate(const MatrixBase<Deri
for (Index i = 0; i < size; i++)
m_transpositions.coeffRef(i) = i;
m_temporary.resize(size);
- m_sign = sigma>=0 ? 1 : -1;
+ m_sign = sigma>=0 ? internal::PositiveSemiDef : internal::NegativeSemiDef;
m_isInitialized = true;
}
diff --git a/test/cholesky.cpp b/test/cholesky.cpp
index 378525a83..b980dc572 100644
--- a/test/cholesky.cpp
+++ b/test/cholesky.cpp
@@ -263,8 +263,8 @@ template<typename MatrixType> void cholesky_bug241(const MatrixType& m)
// 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)
+// See http://forum.kde.org/viewtopic.php?f=74&t=106942 and bug 736
+template<typename MatrixType> void cholesky_definiteness(const MatrixType& m)
{
eigen_assert(m.rows() == 2 && m.cols() == 2);
MatrixType mat;
@@ -280,6 +280,24 @@ template<typename MatrixType> void cholesky_indefinite(const MatrixType& m)
VERIFY(!ldlt.isNegative());
VERIFY(!ldlt.isPositive());
}
+ {
+ mat << 0, 0, 0, 0;
+ LDLT<MatrixType> ldlt(mat);
+ VERIFY(ldlt.isNegative());
+ VERIFY(ldlt.isPositive());
+ }
+ {
+ mat << 0, 0, 0, 1;
+ LDLT<MatrixType> ldlt(mat);
+ VERIFY(!ldlt.isNegative());
+ VERIFY(ldlt.isPositive());
+ }
+ {
+ mat << -1, 0, 0, 0;
+ LDLT<MatrixType> ldlt(mat);
+ VERIFY(ldlt.isNegative());
+ VERIFY(!ldlt.isPositive());
+ }
}
template<typename MatrixType> void cholesky_verify_assert()
@@ -309,7 +327,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_3( cholesky_definiteness(Matrix2d()) );
CALL_SUBTEST_4( cholesky(Matrix3f()) );
CALL_SUBTEST_5( cholesky(Matrix4d()) );
s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE);
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-04 14:24:00 +0000