Replace SVD

Author: termoshtt

src/impl2/mod.rs

@@ -7,5 +7,5 @@ pub use self::opnorm::*;
 pub use self::qr::*;
 pub use self::svd::*;
 
-pub trait LapackScalar: OperatorNorm_ + QR_ {}
-impl<A> LapackScalar for A where A: OperatorNorm_ + QR_ {}
+pub trait LapackScalar: OperatorNorm_ + QR_ + SVD_ {}
+impl<A> LapackScalar for A where A: OperatorNorm_ + QR_ + SVD_ {}

src/impl2/svd.rs

@@ -30,7 +30,7 @@ macro_rules! impl_svd {
 
 impl SVD_ for $scalar {
     fn svd(l: Layout, calc_u: bool, calc_vt: bool, mut a: &mut [Self]) -> Result<SVDOutput<Self>> {
-        let (n, m) = l.size();
+        let (m, n) = l.size();
         let k = ::std::cmp::min(n, m);
         let lda = l.lda();
         let (ju, ldu, mut u) = if calc_u {
@@ -39,7 +39,7 @@ impl SVD_ for $scalar {
             (FlagSVD::No, 0, Vec::new())
         };
         let (jvt, ldvt, mut vt) = if calc_vt {
-            (FlagSVD::All, m, vec![Self::zero(); (m*m) as usize])
+            (FlagSVD::All, n, vec![Self::zero(); (n*n) as usize])
         } else {
             (FlagSVD::No, 0, Vec::new())
         };

src/lib.rs

@@ -49,6 +49,7 @@ pub mod impls;
 pub mod impl2;
 
 pub mod qr;
+pub mod svd;
 pub mod opnorm;
 
 pub mod vector;

src/matrix.rs

@@ -21,8 +21,6 @@ pub trait Matrix: Sized {
     fn size(&self) -> (usize, usize);
     /// Layout (C/Fortran) of matrix
     fn layout(&self) -> Result<Layout, StrideError>;
-    /// singular-value decomposition (SVD)
-    fn svd(self) -> Result<(Self, Self::Vector, Self), LinalgError>;
     /// LU decomposition
     fn lu(self) -> Result<(Self::Permutator, Self, Self), LinalgError>;
     /// permutate matrix (inplace)
@@ -74,18 +72,6 @@ impl<A: MFloat> Matrix for Array<A, Ix2> {
     fn layout(&self) -> Result<Layout, StrideError> {
         check_layout(self.strides())
     }
-    fn svd(self) -> Result<(Self, Self::Vector, Self), LinalgError> {
-        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);
-        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 lu(self) -> Result<(Self::Permutator, Self, Self), LinalgError> {
         let (n, m) = self.size();
         let k = min(n, m);
@@ -129,10 +115,6 @@ impl<A: MFloat> Matrix for RcArray<A, Ix2> {
     fn layout(&self) -> Result<Layout, StrideError> {
         check_layout(self.strides())
     }
-    fn svd(self) -> Result<(Self, Self::Vector, Self), LinalgError> {
-        let (u, s, v) = self.into_owned().svd()?;
-        Ok((u.into_shared(), s.into_shared(), v.into_shared()))
-    }
     fn lu(self) -> Result<(Self::Permutator, Self, Self), LinalgError> {
         let (p, l, u) = self.into_owned().lu()?;
         Ok((p, l.into_shared(), u.into_shared()))

src/prelude.rs

@@ -7,4 +7,5 @@ pub use util::*;
 pub use assert::*;
 
 pub use qr::*;
+pub use svd::*;
 pub use opnorm::*;

src/svd.rs

@@ -0,0 +1,40 @@
+
+use ndarray::*;
+
+use super::error::*;
+use super::layout::{Layout, AllocatedArray, AllocatedArrayMut};
+use impl2::LapackScalar;
+
+pub trait SVD<U, S, VT> {
+    fn svd(self, calc_u: bool, calc_vt: bool) -> Result<(Option<U>, S, Option<VT>)>;
+}
+
+impl<A, S, Su, Svt, Ss> SVD<ArrayBase<Su, Ix2>, ArrayBase<Ss, Ix1>, ArrayBase<Svt, Ix2>> for ArrayBase<S, Ix2>
+    where A: LapackScalar,
+          S: DataMut<Elem = A>,
+          Su: DataOwned<Elem = A>,
+          Svt: DataOwned<Elem = A>,
+          Ss: DataOwned<Elem = A::Real>
+{
+    fn svd(mut self,
+           calc_u: bool,
+           calc_vt: bool)
+           -> Result<(Option<ArrayBase<Su, Ix2>>, ArrayBase<Ss, Ix1>, Option<ArrayBase<Svt, Ix2>>)> {
+        let n = self.rows();
+        let m = self.cols();
+        let l = self.layout()?;
+        let svd_res = A::svd(l, calc_u, calc_vt, self.as_allocated_mut()?)?;
+        let (u, vt) = match l {
+            Layout::C(_) => {
+                (svd_res.u.map(|u| ArrayBase::from_shape_vec((n, n), u).unwrap()),
+                 svd_res.vt.map(|vt| ArrayBase::from_shape_vec((m, m), vt).unwrap()))
+            }
+            Layout::F(_) => {
+                (svd_res.u.map(|u| ArrayBase::from_shape_vec((n, n).f(), u).unwrap()),
+                 svd_res.vt.map(|vt| ArrayBase::from_shape_vec((m, m).f(), vt).unwrap()))
+            }
+        };
+        let s = ArrayBase::from_vec(svd_res.s);
+        Ok((u, s, vt))
+    }
+}

tests/svd.rs

@@ -9,11 +9,13 @@ fn $funcname() {
     use ndarray_linalg::prelude::*;
     let a = $random($n, $m, $t);
     let answer = a.clone();
-    println!("a = \n{}", &a);
-    let (u, s, vt) = a.svd().unwrap();
-    println!("u = \n{}", &u);
-    println!("s = \n{}", &s);
-    println!("v = \n{}", &vt);
+    println!("a = \n{:?}", &a);
+    let (u, s, vt): (_, Array1<_>, _) = a.svd(true, true).unwrap();
+    let u: Array2<_> = u.unwrap();
+    let vt: Array2<_> = vt.unwrap();
+    println!("u = \n{:?}", &u);
+    println!("s = \n{:?}", &s);
+    println!("v = \n{:?}", &vt);
     let mut sm = Array::zeros(($n, $m));
     for i in 0..min($n, $m) {
         sm[(i, i)] = s[i];