add exceptional shift (#1603)

* run cargo fmt

* add comment on exceptional shift

* Update src/linalg/schur.rs

Co-authored-by: Alex H. Room <69592136+alexhroom@users.noreply.github.com>

---------

Co-authored-by: Alex H. Room <69592136+alexhroom@users.noreply.github.com>

Author: Toucan4Life

src/linalg/schur.rs

@@ -135,6 +135,16 @@ where
 
         // Implicit double-shift QR method.
         let mut niter = 0;
+        // Number of QR steps since the active window last shrank due to a deflation.
+        let mut kdefl: usize = 0;
+        // Apply Wilkinson's exceptional shift every kexsh stalled QR steps.
+        let kexsh: usize = 10;
+        // This is used by the exceptional shift algorithm as mentioned in
+        // Wilkinson, J. H., & Reinsch, C. (1971).
+        // Handbook for automatic computation (Vol. 2: Linear algebra)
+        // P.373
+        let exceptional_shift_a: T::RealField = crate::convert(0.75);
+        let exceptional_shift_b: T::RealField = crate::convert(-0.4375);
         let (mut start, mut end) = Self::delimit_subproblem(&mut t, eps.clone(), dim.value() - 1);
 
         while end != start {
@@ -150,10 +160,28 @@ where
                 let h22 = t[(start + 1, start + 1)].clone();
                 let h32 = t[(start + 2, start + 1)].clone();
 
-                let hnn = t[(n, n)].clone();
-                let hmm = t[(m, m)].clone();
-                let hnm = t[(n, m)].clone();
-                let hmn = t[(m, n)].clone();
+                let (hnn, hmm, hnm, hmn) = if kdefl != 0 && kdefl % kexsh == 0 {
+                    let s = t[(start + 1, start)].clone().norm1()
+                        + t[(start + 2, start + 1)].clone().norm1();
+                    let h = t[(start, start)].clone()
+                        + T::from_real(exceptional_shift_a.clone() * s.clone());
+                    // When the implicit QR step stalls, replace the trailing 2x2 Wilkinson shift
+                    // with the exceptional-shift coefficients from Wilkinson & Reinsch. This
+                    // perturbs the bulge start enough to escape the non-convergent cycle.
+                    (
+                        h.clone(),
+                        h,
+                        T::from_real(exceptional_shift_b.clone() * s.clone()),
+                        T::from_real(s),
+                    )
+                } else {
+                    (
+                        t[(n, n)].clone(),
+                        t[(m, m)].clone(),
+                        t[(n, m)].clone(),
+                        t[(m, n)].clone(),
+                    )
+                };
 
                 let tra = hnn.clone() + hmm.clone();
                 let det = hnn * hmm - hnm * hmn;
@@ -256,10 +284,19 @@ where
                 }
             }
 
+            let prev_start = start;
+            let prev_end = end;
             let sub = Self::delimit_subproblem(&mut t, eps.clone(), end);
 
             start = sub.0;
             end = sub.1;
+            if end < prev_end || start > prev_start {
+                // A split or deflation changed the active window, so restart the stall counter.
+                kdefl = 0;
+            } else {
+                // No deflation occurred: count another QR step toward an exceptional shift.
+                kdefl += 1;
+            }
 
             niter += 1;
             if niter == max_niter {

tests/linalg/eigen.rs

@@ -170,6 +170,31 @@ fn eigenvalues_search_should_not_hang_issue_1528() {
     );
 }
 
+// Regression test for #611
+#[test]
+#[rustfmt::skip]
+fn eigenvalues_search_should_not_hang_issue_611() {
+    let m = nalgebra::Matrix4::<f64>::new(
+        0.0, 0.0, 0.0, -0.8286,
+        1.0, 0.0, 0.0,  0.0,
+        0.0, 1.0, 0.0,  1.7094,
+        0.0, 0.0, 1.0,  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