stop fill pivoting LU only if the pivot is exactly 0

(transplanted from f278a3eabaa6553825daafbee0a92515c90b0b74
)
This commit is contained in:
Gael Guennebaud 2011-11-22 09:18:54 +01:00
parent 0ef41ec958
commit 9ed342a30e

View File

@ -443,7 +443,6 @@ FullPivLU<MatrixType>& FullPivLU<MatrixType>::compute(const MatrixType& matrix)
m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
m_maxpivot = RealScalar(0);
RealScalar cutoff(0);
for(Index k = 0; k < size; ++k)
{
@ -458,14 +457,7 @@ FullPivLU<MatrixType>& FullPivLU<MatrixType>::compute(const MatrixType& matrix)
row_of_biggest_in_corner += k; // correct the values! since they were computed in the corner,
col_of_biggest_in_corner += k; // need to add k to them.
// when k==0, biggest_in_corner is the biggest coeff absolute value in the original matrix
if(k == 0) cutoff = biggest_in_corner * NumTraits<Scalar>::epsilon();
// if the pivot (hence the corner) is "zero", terminate to avoid generating nan/inf values.
// Notice that using an exact comparison (biggest_in_corner==0) here, as Golub-van Loan do in
// their pseudo-code, results in numerical instability! The cutoff here has been validated
// by running the unit test 'lu' with many repetitions.
if(biggest_in_corner < cutoff)
if(biggest_in_corner==RealScalar(0))
{
// before exiting, make sure to initialize the still uninitialized transpositions
// in a sane state without destroying what we already have.