From e99163e73230335563006c067b6a3fc9ceca88af Mon Sep 17 00:00:00 2001
From: Mario Rincon-Nigro <mario.rincon.nigro@gmail.com>
Date: Wed, 25 May 2022 15:26:10 +0000
Subject: [PATCH] fix: issue 2481: LDLT produce wrong results with
 AutoDiffScalar

---
 .../src/Core/products/TriangularSolverVector.h  |  4 ++--
 Eigen/src/Core/util/Meta.h                      | 11 +++++++++++
 unsupported/Eigen/src/AutoDiff/AutoDiffScalar.h | 17 +++++++++++++++++
 3 files changed, 30 insertions(+), 2 deletions(-)

diff --git a/Eigen/src/Core/products/TriangularSolverVector.h b/Eigen/src/Core/products/TriangularSolverVector.h
index 57ade287e..b8fbb5bcc 100644
--- a/Eigen/src/Core/products/TriangularSolverVector.h
+++ b/Eigen/src/Core/products/TriangularSolverVector.h
@@ -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);
diff --git a/Eigen/src/Core/util/Meta.h b/Eigen/src/Core/util/Meta.h
index 8d5c2c3ee..32152ac40 100644
--- a/Eigen/src/Core/util/Meta.h
+++ b/Eigen/src/Core/util/Meta.h
@@ -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>
diff --git a/unsupported/Eigen/src/AutoDiff/AutoDiffScalar.h b/unsupported/Eigen/src/AutoDiff/AutoDiffScalar.h
index 47b0b344b..97222d10f 100644
--- a/unsupported/Eigen/src/AutoDiff/AutoDiffScalar.h
+++ b/unsupported/Eigen/src/AutoDiff/AutoDiffScalar.h
@@ -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 {