Fix bug #608: the sign computation in LDLT was broken

(transplanted from a69b4b092beec7caf87f1fe13cbba0781c651852
)
This commit is contained in:
Gael Guennebaud 2013-06-09 23:19:32 +02:00
parent 8f67e02ee2
commit 97c08b43b4

View File

@ -277,15 +277,13 @@ template<> struct ldlt_inplace<Lower>
// are compared; if any diagonal is negligible compared // are compared; if any diagonal is negligible compared
// to the largest overall, the algorithm bails. // to the largest overall, the algorithm bails.
cutoff = abs(NumTraits<Scalar>::epsilon() * biggest_in_corner); cutoff = abs(NumTraits<Scalar>::epsilon() * biggest_in_corner);
if(sign)
*sign = real(mat.diagonal().coeff(index_of_biggest_in_corner)) > 0 ? 1 : -1;
} }
// Finish early if the matrix is not full rank. // Finish early if the matrix is not full rank.
if(biggest_in_corner < cutoff) if(biggest_in_corner < cutoff)
{ {
for(Index i = k; i < size; i++) transpositions.coeffRef(i) = i; for(Index i = k; i < size; i++) transpositions.coeffRef(i) = i;
if(sign) *sign = 0;
break; break;
} }
@ -326,6 +324,16 @@ template<> struct ldlt_inplace<Lower>
} }
if((rs>0) && (abs(mat.coeffRef(k,k)) > cutoff)) if((rs>0) && (abs(mat.coeffRef(k,k)) > cutoff))
A21 /= mat.coeffRef(k,k); A21 /= mat.coeffRef(k,k);
if(sign)
{
// LDLT is not guaranteed to work for indefinite matrices, but let's try to get the sign right
int newSign = real(mat.diagonal().coeff(index_of_biggest_in_corner)) > 0;
if(k == 0)
*sign = newSign;
else if(*sign != newSign)
*sign = 0;
}
} }
return true; return true;