ndarray_linalg/
svddc.rs

1//! Singular-value decomposition (SVD) by divide-and-conquer (?gesdd)
2
3use super::{convert::*, error::*, layout::*, types::*};
4use ndarray::*;
5
6pub use lax::JobSvd;
7
8/// Singular-value decomposition of matrix (copying) by divide-and-conquer
9pub trait SVDDC {
10    type U;
11    type VT;
12    type Sigma;
13    fn svddc(&self, uvt_flag: JobSvd) -> Result<(Option<Self::U>, Self::Sigma, Option<Self::VT>)>;
14}
15
16/// Singular-value decomposition of matrix by divide-and-conquer
17pub trait SVDDCInto {
18    type U;
19    type VT;
20    type Sigma;
21    fn svddc_into(
22        self,
23        uvt_flag: JobSvd,
24    ) -> Result<(Option<Self::U>, Self::Sigma, Option<Self::VT>)>;
25}
26
27/// Singular-value decomposition of matrix reference by divide-and-conquer
28pub trait SVDDCInplace {
29    type U;
30    type VT;
31    type Sigma;
32    fn svddc_inplace(
33        &mut self,
34        uvt_flag: JobSvd,
35    ) -> Result<(Option<Self::U>, Self::Sigma, Option<Self::VT>)>;
36}
37
38impl<A> SVDDC for ArrayRef<A, Ix2>
39where
40    A: Scalar + Lapack,
41{
42    type U = Array2<A>;
43    type VT = Array2<A>;
44    type Sigma = Array1<A::Real>;
45
46    fn svddc(&self, uvt_flag: JobSvd) -> Result<(Option<Self::U>, Self::Sigma, Option<Self::VT>)> {
47        self.to_owned().svddc_into(uvt_flag)
48    }
49}
50
51impl<A, S> SVDDCInto for ArrayBase<S, Ix2>
52where
53    A: Scalar + Lapack,
54    S: DataMut<Elem = A>,
55{
56    type U = Array2<A>;
57    type VT = Array2<A>;
58    type Sigma = Array1<A::Real>;
59
60    fn svddc_into(
61        mut self,
62        uvt_flag: JobSvd,
63    ) -> Result<(Option<Self::U>, Self::Sigma, Option<Self::VT>)> {
64        self.svddc_inplace(uvt_flag)
65    }
66}
67
68impl<A> SVDDCInplace for ArrayRef<A, Ix2>
69where
70    A: Scalar + Lapack,
71{
72    type U = Array2<A>;
73    type VT = Array2<A>;
74    type Sigma = Array1<A::Real>;
75
76    fn svddc_inplace(
77        &mut self,
78        uvt_flag: JobSvd,
79    ) -> Result<(Option<Self::U>, Self::Sigma, Option<Self::VT>)> {
80        let l = self.layout()?;
81        let svd_res = A::svddc(l, uvt_flag, self.as_allocated_mut()?)?;
82        let (m, n) = l.size();
83        let k = m.min(n);
84
85        let (u_col, vt_row) = match uvt_flag {
86            JobSvd::All => (m, n),
87            JobSvd::Some => (k, k),
88            JobSvd::None => (0, 0),
89        };
90
91        let u = svd_res
92            .u
93            .map(|u| into_matrix(l.resized(m, u_col), u).unwrap());
94
95        let vt = svd_res
96            .vt
97            .map(|vt| into_matrix(l.resized(vt_row, n), vt).unwrap());
98
99        let s = ArrayBase::from(svd_res.s);
100        Ok((u, s, vt))
101    }
102}