Eigenvalues: Make dense tridiagonal deflation scale invariant

libeigen/eigen!3000

Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h

diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
index 48dce9a..c79acca 100644
--- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
+++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
@@ -568,13 +568,21 @@
 
   using RealScalar = typename DiagType::RealScalar;
   const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
-  const RealScalar precision_inv = RealScalar(1) / NumTraits<RealScalar>::epsilon();
+  const RealScalar precision = NumTraits<RealScalar>::epsilon();
+  const RealScalar precision_inv = RealScalar(1) / precision;
 
   // Helper lambda for the deflation test.
   auto deflate = [&](Index lo, Index hi) {
     for (Index i = lo; i < hi; ++i) {
-      if (numext::abs(subdiag[i]) < considerAsZero) {
+      const RealScalar absSubdiag = numext::abs(subdiag[i]);
+      if (absSubdiag < considerAsZero) {
         subdiag[i] = RealScalar(0);
+      } else if (!PerBlockScaling) {
+        // Homogeneous in the global scale: both sides are linear in it, so prescaling cannot change which couplings
+        // are discarded (the previous test was quadratic on the left).
+        if (absSubdiag <= precision * numext::maxi(numext::abs(diag[i]), numext::abs(diag[i + 1]))) {
+          subdiag[i] = RealScalar(0);
+        }
       } else {
         const RealScalar scaled_subdiag = precision_inv * subdiag[i];
         if (scaled_subdiag * scaled_subdiag <= (numext::abs(diag[i]) + numext::abs(diag[i + 1]))) {

test/eigensolver_selfadjoint.cpp

diff --git a/test/eigensolver_selfadjoint.cpp b/test/eigensolver_selfadjoint.cpp
index 48c2d71..5f7a32a 100644
--- a/test/eigensolver_selfadjoint.cpp
+++ b/test/eigensolver_selfadjoint.cpp
@@ -486,6 +486,23 @@
   VERIFY_IS_APPROX(eig2.eigenvalues(), eig2v.eigenvalues());
 }
 
+template <typename RealScalar>
+void selfadjointeigensolver_dense_deflation_scale_invariance() {
+  const RealScalar epsilon = NumTraits<RealScalar>::epsilon();
+  const RealScalar scales[] = {RealScalar(1), RealScalar(2)};
+  for (RealScalar scale : scales) {
+    Matrix<RealScalar, 2, 1> diag = Matrix<RealScalar, 2, 1>::Constant(scale);
+    Matrix<RealScalar, Dynamic, 1> subdiag(1);
+    subdiag[0] = RealScalar(1.25) * epsilon * scale;
+    Matrix<std::complex<RealScalar>, Dynamic, Dynamic> eigenvectors =
+        Matrix<std::complex<RealScalar>, Dynamic, Dynamic>::Identity(2, 2);
+
+    const ComputationInfo info = internal::computeFromTridiagonal_impl<false>(diag, subdiag, 30, true, eigenvectors);
+    VERIFY_IS_EQUAL(info, Success);
+    VERIFY(diag[0] < diag[1]);
+  }
+}
+
 // Test computeFromTridiagonal with wide dynamic range across decoupled blocks.
 // This exercises the per-block scaling in computeFromTridiagonal_impl: a zero on the
 // subdiagonal decouples the matrix into blocks with vastly different scales. Global
@@ -967,6 +984,8 @@
   int s = 0;
   CALL_SUBTEST_4(generalizedselfadjointeigensolver_no_malloc<MatrixXd>());
   CALL_SUBTEST_5(generalizedselfadjointeigensolver_no_malloc<MatrixXcd>());
+  CALL_SUBTEST_5(selfadjointeigensolver_dense_deflation_scale_invariance<float>());
+  CALL_SUBTEST_5(selfadjointeigensolver_dense_deflation_scale_invariance<double>());
   CALL_SUBTEST_13(selfadjointeigensolver_subnormal_coefficients());
 
   for (int i = 0; i < g_repeat; i++) {