opensrdk_linear_algebra/matrix/ge/
svd.rs

1use crate::matrix::ge::Matrix;
2use crate::matrix::MatrixError;
3use lapack::dgesvd;
4
5impl Matrix {
6    /// # Singular Value Decomposition
7    ///
8    /// https://en.wikipedia.org/wiki/Singular_value_decomposition
9    ///
10    /// `M = U * Sigma * V^T`
11    /// `(u, sigma, vt)`
12    pub fn gesvd(mut self) -> Result<(Matrix, Matrix, Matrix), MatrixError> {
13        if self.rows != self.cols {
14            return Err(MatrixError::DimensionMismatch);
15        }
16
17        let mut info = 0;
18        let mut u = Matrix::new(self.rows, self.rows);
19        let mut sigma = Matrix::new(self.rows, self.cols);
20        let mut vt = Matrix::new(self.cols, self.cols);
21        let lwork = 1usize.max(5usize * self.rows.min(self.cols));
22
23        unsafe {
24            dgesvd(
25                'A' as u8,
26                'A' as u8,
27                self.rows as i32,
28                self.cols as i32,
29                &mut self.elems,
30                self.rows as i32,
31                &mut sigma.elems,
32                &mut u.elems,
33                self.rows as i32,
34                &mut vt.elems,
35                self.cols as i32,
36                &mut vec![0.0; lwork],
37                lwork as i32,
38                &mut info,
39            );
40        }
41
42        match info {
43            0 => Ok((u, sigma, vt)),
44            _ => Err(MatrixError::LapackRoutineError {
45                routine: "dgesvd".to_owned(),
46                info,
47            }),
48        }
49    }
50}