really fix the LDLt, at the expense of letting isPositiveDefinite() always return true (it was fundamentally broken anyway, especially as in 2.0 we don't even pivot at all).

also, fix compilation
This commit is contained in:
Benoit Jacob 2010-02-23 18:10:31 -05:00
parent 6b3f81b414
commit 5f42104e0a
2 changed files with 6 additions and 10 deletions

View File

@ -96,8 +96,7 @@ void LDLT<MatrixType>::compute(const MatrixType& a)
assert(a.rows()==a.cols()); assert(a.rows()==a.cols());
const int size = a.rows(); const int size = a.rows();
m_matrix.resize(size, size); m_matrix.resize(size, size);
m_isPositiveDefinite = true; m_isPositiveDefinite = true; // always true. This decomposition is not rank-revealing anyway.
const RealScalar eps = precision<Scalar>() * a.cwise().abs().maxCoeff();
if (size<=1) if (size<=1)
{ {
@ -121,12 +120,6 @@ void LDLT<MatrixType>::compute(const MatrixType& a)
RealScalar tmp = ei_real(a.coeff(j,j) - (m_matrix.row(j).start(j) * m_matrix.col(j).start(j).conjugate()).coeff(0,0)); RealScalar tmp = ei_real(a.coeff(j,j) - (m_matrix.row(j).start(j) * m_matrix.col(j).start(j).conjugate()).coeff(0,0));
m_matrix.coeffRef(j,j) = tmp; m_matrix.coeffRef(j,j) = tmp;
if (tmp < eps)
{
m_isPositiveDefinite = false;
return;
}
int endSize = size-j-1; int endSize = size-j-1;
if (endSize>0) if (endSize>0)
{ {
@ -136,7 +129,8 @@ void LDLT<MatrixType>::compute(const MatrixType& a)
m_matrix.row(j).end(endSize) = a.row(j).end(endSize).conjugate() m_matrix.row(j).end(endSize) = a.row(j).end(endSize).conjugate()
- _temporary.end(endSize).transpose(); - _temporary.end(endSize).transpose();
m_matrix.col(j).end(endSize) = m_matrix.row(j).end(endSize) / tmp; if(tmp != RealScalar(0))
m_matrix.col(j).end(endSize) = m_matrix.row(j).end(endSize) / tmp;
} }
} }
} }
@ -192,7 +186,7 @@ template<typename Derived>
inline const LDLT<typename MatrixBase<Derived>::PlainMatrixType> inline const LDLT<typename MatrixBase<Derived>::PlainMatrixType>
MatrixBase<Derived>::ldlt() const MatrixBase<Derived>::ldlt() const
{ {
return LDLT<PlainObject>(derived()); return LDLT<PlainMatrixType>(derived());
} }
#endif // EIGEN_LDLT_H #endif // EIGEN_LDLT_H

View File

@ -100,6 +100,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
VERIFY_IS_APPROX(symm * matX, matB); VERIFY_IS_APPROX(symm * matX, matB);
} }
#if 0 // cholesky is not rank-revealing anyway
// test isPositiveDefinite on non definite matrix // test isPositiveDefinite on non definite matrix
if (rows>4) if (rows>4)
{ {
@ -109,6 +110,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
LDLT<SquareMatrixType> cholnosqrt(symm); LDLT<SquareMatrixType> cholnosqrt(symm);
VERIFY(!cholnosqrt.isPositiveDefinite()); VERIFY(!cholnosqrt.isPositiveDefinite());
} }
#endif
} }
void test_cholesky() void test_cholesky()