mdarray-linalg 0.1.2

Linear algebra operations for mdarray, with multiple exchangeable backends
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111

//! Eigenvalue, eigenvector, and Schur decomposition utilities for general and Hermitian matrices
use mdarray::{DSlice, DTensor, Dense, Layout};
use num_complex::{Complex, ComplexFloat};
use thiserror::Error;

/// Error types related to eigenvalue decomposition
#[derive(Debug, Error)]
pub enum EigError {
    #[error("Backend error code: {0}")]
    BackendError(i32),

    #[error("Backend failed to converge: {iterations} iterations exceeded")]
    BackendDidNotConverge { iterations: i32 },

    #[error("Matrix must be square for eigenvalue decomposition")]
    NotSquareMatrix,
}

/// Holds the results of an eigenvalue decomposition, including
/// eigenvalues (complex) and optionally left and right eigenvectors
pub struct EigDecomp<T: ComplexFloat> {
    pub eigenvalues: DTensor<Complex<T::Real>, 2>,
    pub left_eigenvectors: Option<DTensor<Complex<T::Real>, 2>>,
    pub right_eigenvectors: Option<DTensor<Complex<T::Real>, 2>>,
}

/// Result type for eigenvalue decomposition, returning either an
/// `EigDecomp` or an `EigError`
pub type EigResult<T> = Result<EigDecomp<T>, EigError>;

/// Error types related to Schur decomposition
#[derive(Debug, Error)]
pub enum SchurError {
    #[error("Backend error code: {0}")]
    BackendError(i32),

    #[error("Backend failed to converge: {iterations} iterations exceeded")]
    BackendDidNotConverge { iterations: i32 },

    #[error("Matrix must be square for Schur decomposition")]
    NotSquareMatrix,
}

/// Holds the results of a Schur decomposition: A = Z * T * Z^H
/// where Z is unitary and T is upper-triangular (complex) or quasi-upper triangular (real)
pub struct SchurDecomp<T: ComplexFloat> {
    /// Schur form T (upper-triangular for complex, quasi-upper triangular for real)
    pub t: DTensor<T, 2>,
    /// Unitary Schur transformation matrix Z
    pub z: DTensor<T, 2>,
}

/// Result type for Schur decomposition, returning either a
/// `SchurDecomp` or a `SchurError`
pub type SchurResult<T> = Result<SchurDecomp<T>, SchurError>;

/// Eigenvalue decomposition operations of general and Hermitian/symmetric matrices
pub trait Eig<T: ComplexFloat> {
    /// Compute eigenvalues and right eigenvectors with new allocated matrices
    /// The matrix `A` satisfies: `A * v = λ * v` where v are the right eigenvectors
    fn eig<L: Layout>(&self, a: &mut DSlice<T, 2, L>) -> EigResult<T>;

    /// Compute eigenvalues and both left/right eigenvectors with new allocated matrices
    /// The matrix A satisfies: `A * vr = λ * vr` and `vl^H * A = λ * vl^H`
    /// where `vr` are right eigenvectors and `vl` are left eigenvectors
    fn eig_full<L: Layout>(&self, a: &mut DSlice<T, 2, L>) -> EigResult<T>;

    /// Compute only eigenvalues with new allocated vectors
    fn eig_values<L: Layout>(&self, a: &mut DSlice<T, 2, L>) -> EigResult<T>;

    // Note of october 2025: this method has been
    // temporary removed as it was very hard to provide a coherent
    // interface for the user that deals with all the cases and all
    // the backends.

    // // Compute eigenvalues and right eigenvectors, overwriting
    // existing matrices
    //fn eig_overwrite<L: Layout, Lr: Layout, Li:
    // Layout, Lv: Layout>( &self, a: &mut DSlice<T, 2, L>,
    // eigenvalues: &mut DSlice<Complex<T>, 2, Dense>,
    // right_eigenvectors: &mut DSlice<Complex<T>, 2, Dense>, ) ->
    // Result<(), EigError>;

    /// Compute eigenvalues and eigenvectors of a Hermitian matrix (input should be complex)
    fn eigh<L: Layout>(&self, a: &mut DSlice<T, 2, L>) -> EigResult<T>;

    /// Compute eigenvalues and eigenvectors of a symmetric matrix (input should be real)
    fn eigs<L: Layout>(&self, a: &mut DSlice<T, 2, L>) -> EigResult<T>;

    /// Compute Schur decomposition with new allocated matrices
    fn schur<L: Layout>(&self, a: &mut DSlice<T, 2, L>) -> SchurResult<T>;

    /// Compute Schur decomposition overwriting existing matrices
    fn schur_overwrite<L: Layout>(
        &self,
        a: &mut DSlice<T, 2, L>,
        t: &mut DSlice<T, 2, Dense>,
        z: &mut DSlice<T, 2, Dense>,
    ) -> Result<(), SchurError>;

    /// Compute Schur (complex) decomposition with new allocated matrices
    fn schur_complex<L: Layout>(&self, a: &mut DSlice<T, 2, L>) -> SchurResult<T>;

    /// Compute Schur complex) decomposition overwriting existing matrices
    fn schur_complex_overwrite<L: Layout>(
        &self,
        a: &mut DSlice<T, 2, L>,
        t: &mut DSlice<T, 2, Dense>,
        z: &mut DSlice<T, 2, Dense>,
    ) -> Result<(), SchurError>;
}