Implement absolute threshold for early deflation in Schur algorithm (#1565)

* Implement absolute threshold for early deflation in Schur algorithm

* Implement absolute threshold for early deflation in Schur algorithm

* fix formatting

* added comment

Author: Toucan4Life

src/linalg/schur.rs

@@ -348,12 +348,18 @@ where
         DefaultAllocator: Allocator<DimDiff<D, U1>>,
     {
         let mut n = end;
+        // Equivalent to SMLNUM in LAPACK DLAHQR. Since SMLNUM depends on the machine
+        // precision, we use eps^2 here as best approximation.
+        // This is justified because the relative threshold use eps * ( t[n,n] + t[m,m] )
+        let absolute_threshold = eps.clone() * eps.clone();
 
         while n > 0 {
             let m = n - 1;
+            let off_diag_norm = t[(n, m)].clone().norm1();
 
-            if t[(n, m)].clone().norm1()
-                <= eps.clone() * (t[(n, n)].clone().norm1() + t[(m, m)].clone().norm1())
+            if off_diag_norm <= absolute_threshold
+                || off_diag_norm
+                    <= eps.clone() * (t[(n, n)].clone().norm1() + t[(m, m)].clone().norm1())
             {
                 t[(n, m)] = T::zero();
             } else {
@@ -372,8 +378,11 @@ where
             let m = new_start - 1;
 
             let off_diag = t[(new_start, m)].clone();
+            let off_diag_norm = off_diag.clone().norm1();
+
             if off_diag.is_zero()
-                || off_diag.norm1()
+                || off_diag_norm <= absolute_threshold
+                || off_diag_norm
                     <= eps.clone()
                         * (t[(new_start, new_start)].clone().norm1() + t[(m, m)].clone().norm1())
             {

tests/linalg/eigen.rs

@@ -1,4 +1,4 @@
-use na::{DMatrix, Matrix3};
+use na::{Complex, DMatrix, Matrix3};
 
 #[cfg(feature = "proptest-support")]
 mod proptest_tests {
@@ -141,6 +141,35 @@ fn very_small_deviation_from_identity_issue_1368() {
     }
 }
 
+// Regression test for #1528
+#[test]
+#[rustfmt::skip]
+fn eigenvalues_search_should_not_hang_issue_1528() {
+    let m = DMatrix::from_row_slice(
+        4,
+        4,
+        &[
+            0.0_f64,            0.5773502691896257, 0.0,                0.0,
+            0.5773502691896257, 0.0,                0.5163977794943222, 0.0,
+            0.0,                0.5163977794943222, 0.0,                0.50709255283711,
+            0.0,                0.0,                0.50709255283711,   0.0,
+        ],
+    );
+    let complex_eigenvals = m.complex_eigenvalues();
+
+    assert_relative_eq!(
+        complex_eigenvals.iter().sum::<Complex<f64>>().re,
+        m.trace(),
+        epsilon = 1e-10
+    );
+
+    assert_relative_eq!(
+        complex_eigenvals.iter().product::<Complex<f64>>().re,
+        m.determinant(),
+        epsilon = 1e-10
+    );
+}
+
 //  #[cfg(feature = "arbitrary")]
 //  quickcheck! {
 // TODO: full eigendecomposition is not implemented yet because of its complexity when some