refactor div-conq SVD, impl svd_rand

Cargo.toml

@@ -41,6 +41,11 @@ default-features = false
 version = "0.3"
 default-features = false
 
+[dependencies.sprs]
+git = "https://github.com/pmarks/sprs.git"
+branch = "pmarks/scan-rs"
+optional = true
+
 [dependencies.openblas-src]
 version = "0.6"
 default-features = false

src/lapack/mod.rs

@@ -7,7 +7,6 @@ pub mod qr;
 pub mod solve;
 pub mod solveh;
 pub mod svd;
-pub mod svddc;
 pub mod triangular;
 
 pub use self::cholesky::*;
@@ -17,7 +16,6 @@ pub use self::qr::*;
 pub use self::solve::*;
 pub use self::solveh::*;
 pub use self::svd::*;
-pub use self::svddc::*;
 pub use self::triangular::*;
 
 use super::error::*;
@@ -26,7 +24,7 @@ use super::types::*;
 pub type Pivot = Vec<i32>;
 
 /// Trait for primitive types which implements LAPACK subroutines
-pub trait Lapack: OperatorNorm_ + QR_ + SVD_ + SVDDC_ + Solve_ + Solveh_ + Cholesky_ + Eigh_ + Triangular_ {}
+pub trait Lapack: OperatorNorm_ + QR_ + SVD_ + Solve_ + Solveh_ + Cholesky_ + Eigh_ + Triangular_ {}
 
 impl Lapack for f32 {}
 impl Lapack for f64 {}

src/lapack/svd.rs

@@ -9,13 +9,7 @@ use crate::types::*;
 
 use super::into_result;
 
-#[repr(u8)]
-enum FlagSVD {
-    All = b'A',
-    // OverWrite = b'O',
-    // Separately = b'S',
-    No = b'N',
-}
+use crate::svd::FlagSVD;
 
 /// Result of SVD
 pub struct SVDOutput<A: Scalar> {
@@ -27,35 +21,45 @@ pub struct SVDOutput<A: Scalar> {
     pub vt: Option<Vec<A>>,
 }
 
-/// Wraps `*gesvd`
+/// Wraps `*gesvd` and `*gesdd`
 pub trait SVD_: Scalar {
-    unsafe fn svd(l: MatrixLayout, calc_u: bool, calc_vt: bool, a: &mut [Self]) -> Result<SVDOutput<Self>>;
+    unsafe fn svd(l: MatrixLayout, jobu: FlagSVD, jobvt: FlagSVD, a: &mut [Self]) -> Result<SVDOutput<Self>>;
+    unsafe fn svd_dc(l: MatrixLayout, jobz: FlagSVD, a: &mut [Self]) -> Result<SVDOutput<Self>>;
 }
 
 macro_rules! impl_svd {
-    ($scalar:ty, $gesvd:path) => {
+    ($scalar:ty, $gesvd:path, $gesdd:path) => {
         impl SVD_ for $scalar {
-            unsafe fn svd(l: MatrixLayout, calc_u: bool, calc_vt: bool, mut a: &mut [Self]) -> Result<SVDOutput<Self>> {
+            unsafe fn svd(
+                l: MatrixLayout,
+                jobu: FlagSVD,
+                jobvt: FlagSVD,
+                mut a: &mut [Self],
+            ) -> Result<SVDOutput<Self>> {
                 let (m, n) = l.size();
                 let k = ::std::cmp::min(n, m);
                 let lda = l.lda();
-                let (ju, ldu, mut u) = if calc_u {
-                    (FlagSVD::All, m, vec![Self::zero(); (m * m) as usize])
-                } else {
-                    (FlagSVD::No, 1, Vec::new())
+                let ucol = match jobu {
+                    FlagSVD::All => m,
+                    FlagSVD::Some => k,
+                    FlagSVD::None => 0,
                 };
-                let (jvt, ldvt, mut vt) = if calc_vt {
-                    (FlagSVD::All, n, vec![Self::zero(); (n * n) as usize])
-                } else {
-                    (FlagSVD::No, n, Vec::new())
+                let vtrow = match jobvt {
+                    FlagSVD::All => n,
+                    FlagSVD::Some => k,
+                    FlagSVD::None => 0,
                 };
+                let mut u = vec![Self::zero(); (m * ucol).max(1) as usize];
+                let ldu = l.resized(m, ucol).lda();
+                let mut vt = vec![Self::zero(); (vtrow * n).max(1) as usize];
+                let ldvt = l.resized(vtrow, n).lda();
                 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,
-                    jvt as u8,
+                    jobu as u8,
+                    jobvt as u8,
                     m,
                     n,
                     &mut a,
@@ -71,16 +75,53 @@ macro_rules! impl_svd {
                     info,
                     SVDOutput {
                         s: s,
-                        u: if calc_u { Some(u) } else { None },
-                        vt: if calc_vt { Some(vt) } else { None },
+                        u: if jobu == FlagSVD::None { None } else { Some(u) },
+                        vt: if jobvt == FlagSVD::None { None } else { Some(vt) },
+                    },
+                )
+            }
+
+            unsafe fn svd_dc(l: MatrixLayout, jobz: FlagSVD, mut a: &mut [Self]) -> Result<SVDOutput<Self>> {
+                let (m, n) = l.size();
+                let k = m.min(n);
+                let lda = l.lda();
+                let (ucol, vtrow) = match jobz {
+                    FlagSVD::All => (m, n),
+                    FlagSVD::Some => (k, k),
+                    FlagSVD::None => (0, 0),
+                };
+                let mut s = vec![Self::Real::zero(); k.max(1) as usize];
+                let mut u = vec![Self::zero(); (m * ucol).max(1) as usize];
+                let ldu = l.resized(m, ucol).lda();
+                let mut vt = vec![Self::zero(); (vtrow * n).max(1) as usize];
+                let ldvt = l.resized(vtrow, n).lda();
+                let info = $gesdd(
+                    l.lapacke_layout(),
+                    jobz as u8,
+                    m,
+                    n,
+                    &mut a,
+                    lda,
+                    &mut s,
+                    &mut u,
+                    ldu,
+                    &mut vt,
+                    ldvt,
+                );
+                into_result(
+                    info,
+                    SVDOutput {
+                        s: s,
+                        u: if jobz == FlagSVD::None { None } else { Some(u) },
+                        vt: if jobz == FlagSVD::None { None } else { Some(vt) },
                     },
                 )
             }
         }
     };
 } // impl_svd!
 
-impl_svd!(f64, lapacke::dgesvd);
-impl_svd!(f32, lapacke::sgesvd);
-impl_svd!(c64, lapacke::zgesvd);
-impl_svd!(c32, lapacke::cgesvd);
+impl_svd!(f64, lapacke::dgesvd, lapacke::dgesdd);
+impl_svd!(f32, lapacke::sgesvd, lapacke::sgesdd);
+impl_svd!(c64, lapacke::zgesvd, lapacke::zgesdd);
+impl_svd!(c32, lapacke::cgesvd, lapacke::cgesdd);

src/lib.rs

@@ -57,7 +57,7 @@ pub mod qr;
 pub mod solve;
 pub mod solveh;
 pub mod svd;
-pub mod svddc;
+pub mod svd_rand;
 pub mod trace;
 pub mod triangular;
 pub mod types;
@@ -77,7 +77,7 @@ pub use qr::*;
 pub use solve::*;
 pub use solveh::*;
 pub use svd::*;
-pub use svddc::*;
+pub use svd_rand::*;
 pub use trace::*;
 pub use triangular::*;
 pub use types::*;