SVD

src/binding.rs

@@ -22,6 +22,20 @@ pub trait LapackBinding: Sized {
               lwork: i32,
               info: &mut i32);
     fn _lange(norm: u8, m: i32, n: i32, a: &Vec<Self>, lda: i32, work: &mut Vec<Self>) -> Self;
+    fn _gesvd(jobu: u8,
+              jobvt: u8,
+              m: i32,
+              n: i32,
+              a: &mut [Self],
+              lda: i32,
+              s: &mut [Self],
+              u: &mut [Self],
+              ldu: i32,
+              vt: &mut [Self],
+              ldvt: i32,
+              work: &mut [Self],
+              lwork: i32,
+              info: &mut i32);
 }
 
 impl LapackBinding for f64 {
@@ -51,6 +65,35 @@ impl LapackBinding for f64 {
     fn _lange(norm: u8, m: i32, n: i32, a: &Vec<Self>, lda: i32, work: &mut Vec<Self>) -> Self {
         dlange(norm, m, n, a, lda, work)
     }
+    fn _gesvd(jobu: u8,
+              jobvt: u8,
+              m: i32,
+              n: i32,
+              a: &mut [Self],
+              lda: i32,
+              s: &mut [Self],
+              u: &mut [Self],
+              ldu: i32,
+              vt: &mut [Self],
+              ldvt: i32,
+              work: &mut [Self],
+              lwork: i32,
+              info: &mut i32) {
+        dgesvd(jobu,
+               jobvt,
+               m,
+               n,
+               a,
+               lda,
+               s,
+               u,
+               ldu,
+               vt,
+               ldvt,
+               work,
+               lwork,
+               info);
+    }
 }
 
 impl LapackBinding for f32 {
@@ -80,4 +123,33 @@ impl LapackBinding for f32 {
     fn _lange(norm: u8, m: i32, n: i32, a: &Vec<Self>, lda: i32, work: &mut Vec<Self>) -> Self {
         slange(norm, m, n, a, lda, work)
     }
+    fn _gesvd(jobu: u8,
+              jobvt: u8,
+              m: i32,
+              n: i32,
+              a: &mut [Self],
+              lda: i32,
+              s: &mut [Self],
+              u: &mut [Self],
+              ldu: i32,
+              vt: &mut [Self],
+              ldvt: i32,
+              work: &mut [Self],
+              lwork: i32,
+              info: &mut i32) {
+        sgesvd(jobu,
+               jobvt,
+               m,
+               n,
+               a,
+               lda,
+               s,
+               u,
+               ldu,
+               vt,
+               ldvt,
+               work,
+               lwork,
+               info);
+    }
 }

src/matrix.rs

@@ -1,6 +1,7 @@
 
 use ndarray::prelude::*;
 
+use error::LapackError;
 use scalar::LapackScalar;
 
 pub trait Matrix: Sized {
@@ -11,7 +12,7 @@ pub trait Matrix: Sized {
     fn norm_1(&self) -> Self::Scalar;
     fn norm_i(&self) -> Self::Scalar;
     fn norm_f(&self) -> Self::Scalar;
-    // fn svd(self) -> (Self, Self::Vector, Self);
+    fn svd(self) -> Result<(Self, Self::Vector, Self), LapackError>;
 }
 
 impl<A: LapackScalar> Matrix for Array<A, (Ix, Ix)> {
@@ -42,4 +43,24 @@ impl<A: LapackScalar> Matrix for Array<A, (Ix, Ix)> {
         let (m, n) = self.size();
         LapackScalar::norm_f(m, n, self.clone().into_raw_vec())
     }
+    fn svd(self) -> Result<(Self, Self::Vector, Self), LapackError> {
+        let strides = self.strides();
+        let (m, n) = if strides[0] > strides[1] {
+            self.size()
+        } else {
+            let (n, m) = self.size();
+            (m, n)
+        };
+        let (u, s, vt) = try!(LapackScalar::svd(m, n, self.clone().into_raw_vec()));
+        let sv = Array::from_vec(s);
+        if strides[0] > strides[1] {
+            let ua = Array::from_vec(u).into_shape((n, n)).unwrap();
+            let va = Array::from_vec(vt).into_shape((m, m)).unwrap();
+            Ok((va, sv, ua))
+        } else {
+            let ua = Array::from_vec(u).into_shape((n, n)).unwrap().reversed_axes();
+            let va = Array::from_vec(vt).into_shape((m, m)).unwrap().reversed_axes();
+            Ok((ua, sv, va))
+        }
+    }
 }

src/scalar.rs

@@ -1,10 +1,29 @@
 
+use std::cmp::min;
+use std::fmt::Debug;
 use ndarray::LinalgScalar;
+use num_traits::float::Float;
 
 use error::LapackError;
 use binding;
 
-pub trait LapackScalar: LinalgScalar + binding::LapackBinding {
+pub trait TruncatableFloat: Float {
+    fn to_int(self) -> i32;
+}
+
+impl TruncatableFloat for f64 {
+    fn to_int(self) -> i32 {
+        self as i32
+    }
+}
+impl TruncatableFloat for f32 {
+    fn to_int(self) -> i32 {
+        self as i32
+    }
+}
+
+pub trait LapackScalar
+    : Debug + TruncatableFloat + LinalgScalar + binding::LapackBinding {
     fn eigh(n: usize, mut a: Vec<Self>) -> Result<(Vec<Self>, Vec<Self>), LapackError> {
         let mut w = vec![Self::zero(); n];
         let mut work = vec![Self::zero(); 4 * n];
@@ -54,6 +73,57 @@ pub trait LapackScalar: LinalgScalar + binding::LapackBinding {
         let mut work = Vec::<Self>::new();
         Self::_lange(b'f', m as i32, n as i32, &mut a, m as i32, &mut work)
     }
+    fn svd(n: usize,
+           m: usize,
+           mut a: Vec<Self>)
+           -> Result<(Vec<Self>, Vec<Self>, Vec<Self>), LapackError> {
+        let mut info = 0;
+        let n = n as i32;
+        let m = m as i32;
+        let lda = m;
+        let ldu = m;
+        let ldvt = n;
+        let lwmax = 1000; // XXX
+        let lwork = -1;
+        let mut u = vec![Self::zero(); (ldu * m) as usize];
+        let mut vt = vec![Self::zero(); (ldvt * n) as usize];
+        let mut s = vec![Self::zero(); n as usize];
+        let mut work = vec![Self::zero(); lwmax];
+        Self::_gesvd('A' as u8,
+                     'A' as u8,
+                     m,
+                     n,
+                     &mut a,
+                     lda,
+                     &mut s,
+                     &mut u,
+                     ldu,
+                     &mut vt,
+                     ldvt,
+                     &mut work,
+                     lwork,
+                     &mut info); // calc optimal work
+        let lwork = min(lwmax as i32, work[0].to_int());
+        Self::_gesvd('A' as u8,
+                     'A' as u8,
+                     m,
+                     n,
+                     &mut a,
+                     lda,
+                     &mut s,
+                     &mut u,
+                     ldu,
+                     &mut vt,
+                     ldvt,
+                     &mut work,
+                     lwork,
+                     &mut info);
+        if info == 0 {
+            Ok((u, s, vt))
+        } else {
+            Err(From::from(info))
+        }
+    }
 }
 
 impl LapackScalar for f64 {}

tests/svd.rs

@@ -0,0 +1,87 @@
+
+extern crate rand;
+extern crate ndarray;
+extern crate ndarray_rand;
+extern crate ndarray_linalg;
+
+use ndarray::prelude::*;
+use ndarray_linalg::prelude::*;
+use rand::distributions::*;
+use ndarray_rand::RandomExt;
+
+fn all_close(a: Array<f64, (Ix, Ix)>, b: Array<f64, (Ix, Ix)>) {
+    if !a.all_close(&b, 1.0e-7) {
+        panic!("\nTwo matrices are not equal:\na = \n{:?}\nb = \n{:?}\n",
+               a,
+               b);
+    }
+}
+
+#[test]
+fn svd_square() {
+    let r_dist = Range::new(0., 1.);
+    let a = Array::<f64, _>::random((3, 3), r_dist);
+    let (u, s, vt) = a.clone().svd().unwrap();
+    let mut sm = Array::eye(3);
+    for i in 0..3 {
+        sm[(i, i)] = s[i];
+    }
+    all_close(u.dot(&sm).dot(&vt), a);
+}
+#[test]
+fn svd_square_t() {
+    let r_dist = Range::new(0., 1.);
+    let a = Array::<f64, _>::random((3, 3), r_dist).reversed_axes();
+    let (u, s, vt) = a.clone().svd().unwrap();
+    let mut sm = Array::eye(3);
+    for i in 0..3 {
+        sm[(i, i)] = s[i];
+    }
+    all_close(u.dot(&sm).dot(&vt), a);
+}
+
+#[test]
+fn svd_4x3() {
+    let r_dist = Range::new(0., 1.);
+    let a = Array::<f64, _>::random((4, 3), r_dist);
+    let (u, s, vt) = a.clone().svd().unwrap();
+    let mut sm = Array::zeros((4, 3));
+    for i in 0..3 {
+        sm[(i, i)] = s[i];
+    }
+    all_close(u.dot(&sm).dot(&vt), a);
+}
+#[test]
+fn svd_4x3_t() {
+    let r_dist = Range::new(0., 1.);
+    let a = Array::<f64, _>::random((3, 4), r_dist).reversed_axes();
+    let (u, s, vt) = a.clone().svd().unwrap();
+    let mut sm = Array::zeros((4, 3));
+    for i in 0..3 {
+        sm[(i, i)] = s[i];
+    }
+    all_close(u.dot(&sm).dot(&vt), a);
+}
+
+#[test]
+fn svd_3x4() {
+    let r_dist = Range::new(0., 1.);
+    let a = Array::<f64, _>::random((3, 4), r_dist);
+    let (u, s, vt) = a.clone().svd().unwrap();
+    let mut sm = Array::zeros((3, 4));
+    for i in 0..3 {
+        sm[(i, i)] = s[i];
+    }
+    all_close(u.dot(&sm).dot(&vt), a);
+}
+#[test]
+fn svd_3x4_t() {
+    let r_dist = Range::new(0., 1.);
+    let a = Array::<f64, _>::random((4, 3), r_dist).reversed_axes();
+    let (u, s, vt) = a.clone().svd().unwrap();
+    let mut sm = Array::zeros((3, 4));
+    for i in 0..3 {
+        sm[(i, i)] = s[i];
+    }
+    all_close(u.dot(&sm).dot(&vt), a);
+}