diff options
| -rw-r--r-- | Eigen/src/Cholesky/LDLT.h | 14 | ||||
| -rw-r--r-- | test/cholesky.cpp | 32 |
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 |