mirror of
https://gitlab.com/libeigen/eigen.git
synced 2026-04-10 11:34:33 +08:00
Remove early termination in LDLT: the zero on the diagonal of the input matrix does not mean the matrix is not full rank. Typical examples are matrices coming from LS with linear equality constraints.
This commit is contained in:
@@ -179,6 +179,38 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
|
||||
// restore
|
||||
if(sign == -1)
|
||||
symm = -symm;
|
||||
|
||||
// check matrices coming from linear constraints with Lagrange multipliers
|
||||
if(rows>=3)
|
||||
{
|
||||
SquareMatrixType A = symm;
|
||||
int c = internal::random<int>(0,rows-2);
|
||||
A.bottomRightCorner(c,c).setZero();
|
||||
// Make sure a solution exists:
|
||||
vecX.setRandom();
|
||||
vecB = A * vecX;
|
||||
vecX.setZero();
|
||||
ldltlo.compute(A);
|
||||
VERIFY_IS_APPROX(A, ldltlo.reconstructedMatrix());
|
||||
vecX = ldltlo.solve(vecB);
|
||||
VERIFY_IS_APPROX(A * vecX, vecB);
|
||||
}
|
||||
|
||||
// check non-full rank matrices
|
||||
if(rows>=3)
|
||||
{
|
||||
int r = internal::random<int>(1,rows-1);
|
||||
Matrix<Scalar,Dynamic,Dynamic> a = Matrix<Scalar,Dynamic,Dynamic>::Random(rows,r);
|
||||
SquareMatrixType A = a * a.adjoint();
|
||||
// Make sure a solution exists:
|
||||
vecX.setRandom();
|
||||
vecB = A * vecX;
|
||||
vecX.setZero();
|
||||
ldltlo.compute(A);
|
||||
VERIFY_IS_APPROX(A, ldltlo.reconstructedMatrix());
|
||||
vecX = ldltlo.solve(vecB);
|
||||
VERIFY_IS_APPROX(A * vecX, vecB);
|
||||
}
|
||||
}
|
||||
|
||||
// update/downdate
|
||||
|
||||
Reference in New Issue
Block a user