Merge pull request #163 from rust-ndarray/bugfix_svd_162

Fix SVD without u or vt

Author: termoshtt

Cargo.toml

@@ -46,3 +46,6 @@ version = "0.6"
 default-features = false
 features = ["static"]
 optional = true
+
+[dev-dependencies]
+paste = "0.1"

src/lapack/svd.rs

@@ -42,15 +42,16 @@ macro_rules! impl_svd {
                 let (ju, ldu, mut u) = if calc_u {
                     (FlagSVD::All, m, vec![Self::zero(); (m * m) as usize])
                 } else {
-                    (FlagSVD::No, 0, Vec::new())
+                    (FlagSVD::No, 1, Vec::new())
                 };
                 let (jvt, ldvt, mut vt) = if calc_vt {
                     (FlagSVD::All, n, vec![Self::zero(); (n * n) as usize])
                 } else {
-                    (FlagSVD::No, 0, Vec::new())
+                    (FlagSVD::No, n, Vec::new())
                 };
                 let mut s = vec![Self::Real::zero(); k as usize];
                 let mut superb = vec![Self::Real::zero(); (k - 1) as usize];
+                dbg!(ldvt);
                 let info = $gesvd(
                     l.lapacke_layout(),
                     ju as u8,
@@ -70,8 +71,8 @@ macro_rules! impl_svd {
                     info,
                     SVDOutput {
                         s: s,
-                        u: if ldu > 0 { Some(u) } else { None },
-                        vt: if ldvt > 0 { Some(vt) } else { None },
+                        u: if calc_u { Some(u) } else { None },
+                        vt: if calc_vt { Some(vt) } else { None },
                     },
                 )
             }

src/svd.rs

@@ -75,8 +75,12 @@ where
         let l = self.layout()?;
         let svd_res = unsafe { A::svd(l, calc_u, calc_vt, self.as_allocated_mut()?)? };
         let (n, m) = l.size();
-        let u = svd_res.u.map(|u| into_matrix(l.resized(n, n), u).unwrap());
-        let vt = svd_res.vt.map(|vt| into_matrix(l.resized(m, m), vt).unwrap());
+        let u = svd_res
+            .u
+            .map(|u| into_matrix(l.resized(n, n), u).expect("Size of U mismatches"));
+        let vt = svd_res
+            .vt
+            .map(|vt| into_matrix(l.resized(m, m), vt).expect("Size of VT mismatches"));
         let s = ArrayBase::from_vec(svd_res.s);
         Ok((u, s, vt))
     }

tests/svd.rs

@@ -2,7 +2,8 @@ use ndarray::*;
 use ndarray_linalg::*;
 use std::cmp::min;
 
-fn test(a: &Array2<f64>, n: usize, m: usize) {
+fn test(a: &Array2<f64>) {
+    let (n, m) = a.dim();
     let answer = a.clone();
     println!("a = \n{:?}", a);
     let (u, s, vt): (_, Array1<_>, _) = a.svd(true, true).unwrap();
@@ -18,38 +19,62 @@ fn test(a: &Array2<f64>, n: usize, m: usize) {
     assert_close_l2!(&u.dot(&sm).dot(&vt), &answer, 1e-7);
 }
 
-#[test]
-fn svd_square() {
-    let a = random((3, 3));
-    test(&a, 3, 3);
+fn test_no_vt(a: &Array2<f64>) {
+    let (n, _m) = a.dim();
+    println!("a = \n{:?}", a);
+    let (u, _s, vt): (_, Array1<_>, _) = a.svd(true, false).unwrap();
+    assert!(u.is_some());
+    assert!(vt.is_none());
+    let u = u.unwrap();
+    assert_eq!(u.dim().0, n);
+    assert_eq!(u.dim().1, n);
 }
 
-#[test]
-fn svd_square_t() {
-    let a = random((3, 3).f());
-    test(&a, 3, 3);
+fn test_no_u(a: &Array2<f64>) {
+    let (_n, m) = a.dim();
+    println!("a = \n{:?}", a);
+    let (u, _s, vt): (_, Array1<_>, _) = a.svd(false, true).unwrap();
+    assert!(u.is_none());
+    assert!(vt.is_some());
+    let vt = vt.unwrap();
+    assert_eq!(vt.dim().0, m);
+    assert_eq!(vt.dim().1, m);
 }
 
-#[test]
-fn svd_3x4() {
-    let a = random((3, 4));
-    test(&a, 3, 4);
+fn test_diag_only(a: &Array2<f64>) {
+    println!("a = \n{:?}", a);
+    let (u, _s, vt): (_, Array1<_>, _) = a.svd(false, false).unwrap();
+    assert!(u.is_none());
+    assert!(vt.is_none());
 }
 
-#[test]
-fn svd_3x4_t() {
-    let a = random((3, 4).f());
-    test(&a, 3, 4);
-}
+macro_rules! test_svd_impl {
+    ($test:ident, $n:expr, $m:expr) => {
+        paste::item! {
+            #[test]
+            fn [<svd_ $test _ $n x $m>]() {
+                let a = random(($n, $m));
+                $test(&a);
+            }
 
-#[test]
-fn svd_4x3() {
-    let a = random((4, 3));
-    test(&a, 4, 3);
+            #[test]
+            fn [<svd_ $test _ $n x $m _t>]() {
+                let a = random(($n, $m).f());
+                $test(&a);
+            }
+        }
+    };
 }
 
-#[test]
-fn svd_4x3_t() {
-    let a = random((4, 3).f());
-    test(&a, 4, 3);
-}
+test_svd_impl!(test, 3, 3);
+test_svd_impl!(test_no_vt, 3, 3);
+test_svd_impl!(test_no_u, 3, 3);
+test_svd_impl!(test_diag_only, 3, 3);
+test_svd_impl!(test, 4, 3);
+test_svd_impl!(test_no_vt, 4, 3);
+test_svd_impl!(test_no_u, 4, 3);
+test_svd_impl!(test_diag_only, 4, 3);
+test_svd_impl!(test, 3, 4);
+test_svd_impl!(test_no_vt, 3, 4);
+test_svd_impl!(test_no_u, 3, 4);
+test_svd_impl!(test_diag_only, 3, 4);