array_matrix/matrix/
eig.rs

1use num_complex::{Complex};
2use num_traits::{Float};
3
4use crate::{MMul, Diag, QRHouseholder, Matrix};
5
6const ITERATIONS: usize = 1000;
7
8pub trait Eig: Matrix
9{
10    type Output;
11
12    /// Returns an approximation of the eigenvalues of a matrix
13    /// 
14    /// eig(A)
15    /// 
16    /// # Examples
17    /// 
18    /// ```rust
19    /// let a = [
20    ///     [1.0, 2.0],
21    ///     [3.0, 4.0]
22    /// ];
23    /// let eig_a = eig(a);
24    /// ```
25    fn eig(&self) -> Self::Output;
26}
27
28
29impl<F: Float, const L: usize, const H: usize> Eig for [[Complex<F>; L]; H]
30where
31    [[Complex<F>; L]; H]: QRHouseholder,
32    <[[Complex<F>; L]; H] as QRHouseholder>::OutputR: Diag
33        + MMul<<[[Complex<F>; L]; H] as QRHouseholder>::OutputQ, Output = [[Complex<F>; L]; H]>
34{
35    type Output = <<[[Complex<F>; L]; H] as QRHouseholder>::OutputR as Diag>::Output;
36    fn eig(&self) -> Self::Output
37    {
38        let mut a = self.clone();
39        for _ in 0..ITERATIONS
40        {
41            let (q, r) = a.qr_householder();
42            a = r.mul(q)
43        }
44        let r = a.qr_householder().1;
45        r.diag()
46    }
47}
48
49impl<const L: usize, const H: usize> Eig for [[f32; L]; H]
50where
51    Self: Matrix,
52    [[Complex<f32>; L]; H]: Eig
53{
54    type Output = <[[Complex<f32>; L]; H] as Eig>::Output;
55    fn eig(&self) -> Self::Output
56    {
57        self.map(|ar| ar.map(|arc| Complex::from(arc))).eig()
58    }
59}
60
61impl<const L: usize, const H: usize> Eig for [[f64; L]; H]
62where
63    Self: Matrix,
64    [[Complex<f64>; L]; H]: Eig
65{
66    type Output = <[[Complex<f64>; L]; H] as Eig>::Output;
67    fn eig(&self) -> Self::Output
68    {
69        self.map(|ar| ar.map(|arc| Complex::from(arc))).eig()
70    }
71}