diff --git a/src/linalg/schur.rs b/src/linalg/schur.rs index 9dee67a49..bbba07228 100644 --- a/src/linalg/schur.rs +++ b/src/linalg/schur.rs @@ -348,12 +348,15 @@ where DefaultAllocator: Allocator>, { let mut n = end; + 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 +375,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()) { diff --git a/tests/linalg/eigen.rs b/tests/linalg/eigen.rs index 2c43e11ec..ac775b35e 100644 --- a/tests/linalg/eigen.rs +++ b/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::>().re, + m.trace(), + epsilon = 1e-10 + ); + + assert_relative_eq!( + complex_eigenvals.iter().product::>().re, + m.determinant(), + epsilon = 1e-10 + ); +} + // #[cfg(feature = "arbitrary")] // quickcheck! { // TODO: full eigendecomposition is not implemented yet because of its complexity when some @@ -225,4 +254,4 @@ fn very_small_deviation_from_identity_issue_1368() { // println!("{}{:.5}{:.5}", m, mv, eig.eigenvectors); // // relative_eq!(eig.eigenvectors, mv, epsilon = 1.0e-5) -// } +// } \ No newline at end of file