pub fn iter(&mut self) {
        let a = &mut self.a;
        let n = self.n;
        let v = &mut self.v;
        let mut b = a.diag();
        let d = &mut self.d;
        let mut z = vec![0f64; n];
        let mut h: f64;
        let mut _n_rot = self.n_rot;

        for i in 1..51 {
            let mut sm = 0f64;
            for ip in 0..n - 1 {
                for iq in ip + 1..n {
                    sm += a[(ip, iq)].abs();
                }
            }
            if sm == 0f64 {
                eigsrt(d, v);
                return ();
            }
            let tresh = if i < 4 {
                0.2 * sm / (n.pow(2) as f64)
            } else {
                0f64
            };
            for ip in 0..n - 1 {
                for iq in ip + 1..n {
                    let g = 100f64 * a[(ip, iq)].abs();
                    if i > 4 && g <= f64::EPSILON * d[ip].abs() && g <= f64::EPSILON * d[iq].abs() {
                        a[(ip, iq)] = 0f64;
                    } else if a[(ip, iq)].abs() > tresh {
                        h = d[iq] - d[ip];
                        let t = if g <= f64::EPSILON * h.abs() {
                            a[(ip, iq)] / h
                        } else {
                            let theta = 0.5 * h / a[(ip, iq)];
                            let mut temp = 1f64 / (theta.abs() + (1f64 + theta.powi(2)).sqrt());
                            if theta < 0f64 {
                                temp = -temp;
                            }
                            temp
                        };
                        let c = 1f64 / (1f64 + t.powi(2)).sqrt();
                        let s = t * c;
                        let tau = s / (1f64 + c);
                        h = t * a[(ip, iq)];
                        z[ip] -= h;
                        z[iq] += h;
                        d[ip] -= h;
                        d[iq] += h;
                        a[(ip, iq)] = 0f64;
                        for j in 0..ip {
                            rot(a, s, tau, j, ip, j, iq);
                        }
                        for j in ip + 1..iq {
                            rot(a, s, tau, ip, j, j, iq);
                        }
                        for j in iq + 1..n {
                            rot(a, s, tau, ip, j, iq, j);
                        }
                        for j in 0..n {
                            rot(v, s, tau, j, ip, j, iq);
                        }
                        _n_rot += 1;
                    }
                }
            }
            for ip in 0..n {
                b[ip] += z[ip];
                d[ip] = b[ip];
                z[ip] = 0f64;
            }
        }
        assert!(false, "Too many iterations in routine jacobi");
    }
Main Jacobi traits

* Reference : Press, William H., and William T. Vetterling. *Numerical Recipes.* Cambridge: Cambridge Univ. Press, 2007.
src/numerical/eigen.rs