fix: issue 2481: LDLT produce wrong results with AutoDiffScalar

This commit is contained in:
Mario Rincon-Nigro 2022-05-25 15:26:10 +00:00 committed by Antonio Sánchez
parent 477eb7f630
commit e99163e732
3 changed files with 30 additions and 2 deletions

View File

@ -78,7 +78,7 @@ struct triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheLeft, Mode, Con
if (k>0)
rhs[i] -= (cjLhs.row(i).segment(s,k).transpose().cwiseProduct(Map<const Matrix<RhsScalar,Dynamic,1> >(rhs+s,k))).sum();
if((!(Mode & UnitDiag)) && numext::not_equal_strict(rhs[i],RhsScalar(0)))
if((!(Mode & UnitDiag)) && !is_identically_zero(rhs[i]))
rhs[i] /= cjLhs(i,i);
}
}
@ -115,7 +115,7 @@ struct triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheLeft, Mode, Con
for(Index k=0; k<actualPanelWidth; ++k)
{
Index i = IsLower ? pi+k : pi-k-1;
if(numext::not_equal_strict(rhs[i],RhsScalar(0)))
if(!is_identically_zero(rhs[i]))
{
if(!(Mode & UnitDiag))
rhs[i] /= cjLhs.coeff(i,i);

View File

@ -441,6 +441,17 @@ bool not_equal_strict(const double& x,const double& y) { return std::not_equal_t
} // end namespace numext
namespace internal {
template<typename Scalar>
struct is_identically_zero_impl {
static inline bool run(const Scalar& s) {
return numext::is_exactly_zero(s);
}
};
template<typename Scalar> EIGEN_STRONG_INLINE
bool is_identically_zero(const Scalar& s) { return is_identically_zero_impl<Scalar>::run(s); }
/// \internal Returns true if its argument is of integer or enum type.
/// FIXME this has the same purpose as `is_valid_index_type` in XprHelper.h
template<typename A>

View File

@ -715,6 +715,23 @@ template<typename DerType> struct NumTraits<AutoDiffScalar<DerType> >
};
};
namespace internal {
template<typename DerivativeType>
struct is_identically_zero_impl<AutoDiffScalar<DerivativeType>> {
static inline bool run(const AutoDiffScalar<DerivativeType>& s)
{
const DerivativeType& derivatives = s.derivatives();
for(int i=0; i<derivatives.size(); ++i)
{
if(!numext::is_exactly_zero(derivatives[i]))
{
return false;
}
}
return numext::is_exactly_zero(s.value());
}
};
}
}
namespace std {