Document pseudo-inverse error conditions

src/base/matrix_view.rs

@@ -33,14 +33,6 @@ macro_rules! view_storage_impl (
         #[deprecated = "Use ViewStorage(Mut) instead."]
         pub type $legacy_name<'a, T, R, C, RStride, CStride> = $T<'a, T, R, C, RStride, CStride>;
 
-        unsafe impl<'a, T: Send, R: Dim, C: Dim, RStride: Dim, CStride: Dim> Send
-            for $T<'a, T, R, C, RStride, CStride>
-        {}
-
-        unsafe impl<'a, T: Sync, R: Dim, C: Dim, RStride: Dim, CStride: Dim> Sync
-            for $T<'a, T, R, C, RStride, CStride>
-        {}
-
         impl<'a, T, R: Dim, C: Dim, RStride: Dim, CStride: Dim> $T<'a, T, R, C, RStride, CStride> {
             /// Create a new matrix view without bounds checking and from a raw pointer.
             ///
@@ -136,6 +128,20 @@ impl<T: Scalar, R: Dim, C: Dim, RStride: Dim, CStride: Dim> Copy
 {
 }
 
+/// Safety: Equivalent to a shared reference to `T`. All `Dim` type arguments are `Send + Sync`. A
+/// shared reference can be sent iff `T: Sync`.
+unsafe impl<'a, T: Sync, R: Dim, C: Dim, RStride: Dim, CStride: Dim> Send
+    for ViewStorage<'a, T, R, C, RStride, CStride>
+{
+}
+
+/// Safety: Equivalent to a shared reference to `T`. All `Dim` type arguments are `Send + Sync`. A
+/// shared reference is `Sync` iff `T: Sync`.
+unsafe impl<'a, T: Sync, R: Dim, C: Dim, RStride: Dim, CStride: Dim> Sync
+    for ViewStorage<'a, T, R, C, RStride, CStride>
+{
+}
+
 impl<T: Scalar, R: Dim, C: Dim, RStride: Dim, CStride: Dim> Clone
     for ViewStorage<'_, T, R, C, RStride, CStride>
 {
@@ -145,6 +151,20 @@ impl<T: Scalar, R: Dim, C: Dim, RStride: Dim, CStride: Dim> Clone
     }
 }
 
+/// Safety: Equivalent to a unique reference to `T`. All `Dim` type arguments are `Send + Sync`. A
+/// unique reference is `Send` iff `T: Send`.
+unsafe impl<'a, T: Send, R: Dim, C: Dim, RStride: Dim, CStride: Dim> Send
+    for ViewStorageMut<'a, T, R, C, RStride, CStride>
+{
+}
+
+/// Safety: Equivalent to a unique reference to `T`. All `Dim` type arguments are `Send + Sync`. A
+/// unique reference is `Sync` iff `T: Sync`.
+unsafe impl<'a, T: Sync, R: Dim, C: Dim, RStride: Dim, CStride: Dim> Sync
+    for ViewStorageMut<'a, T, R, C, RStride, CStride>
+{
+}
+
 impl<'a, T: Scalar, R: Dim, C: Dim, RStride: Dim, CStride: Dim>
     ViewStorageMut<'a, T, R, C, RStride, CStride>
 where

src/linalg/svd.rs

@@ -584,7 +584,7 @@ where
     ///
     /// Any singular value smaller than `eps` is assumed to be zero.
     /// Returns `Err` if the right- and left- singular vectors have not
-    /// been computed at construction-time.
+    /// been computed at construction-time, or if `eps` is less than zero.
     pub fn pseudo_inverse(mut self, eps: T::RealField) -> Result<OMatrix<T, C, R>, &'static str>
     where
         DefaultAllocator: Allocator<C, R>,
@@ -818,6 +818,8 @@ where
     /// Computes the pseudo-inverse of this matrix.
     ///
     /// All singular values below `eps` are considered equal to 0.
+    ///
+    /// Returns `Err` if `eps` is less than zero.
     pub fn pseudo_inverse(self, eps: T::RealField) -> Result<OMatrix<T, C, R>, &'static str>
     where
         DefaultAllocator: Allocator<C, R>,

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();