aboutsummaryrefslogtreecommitdiffhomepage
diff options
context:
space:
mode:
-rw-r--r--Eigen/src/Cholesky/LDLT.h14
-rw-r--r--test/cholesky.cpp32
2 files changed, 36 insertions, 10 deletions
diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h
index f34d26465..b43e85e7f 100644
--- a/Eigen/src/Cholesky/LDLT.h
+++ b/Eigen/src/Cholesky/LDLT.h
@@ -291,13 +291,6 @@ template<> struct ldlt_inplace<Lower>
cutoff = abs(NumTraits<Scalar>::epsilon() * biggest_in_corner);
}
- // Finish early if the matrix is not full rank.
- if(biggest_in_corner < cutoff)
- {
- for(Index i = k; i < size; i++) transpositions.coeffRef(i) = i;
- break;
- }
-
transpositions.coeffRef(k) = index_of_biggest_in_corner;
if(k != index_of_biggest_in_corner)
{
@@ -333,6 +326,7 @@ template<> struct ldlt_inplace<Lower>
if(rs>0)
A21.noalias() -= A20 * temp.head(k);
}
+
if((rs>0) && (abs(mat.coeffRef(k,k)) > cutoff))
A21 /= mat.coeffRef(k,k);
@@ -518,12 +512,12 @@ struct solve_retval<LDLT<_MatrixType,_UpLo>, Rhs>
typedef typename LDLTType::RealScalar RealScalar;
const Diagonal<const MatrixType> vectorD = dec().vectorD();
RealScalar tolerance = (max)(vectorD.array().abs().maxCoeff() * NumTraits<Scalar>::epsilon(),
- RealScalar(1) / NumTraits<RealScalar>::highest()); // motivated by LAPACK's xGELSS
+ RealScalar(1) / NumTraits<RealScalar>::highest()); // motivated by LAPACK's xGELSS
for (Index i = 0; i < vectorD.size(); ++i) {
if(abs(vectorD(i)) > tolerance)
- dst.row(i) /= vectorD(i);
+ dst.row(i) /= vectorD(i);
else
- dst.row(i).setZero();
+ dst.row(i).setZero();
}
// dst = L^-T (D^-1 L^-1 P b)
diff --git a/test/cholesky.cpp b/test/cholesky.cpp
index b980dc572..d4d90e467 100644
--- a/test/cholesky.cpp
+++ b/test/cholesky.cpp
@@ -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