rusolver/
eigen.rs

1// SPDX-License-Identifier: Apache-2.0
2use crate::{
3    Matrix, MatrixView, SolverError, Tolerance
4};
5use crate::numerics::{
6    checked, norm_iter, symmetric, zeros
7};
8#[derive(Clone, Copy, Debug)]
9pub struct EigenOptions {
10    pub tolerance: Tolerance,
11    pub symmetry_tolerance: Tolerance,
12    pub max_sweeps: usize,
13}
14impl Default for EigenOptions {
15    fn default()->Self {
16        Self{
17            tolerance: Tolerance{
18                absolute: 0.0, relative: 1e-12
19            },
20            symmetry_tolerance: Tolerance::default(), max_sweeps: 64
21        }
22    }
23}
24#[derive(Clone, Debug)]
25pub struct SymmetricEigen {
26    /// Ascending eigenvalues. Eigenvectors are corresponding columns.
27    pub values: Vec<f64>, pub vectors: Matrix, pub sweeps: usize,
28    /// Off-diagonal Frobenius norm in original units; NOT a rigorous error bound.
29    pub off_diagonal_norm: f64,
30}
31/// Cyclic Jacobi rotations for a real symmetric, small/moderate dense matrix.
32/// Scales A by max(abs(A)); no general-complex eigensolver or cuSOLVER dispatch.
33/// A zero convergence tolerance and zero sweep budget are explicitly rejected.
34pub fn symmetric_eigen(a: MatrixView<'_>, options: EigenOptions)->Result<SymmetricEigen, SolverError> {
35    let n=a.square()?;
36    symmetric(a, options.symmetry_tolerance)?;
37    options.tolerance.validate()?;
38    if options.max_sweeps==0 || (options.tolerance.absolute==0.0 && options.tolerance.relative==0.0) {
39        return Err(SolverError::InvalidOption("eigen tolerance and sweep budget must be positive"));
40    }
41    let scale=a.max_abs();
42    let mut v=Matrix::identity(n)?;
43    if scale==0.0 {
44        return Ok(SymmetricEigen{
45            values: zeros(n)?, vectors: v, sweeps: 0, off_diagonal_norm: 0.0
46        });
47    }
48    let mut w=a.to_owned()?;
49    for x in &mut w.data {
50        *x/=scale;
51    }
52    // The lower triangle defines the accepted symmetric matrix.
53    for i in 0..n {
54        for j in 0..i {
55            w.data[j*n+i]=w.data[i*n+j];
56        }
57    }
58    let norm=norm_iter(w.data.iter().copied())?;
59    let cutoff=(options.tolerance.absolute/scale).max(options.tolerance.relative*norm);
60    for sweep in 0..=options.max_sweeps {
61        let off=norm_iter((0..n).flat_map(|i|(0..i).map(move|j|(i, j)))
62        .map(|(i, j)|w.data[i*n+j]))?*2.0f64.sqrt();
63        if off<=cutoff {
64            let mut order: Vec<usize>=(0..n).collect();
65            order.sort_by(|&i, &j|w.data[i*n+i].total_cmp(&w.data[j*n+j]));
66            let mut values=zeros(n)?;
67            let mut vectors=Matrix::zeros(n, n)?;
68            for (column, &old) in order.iter().enumerate() {
69                values[column]=checked(w.data[old*n+old]*scale, "eigenvalue rescale")?;
70                for i in 0..n {
71                    vectors.data[i*n+column]=v.data[i*n+old];
72                }
73            }
74            return Ok(SymmetricEigen{
75                values, vectors, sweeps: sweep,
76                off_diagonal_norm: checked(off*scale, "eigen residual rescale")?
77            });
78        }
79        if sweep==options.max_sweeps {
80            return Err(SolverError::NonConvergence{
81                iterations: sweep, residual: checked(off*scale, "eigen residual rescale")?
82            });
83        }
84        for p in 0..n {
85            for q in p+1..n {
86                let apq=w.data[p*n+q];
87                if apq==0.0 {
88                    continue;
89                }
90                let delta=(w.data[q*n+q]-w.data[p*n+p])*0.5;
91                let t=if delta==0.0 {
92                    1.0
93                } else {
94                    apq/(delta+delta.hypot(apq).copysign(delta))
95                };
96                let c=1.0/(1.0+t*t).sqrt();
97                let s=t*c;
98                let app=w.data[p*n+p];
99                let aqq=w.data[q*n+q];
100                w.data[p*n+p]=checked(app-t*apq, "Jacobi diagonal")?;
101                w.data[q*n+q]=checked(aqq+t*apq, "Jacobi diagonal")?;
102                w.data[p*n+q]=0.0;
103                w.data[q*n+p]=0.0;
104                for k in 0..n {
105                    if k!=p && k!=q {
106                        let x=w.data[k*n+p];
107                        let y=w.data[k*n+q];
108                        let xp=checked(c*x-s*y, "Jacobi rotation")?;
109                        let yq=checked(s*x+c*y, "Jacobi rotation")?;
110                        w.data[k*n+p]=xp;
111                        w.data[p*n+k]=xp;
112                        w.data[k*n+q]=yq;
113                        w.data[q*n+k]=yq;
114                    }
115                }
116                for k in 0..n {
117                    let x=v.data[k*n+p];
118                    let y=v.data[k*n+q];
119                    v.data[k*n+p]=c*x-s*y;
120                    v.data[k*n+q]=s*x+c*y;
121                }
122            }
123        }
124    }
125    Err(SolverError::Arithmetic("unreachable eigen iteration state"))
126}
rusolver/eigen.rs

rusolver/
eigen.rs
