From c44bbabdccc8269cfe395539c3aea85415583fab Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 12 Jul 2010 16:42:38 +0200 Subject: [PATCH] fix LU and QR solve when rank==0, fix LLT when the matrix is purely 0 --- Eigen/src/Cholesky/LLT.h | 2 +- Eigen/src/LU/LU.h | 7 ++++--- Eigen/src/QR/QR.h | 3 +++ 3 files changed, 8 insertions(+), 4 deletions(-) diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h index 42c959f83..b7bf58e63 100644 --- a/Eigen/src/Cholesky/LLT.h +++ b/Eigen/src/Cholesky/LLT.h @@ -136,7 +136,7 @@ void LLT::compute(const MatrixType& a) for (int j = 1; j < size; ++j) { x = ei_real(a.coeff(j,j)) - m_matrix.row(j).start(j).squaredNorm(); - if (x < cutoff) + if (x <= cutoff) { m_isPositiveDefinite = false; continue; diff --git a/Eigen/src/LU/LU.h b/Eigen/src/LU/LU.h index 862fa066c..3fb57ae25 100644 --- a/Eigen/src/LU/LU.h +++ b/Eigen/src/LU/LU.h @@ -515,9 +515,10 @@ bool LU::solve( if(!ei_isMuchSmallerThan(c.coeff(row,col), biggest_in_c, m_precision)) return false; } - m_lu.corner(TopLeft, m_rank, m_rank) - .template marked() - .solveTriangularInPlace(c.corner(TopLeft, m_rank, c.cols())); + if(m_rank>0) + m_lu.corner(TopLeft, m_rank, m_rank) + .template marked() + .solveTriangularInPlace(c.corner(TopLeft, m_rank, c.cols())); // Step 4 result->resize(m_lu.cols(), b.cols()); diff --git a/Eigen/src/QR/QR.h b/Eigen/src/QR/QR.h index 90751dd42..d96f4971a 100644 --- a/Eigen/src/QR/QR.h +++ b/Eigen/src/QR/QR.h @@ -276,6 +276,9 @@ bool QR::solve( // Q^T without explicitly forming matrixQ(). Investigate. *result = matrixQ().transpose()*b; + if(m_rank==0) + return result.isZero(); + if(!isSurjective()) { // is result is in the image of R ?