Use layout at SVD

Author: termoshtt

src/matrix.rs

@@ -89,23 +89,15 @@ impl<A> Matrix for Array<A, Ix2>
         ImplNorm::norm_f(m, n, self.clone().into_raw_vec())
     }
     fn svd(self) -> Result<(Self, Self::Vector, Self), LinalgError> {
-        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!(ImplSVD::svd(m, n, self.clone().into_raw_vec()));
+        let (n, m) = self.size();
+        let layout = self.layout()?;
+        let (u, s, vt) = ImplSVD::svd(layout, 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))
+        let ua = Array::from_vec(u).into_shape((n, n)).unwrap();
+        let va = Array::from_vec(vt).into_shape((m, m)).unwrap();
+        match layout {
+            Layout::RowMajor => Ok((ua, sv, va)), 
+            Layout::ColumnMajor => Ok((ua.reversed_axes(), sv, va.reversed_axes())),
         }
     }
     fn qr(self) -> Result<(Self, Self), LinalgError> {

src/svd.rs

@@ -1,65 +1,37 @@
 //! Implement SVD
 
-use lapack::fortran::*;
+use std::cmp::min;
+use lapack::c::*;
 use num_traits::Zero;
 
 use error::LapackError;
 
 pub trait ImplSVD: Sized {
-    fn svd(n: usize, m: usize, mut a: Vec<Self>) -> Result<(Vec<Self>, Vec<Self>, Vec<Self>), LapackError>;
+    fn svd(layout: Layout,
+           n: usize,
+           m: usize,
+           mut a: Vec<Self>)
+           -> Result<(Vec<Self>, Vec<Self>, Vec<Self>), LapackError>;
 }
 
 macro_rules! impl_svd {
     ($scalar:ty, $gesvd:path) => {
 impl ImplSVD for $scalar {
-    fn svd(n: usize,
-           m: usize,
-           mut a: Vec<Self>)
-           -> Result<(Vec<Self>, Vec<Self>, Vec<Self>), LapackError> {
-        let mut info = 0;
+    fn svd(layout: Layout, n: usize, m: usize, mut a: Vec<Self>) -> Result<(Vec<Self>, Vec<Self>, Vec<Self>), LapackError> {
+        let k = min(n, m);
         let n = n as i32;
         let m = m as i32;
-        let lda = m;
+        let lda = match layout {
+            Layout::RowMajor => n,
+            Layout::ColumnMajor => m,
+        };
         let ldu = m;
         let ldvt = n;
-        let lwork = -1;
-        let lw_default = 1000;
         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(); lw_default];
-        $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 = work[0] as i32;
-        if lwork > lw_default as i32 {
-            work = vec![Self::zero(); lwork as usize];
-        }
-        $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);
+        let mut superb = vec![Self::zero(); k-2];
+        let info = $gesvd(layout, 'A' as u8, 'A' as u8, m, n, &mut a, lda, &mut s, &mut u, ldu, &mut vt, ldvt, &mut superb);
         if info == 0 {
             Ok((u, s, vt))
         } else {