1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50

use crate::matrix::Matrix;
use crate::matrix::MatrixError;
use lapack::dgesvd;

impl Matrix {
  /// # Singular Value Decomposition
  ///
  /// https://en.wikipedia.org/wiki/Singular_value_decomposition
  ///
  /// `M = U * Sigma * V^T`
  /// `(u, sigma, vt)`
  pub fn gesvd(mut self) -> Result<(Matrix, Matrix, Matrix), MatrixError> {
    if self.rows != self.cols {
      return Err(MatrixError::DimensionMismatch);
    }

    let mut info = 0;
    let mut u = Matrix::new(self.rows, self.rows);
    let mut sigma = Matrix::new(self.rows, self.cols);
    let mut vt = Matrix::new(self.cols, self.cols);
    let lwork = 1usize.max(5usize * self.rows.min(self.cols));

    unsafe {
      dgesvd(
        'A' as u8,
        'A' as u8,
        self.rows as i32,
        self.cols as i32,
        &mut self.elems,
        self.rows as i32,
        &mut sigma.elems,
        &mut u.elems,
        self.rows as i32,
        &mut vt.elems,
        self.cols as i32,
        &mut vec![0.0; lwork],
        lwork as i32,
        &mut info,
      );
    }

    match info {
      0 => Ok((u, sigma, vt)),
      _ => Err(MatrixError::LapackRoutineError {
        routine: "dgesvd".to_owned(),
        info,
      }),
    }
  }
}