Fix SymmetricEigen Routine (#1210)

* Add known-failing unit test case for SymmetricEigen

* Fix eigenvector calculation for subdim=2 case in SymmetricEigen

Author: CattleProdigy

src/linalg/symmetric_eigen.rs

@@ -205,10 +205,21 @@ where
                     diag[start + 1].clone(),
                 );
                 let eigvals = m.eigenvalues().unwrap();
-                let basis = Vector2::new(
-                    eigvals.x.clone() - diag[start + 1].clone(),
-                    off_diag[start].clone(),
-                );
+
+                // Choose the basis least likely to experience cancellation
+                let basis = if (eigvals.x.clone() - diag[start + 1].clone()).abs()
+                    > (eigvals.x.clone() - diag[start].clone()).abs()
+                {
+                    Vector2::new(
+                        eigvals.x.clone() - diag[start + 1].clone(),
+                        off_diag[start].clone(),
+                    )
+                } else {
+                    Vector2::new(
+                        off_diag[start].clone(),
+                        eigvals.x.clone() - diag[start].clone(),
+                    )
+                };
 
                 diag[start] = eigvals[0].clone();
                 diag[start + 1] = eigvals[1].clone();
@@ -348,7 +359,36 @@ where
 
 #[cfg(test)]
 mod test {
-    use crate::base::Matrix2;
+    use crate::base::{Matrix2, Matrix4};
+
+    /// Exercises bug reported in issue #1109 of nalgebra (https://github.com/dimforge/nalgebra/issues/1109)
+    #[test]
+    fn symmetric_eigen_regression_issue_1109() {
+        let m = Matrix4::new(
+            -19884.07f64,
+            -10.07188,
+            11.277279,
+            -188560.63,
+            -10.07188,
+            12.518197,
+            1.3770627,
+            -102.97504,
+            11.277279,
+            1.3770627,
+            14.587362,
+            113.26099,
+            -188560.63,
+            -102.97504,
+            113.26099,
+            -1788112.3,
+        );
+        let eig = m.symmetric_eigen();
+        assert!(relative_eq!(
+            m.lower_triangle(),
+            eig.recompose().lower_triangle(),
+            epsilon = 1.0e-5
+        ));
+    }
 
     fn expected_shift(m: Matrix2<f64>) -> f64 {
         let vals = m.eigenvalues().unwrap();