opensrdk_linear_algebra/matrix/ci/
evd.rs

1use crate::matrix::ci::CirculantMatrix;
2use crate::{number::c64, Matrix};
3use rayon::prelude::*;
4use rustfft::FftPlanner;
5use std::f64::consts::PI;
6
7impl CirculantMatrix<f64> {
8    pub fn cievd(&self) -> (Matrix<c64>, Vec<c64>) {
9        let row_elems = self.row_elems();
10        let n = row_elems.len();
11
12        let mut fourier_matrix: Matrix<c64> = Matrix::<c64>::new(n, n);
13        let omega = c64::new(0.0, 2.0 * PI / (n as f64)).exp();
14
15        // for i in 0..n {
16        //     for j in 0..i {
17        //         fourier_matrix[i][j] = fourier_matrix[j][i];
18        //     }
19        //     for j in i..n {
20        //         fourier_matrix[i][j] = omega.powi((i * j) as i32);
21        //     }
22        // }
23
24        fourier_matrix
25            .elems_mut()
26            .par_iter_mut()
27            .enumerate()
28            .map(|(k, elem)| ((k / n, k % n), elem))
29            .for_each(|((i, j), elem)| {
30                *elem = omega.powi((i * j) as i32);
31            });
32
33        let mut buffer = row_elems
34            .par_iter()
35            .map(|&e| c64::new(e, 0.0))
36            .collect::<Vec<c64>>();
37
38        let mut planner = FftPlanner::new();
39        let fft = planner.plan_fft_forward(n);
40        fft.process(&mut buffer);
41
42        (fourier_matrix, buffer)
43    }
44}
45
46#[cfg(test)]
47mod tests {
48    use super::*;
49    #[test]
50    fn it_works() {
51        let a = CirculantMatrix::new(vec![1.0, 2.0, 3.0]);
52        let diagonalized = a.cievd();
53
54        assert_eq!(diagonalized.1[0].re, 6.0);
55    }
56}