peroxide/numerical/
eigen.rs

1//! To find Eigenvalue & Eigenvector
2//!
3//! * Reference : Press, William H., and William T. Vetterling. *Numerical Recipes.* Cambridge: Cambridge Univ. Press, 2007.
4
5pub use self::EigenMethod::*;
6use crate::structure::matrix::Matrix;
7use crate::traits::matrix::MatrixTrait;
8use crate::util::non_macro::eye_shape;
9
10#[derive(Debug, Copy, Clone)]
11pub enum EigenMethod {
12    Jacobi,
13}
14
15#[derive(Debug, Clone)]
16pub struct Eigen {
17    pub eigenvalue: Vec<f64>,
18    pub eigenvector: Matrix,
19}
20
21impl Eigen {
22    pub fn extract(self) -> (Vec<f64>, Matrix) {
23        (self.eigenvalue, self.eigenvector)
24    }
25}
26
27pub fn eigen(m: &Matrix, em: EigenMethod) -> Eigen {
28    match em {
29        Jacobi => {
30            let mat = m.clone();
31            let mut j = jacobi(mat);
32            j.iter();
33            j.extract()
34        }
35    }
36}
37
38// =============================================================================
39// Jacobi Method
40// =============================================================================
41/// To do Jacobi method
42///
43/// * Reference : Press, William H., and William T. Vetterling. *Numerical Recipes.* Cambridge: Cambridge Univ. Press, 2007.
44#[derive(Debug)]
45pub struct JacobiTemp {
46    pub n: usize,
47    pub a: Matrix,
48    pub v: Matrix,
49    pub d: Vec<f64>,
50    pub n_rot: usize,
51}
52
53pub fn jacobi(m: Matrix) -> JacobiTemp {
54    let n = m.row;
55    let v = eye_shape(n, m.shape);
56    let d = m.diag();
57    let a = m;
58
59    JacobiTemp {
60        n,
61        a,
62        v,
63        d,
64        n_rot: 0,
65    }
66}
67
68impl JacobiTemp {
69    pub fn new(m: Matrix) -> Self {
70        jacobi(m)
71    }
72
73    pub fn extract(self) -> Eigen {
74        let v = self.v;
75        let d = self.d;
76        Eigen {
77            eigenvalue: d,
78            eigenvector: v,
79        }
80    }
81
82    /// Main Jacobi traits
83    ///
84    /// * Reference : Press, William H., and William T. Vetterling. *Numerical Recipes.* Cambridge: Cambridge Univ. Press, 2007.
85    pub fn iter(&mut self) {
86        let a = &mut self.a;
87        let n = self.n;
88        let v = &mut self.v;
89        let mut b = a.diag();
90        let d = &mut self.d;
91        let mut z = vec![0f64; n];
92        let mut h: f64;
93        let mut _n_rot = self.n_rot;
94
95        for i in 1..51 {
96            let mut sm = 0f64;
97            for ip in 0..n - 1 {
98                for iq in ip + 1..n {
99                    sm += a[(ip, iq)].abs();
100                }
101            }
102            if sm == 0f64 {
103                eigsrt(d, v);
104                return ();
105            }
106            let tresh = if i < 4 {
107                0.2 * sm / (n.pow(2) as f64)
108            } else {
109                0f64
110            };
111            for ip in 0..n - 1 {
112                for iq in ip + 1..n {
113                    let g = 100f64 * a[(ip, iq)].abs();
114                    if i > 4 && g <= f64::EPSILON * d[ip].abs() && g <= f64::EPSILON * d[iq].abs() {
115                        a[(ip, iq)] = 0f64;
116                    } else if a[(ip, iq)].abs() > tresh {
117                        h = d[iq] - d[ip];
118                        let t = if g <= f64::EPSILON * h.abs() {
119                            a[(ip, iq)] / h
120                        } else {
121                            let theta = 0.5 * h / a[(ip, iq)];
122                            let mut temp = 1f64 / (theta.abs() + (1f64 + theta.powi(2)).sqrt());
123                            if theta < 0f64 {
124                                temp = -temp;
125                            }
126                            temp
127                        };
128                        let c = 1f64 / (1f64 + t.powi(2)).sqrt();
129                        let s = t * c;
130                        let tau = s / (1f64 + c);
131                        h = t * a[(ip, iq)];
132                        z[ip] -= h;
133                        z[iq] += h;
134                        d[ip] -= h;
135                        d[iq] += h;
136                        a[(ip, iq)] = 0f64;
137                        for j in 0..ip {
138                            rot(a, s, tau, j, ip, j, iq);
139                        }
140                        for j in ip + 1..iq {
141                            rot(a, s, tau, ip, j, j, iq);
142                        }
143                        for j in iq + 1..n {
144                            rot(a, s, tau, ip, j, iq, j);
145                        }
146                        for j in 0..n {
147                            rot(v, s, tau, j, ip, j, iq);
148                        }
149                        _n_rot += 1;
150                    }
151                }
152            }
153            for ip in 0..n {
154                b[ip] += z[ip];
155                d[ip] = b[ip];
156                z[ip] = 0f64;
157            }
158        }
159        assert!(false, "Too many iterations in routine jacobi");
160    }
161}
162
163fn rot(a: &mut Matrix, s: f64, tau: f64, i: usize, j: usize, k: usize, l: usize) {
164    let g = a[(i, j)];
165    let h = a[(k, l)];
166    a[(i, j)] = g - s * (h + g * tau);
167    a[(k, l)] = h + s * (g - h * tau);
168}
169
170/// Given eigenvalue & eigenvector, sorts thod eigenvalues into descending order
171///
172/// * Reference : Press, William H., and William T. Vetterling. *Numerical Recipes.* Cambridge: Cambridge Univ. Press, 2007.
173fn eigsrt(d: &mut Vec<f64>, v: &mut Matrix) {
174    let mut k: usize;
175    let n = d.len();
176    for i in 0..n - 1 {
177        k = i;
178        let mut p = d[k];
179        for j in i..n {
180            if d[j] >= p {
181                k = j;
182                p = d[k];
183            }
184        }
185        if k != i {
186            d[k] = d[i];
187            d[i] = p;
188            for j in 0..n {
189                p = v[(j, i)];
190                v[(j, i)] = v[(j, k)];
191                v[(j, k)] = p;
192            }
193        }
194    }
195}
peroxide/numerical/eigen.rs

peroxide/numerical/
eigen.rs
