qualia_core_db/solvers/linear_algebra/
spectral.rs

1//! Matrix-spectral routines that bridge linear algebra and polynomial algebra:
2//! the characteristic polynomial and general (non-symmetric) eigenvalues.
3//!
4//! `characteristic_polynomial` uses Faddeev–LeVerrier; `eigenvalues_general` factors it
5//! with the engine's [`crate::solvers::polynomial::polynomial_roots`]. For symmetric
6//! matrices prefer [`super::eigen::symmetric_eigen`] (eigenvectors + better conditioning).
7
8use crate::solvers::polynomial::{polynomial_roots, Complex};
9use crate::solvers::SolversError;
10
11/// Characteristic polynomial of a row-major `n×n` matrix via the Faddeev–LeVerrier
12/// algorithm. Returns DESCENDING coefficients `[1, c₁, …, cₙ]` of
13/// `p(λ) = λⁿ + c₁λⁿ⁻¹ + … + cₙ` (so `det(A) = (-1)ⁿ·cₙ`). Exact for integer matrices;
14/// for large/ill-conditioned matrices prefer an iterative eigensolver.
15pub fn characteristic_polynomial(n: usize, data: &[f64]) -> Result<Vec<f64>, SolversError> {
16    if n == 0 || data.len() != n * n {
17        return Err(SolversError::InvalidDimension);
18    }
19    // Each Faddeev–LeVerrier step needs the dense product `A·M` (`n×n · n×n`). Route it
20    // through the engine's one GEMM (`super::gemm::matmul`), which itself picks the
21    // best path on this machine — offloading to the forge dispatcher above
22    // `GEMM_GPU_THRESHOLD` on an accelerator, and otherwise running its exact f64 CPU
23    // floor (same increasing-`k` accumulation, so byte-identical off-accelerator). The
24    // whole algorithm is `n` such products (O(n⁴)), so a large characteristic-polynomial
25    // is exactly the case the accelerator earns its keep.
26    let mul = |x: &[f64], y: &[f64]| -> Result<Vec<f64>, SolversError> {
27        let mut out = vec![0.0_f64; n * n];
28        super::gemm::matmul(n, n, n, x, y, &mut out)?;
29        Ok(out)
30    };
31    let trace = |x: &[f64]| -> f64 { (0..n).map(|i| x[i * n + i]).sum() };
32
33    let mut coeffs = vec![0.0_f64; n + 1];
34    coeffs[0] = 1.0;
35    // M starts as the identity.
36    let mut m = vec![0.0_f64; n * n];
37    for i in 0..n {
38        m[i * n + i] = 1.0;
39    }
40    for k in 1..=n {
41        let am = mul(data, &m)?;
42        let ck = -trace(&am) / (k as f64);
43        coeffs[k] = ck;
44        // M ← A·M + ck·I  (not needed after the final iteration).
45        m = am;
46        for i in 0..n {
47            m[i * n + i] += ck;
48        }
49    }
50    Ok(coeffs)
51}
52
53/// Eigenvalues of a GENERAL (not necessarily symmetric) row-major `n×n` matrix, as
54/// complex numbers. Computes the characteristic polynomial (Faddeev–LeVerrier) and finds
55/// its roots. Returns all `n` eigenvalues; real ones have `im ≈ 0`.
56pub fn eigenvalues_general(n: usize, data: &[f64]) -> Result<Vec<Complex>, SolversError> {
57    let charpoly = characteristic_polynomial(n, data)?;
58    polynomial_roots(&charpoly)
59}
60
61#[cfg(test)]
62mod tests {
63    use super::*;
64
65    #[test]
66    fn charpoly_of_2x2() {
67        // [[2,0],[0,3]] → (λ−2)(λ−3) = λ² − 5λ + 6 → [1, −5, 6]
68        let c = characteristic_polynomial(2, &[2.0, 0.0, 0.0, 3.0]).unwrap();
69        assert!((c[0] - 1.0).abs() < 1e-12);
70        assert!((c[1] + 5.0).abs() < 1e-9);
71        assert!((c[2] - 6.0).abs() < 1e-9);
72    }
73
74    #[test]
75    fn general_eigenvalues_real() {
76        // [[2,0],[0,3]] → eigenvalues {2,3}.
77        let ev = eigenvalues_general(2, &[2.0, 0.0, 0.0, 3.0]).unwrap();
78        assert_eq!(ev.len(), 2);
79        let mut reals: Vec<f64> = ev.iter().map(|z| z.re).collect();
80        reals.sort_by(|a, b| a.partial_cmp(b).unwrap());
81        assert!((reals[0] - 2.0).abs() < 1e-6 && (reals[1] - 3.0).abs() < 1e-6);
82    }
83
84    #[test]
85    fn general_eigenvalues_complex_rotation() {
86        // 90° rotation [[0,-1],[1,0]] → eigenvalues ±i.
87        let ev = eigenvalues_general(2, &[0.0, -1.0, 1.0, 0.0]).unwrap();
88        assert_eq!(ev.len(), 2);
89        assert!(ev.iter().any(|z| (z.im.abs() - 1.0).abs() < 1e-6));
90    }
91
92    #[test]
93    fn rejects_bad_dims() {
94        assert!(matches!(
95            characteristic_polynomial(2, &[1.0, 2.0, 3.0]),
96            Err(SolversError::InvalidDimension)
97        ));
98    }
99}
qualia_core_db/solvers/linear_algebra/spectral.rs

qualia_core_db/solvers/linear_algebra/
spectral.rs
