mdarray_linalg/
eig.rs

1//! Eigenvalue, eigenvector, and Schur decomposition utilities for general and Hermitian matrices
2use mdarray::{DSlice, DTensor, Dense, Layout};
3use num_complex::{Complex, ComplexFloat};
4use thiserror::Error;
5
6/// Error types related to eigenvalue decomposition
7#[derive(Debug, Error)]
8pub enum EigError {
9    #[error("Backend error code: {0}")]
10    BackendError(i32),
11
12    #[error("Backend failed to converge: {iterations} iterations exceeded")]
13    BackendDidNotConverge { iterations: i32 },
14
15    #[error("Matrix must be square for eigenvalue decomposition")]
16    NotSquareMatrix,
17}
18
19/// Holds the results of an eigenvalue decomposition, including
20/// eigenvalues (complex) and optionally left and right eigenvectors
21pub struct EigDecomp<T: ComplexFloat> {
22    pub eigenvalues: DTensor<Complex<T::Real>, 2>,
23    pub left_eigenvectors: Option<DTensor<Complex<T::Real>, 2>>,
24    pub right_eigenvectors: Option<DTensor<Complex<T::Real>, 2>>,
25}
26
27/// Result type for eigenvalue decomposition, returning either an
28/// `EigDecomp` or an `EigError`
29pub type EigResult<T> = Result<EigDecomp<T>, EigError>;
30
31/// Error types related to Schur decomposition
32#[derive(Debug, Error)]
33pub enum SchurError {
34    #[error("Backend error code: {0}")]
35    BackendError(i32),
36
37    #[error("Backend failed to converge: {iterations} iterations exceeded")]
38    BackendDidNotConverge { iterations: i32 },
39
40    #[error("Matrix must be square for Schur decomposition")]
41    NotSquareMatrix,
42}
43
44/// Holds the results of a Schur decomposition: A = Z * T * Z^H
45/// where Z is unitary and T is upper-triangular (complex) or quasi-upper triangular (real)
46pub struct SchurDecomp<T: ComplexFloat> {
47    /// Schur form T (upper-triangular for complex, quasi-upper triangular for real)
48    pub t: DTensor<T, 2>,
49    /// Unitary Schur transformation matrix Z
50    pub z: DTensor<T, 2>,
51}
52
53/// Result type for Schur decomposition, returning either a
54/// `SchurDecomp` or a `SchurError`
55pub type SchurResult<T> = Result<SchurDecomp<T>, SchurError>;
56
57/// Eigenvalue decomposition operations of general and Hermitian/symmetric matrices
58pub trait Eig<T: ComplexFloat> {
59    /// Compute eigenvalues and right eigenvectors with new allocated matrices
60    /// The matrix `A` satisfies: `A * v = λ * v` where v are the right eigenvectors
61    fn eig<L: Layout>(&self, a: &mut DSlice<T, 2, L>) -> EigResult<T>;
62
63    /// Compute eigenvalues and both left/right eigenvectors with new allocated matrices
64    /// The matrix A satisfies: `A * vr = λ * vr` and `vl^H * A = λ * vl^H`
65    /// where `vr` are right eigenvectors and `vl` are left eigenvectors
66    fn eig_full<L: Layout>(&self, a: &mut DSlice<T, 2, L>) -> EigResult<T>;
67
68    /// Compute only eigenvalues with new allocated vectors
69    fn eig_values<L: Layout>(&self, a: &mut DSlice<T, 2, L>) -> EigResult<T>;
70
71    // Note of october 2025: this method has been
72    // temporary removed as it was very hard to provide a coherent
73    // interface for the user that deals with all the cases and all
74    // the backends.
75
76    // // Compute eigenvalues and right eigenvectors, overwriting
77    // existing matrices
78    //fn eig_overwrite<L: Layout, Lr: Layout, Li:
79    // Layout, Lv: Layout>( &self, a: &mut DSlice<T, 2, L>,
80    // eigenvalues: &mut DSlice<Complex<T>, 2, Dense>,
81    // right_eigenvectors: &mut DSlice<Complex<T>, 2, Dense>, ) ->
82    // Result<(), EigError>;
83
84    /// Compute eigenvalues and eigenvectors of a Hermitian matrix (input should be complex)
85    fn eigh<L: Layout>(&self, a: &mut DSlice<T, 2, L>) -> EigResult<T>;
86
87    /// Compute eigenvalues and eigenvectors of a symmetric matrix (input should be real)
88    fn eigs<L: Layout>(&self, a: &mut DSlice<T, 2, L>) -> EigResult<T>;
89
90    /// Compute Schur decomposition with new allocated matrices
91    fn schur<L: Layout>(&self, a: &mut DSlice<T, 2, L>) -> SchurResult<T>;
92
93    /// Compute Schur decomposition overwriting existing matrices
94    fn schur_overwrite<L: Layout>(
95        &self,
96        a: &mut DSlice<T, 2, L>,
97        t: &mut DSlice<T, 2, Dense>,
98        z: &mut DSlice<T, 2, Dense>,
99    ) -> Result<(), SchurError>;
100
101    /// Compute Schur (complex) decomposition with new allocated matrices
102    fn schur_complex<L: Layout>(&self, a: &mut DSlice<T, 2, L>) -> SchurResult<T>;
103
104    /// Compute Schur complex) decomposition overwriting existing matrices
105    fn schur_complex_overwrite<L: Layout>(
106        &self,
107        a: &mut DSlice<T, 2, L>,
108        t: &mut DSlice<T, 2, Dense>,
109        z: &mut DSlice<T, 2, Dense>,
110    ) -> Result<(), SchurError>;
111}