mirror of
https://gitlab.com/libeigen/eigen.git
synced 2025-06-04 18:54:00 +08:00
Fix bug #608: the sign computation in LDLT was broken
This commit is contained in:
parent
c98fd7a6ca
commit
a69b4b092b
@ -278,22 +278,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;
|
|
||||||
}
|
|
||||||
else 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(newSign != *sign)
|
|
||||||
*sign = 0;
|
|
||||||
}
|
}
|
||||||
|
|
||||||
// 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;
|
||||||
}
|
}
|
||||||
|
|
||||||
@ -334,6 +325,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;
|
||||||
|
Loading…
x
Reference in New Issue
Block a user