solver/linalg/
eigen.rs

1extern crate lapack_src;
2use lapack::dstev;
3
4#[doc = include_str!("../../docs/math/eigen_symmetric_tridiagonal.md")]
5///
6/// # Arguments
7/// * `d_in` - The diagonal elements of the matrix.
8/// * `e_in` - The off-diagonal elements of the matrix. `e[i]` is the element at `(i, i+1)`.
9///
10/// # Returns
11/// A tuple containing:
12/// * A vector of eigenvalues.
13/// * A flattened vector of eigenvectors (column-major). The element at row `i` and column `j` is at index `i + j * n`.
14///
15/// # Examples
16///
17/// ```
18/// use solver::linalg::eigen::eigen_symmetric_tridiagonal;
19///
20/// // Matrix:
21/// // [ 2 -1 ]
22/// // [-1  2 ]
23/// // Eigenvalues should be 1 and 3.
24///
25/// let d = vec![2.0, 2.0];
26/// let e = vec![-1.0]; // Off-diagonal
27///
28/// let (evals, evecs) = eigen_symmetric_tridiagonal(d, e);
29///
30/// assert!((evals[0] - 1.0).abs() < 1e-6);
31/// assert!((evals[1] - 3.0).abs() < 1e-6);
32/// ```
33pub fn eigen_symmetric_tridiagonal(d_in: Vec<f64>, e_in: Vec<f64>) -> (Vec<f64>, Vec<f64>) {
34    let n = d_in.len();
35    let mut d = d_in;
36    // LAPACK expects e to be of length n-1, but our input might be n (with last 0) or n-1.
37    // We need to ensure it matches what dstev expects.
38    // dstev: E (input/output) REAL array, dimension (N-1)
39    let mut e = e_in;
40    let e_slice = if e.len() == n {
41        &mut e[0..n - 1]
42    } else {
43        &mut e[..]
44    };
45
46    let mut z = vec![0.0; n * n]; // Output eigenvectors (column-major by default in LAPACK)
47    let ldz = n as i32;
48    let mut work = vec![0.0; std::cmp::max(1, 2 * n - 2)];
49    let mut info = 0;
50
51    unsafe {
52        dstev(
53            b'V', n as i32, &mut d, e_slice, &mut z, ldz, &mut work, &mut info,
54        );
55    }
56
57    if info != 0 {
58        panic!("LAPACK dstev failed with info = {}", info);
59    }
60
61    (d, z)
62}
solver/linalg/eigen.rs

solver/linalg/
eigen.rs
