ndarray_linalg/
eig.rs

1//! Eigenvalue decomposition for non-symmetric square matrices
2
3use crate::error::*;
4use crate::layout::*;
5use crate::types::*;
6use ndarray::*;
7
8#[cfg_attr(doc, katexit::katexit)]
9/// Eigenvalue decomposition of general matrix reference
10pub trait Eig {
11    /// EigVec is the right eivenvector
12    type EigVal;
13    type EigVec;
14    /// Calculate eigenvalues with the right eigenvector
15    ///
16    /// $$ A u_i = \lambda_i u_i $$
17    ///
18    /// ```
19    /// use ndarray::*;
20    /// use ndarray_linalg::*;
21    ///
22    /// let a: Array2<f64> = array![
23    ///     [-1.01,  0.86, -4.60,  3.31, -4.81],
24    ///     [ 3.98,  0.53, -7.04,  5.29,  3.55],
25    ///     [ 3.30,  8.26, -3.89,  8.20, -1.51],
26    ///     [ 4.43,  4.96, -7.66, -7.33,  6.18],
27    ///     [ 7.31, -6.43, -6.16,  2.47,  5.58],
28    /// ];
29    /// let (eigs, vecs) = a.eig().unwrap();
30    ///
31    /// let a = a.map(|v| v.as_c());
32    /// for (&e, vec) in eigs.iter().zip(vecs.axis_iter(Axis(1))) {
33    ///     let ev = vec.map(|v| v * e);
34    ///     let av = a.dot(&vec);
35    ///     assert_close_l2!(&av, &ev, 1e-5);
36    /// }
37    /// ```
38    fn eig(&self) -> Result<(Self::EigVal, Self::EigVec)>;
39}
40
41impl<A, S> Eig for ArrayBase<S, Ix2>
42where
43    A: Scalar + Lapack,
44    S: Data<Elem = A>,
45{
46    type EigVal = Array1<A::Complex>;
47    type EigVec = Array2<A::Complex>;
48
49    fn eig(&self) -> Result<(Self::EigVal, Self::EigVec)> {
50        let mut a = self.to_owned();
51        let layout = a.square_layout()?;
52        let (s, t) = A::eig(true, layout, a.as_allocated_mut()?)?;
53        let n = layout.len() as usize;
54        Ok((
55            ArrayBase::from(s),
56            Array2::from_shape_vec((n, n).f(), t).unwrap(),
57        ))
58    }
59}
60
61/// Calculate eigenvalues without eigenvectors
62pub trait EigVals {
63    type EigVal;
64    fn eigvals(&self) -> Result<Self::EigVal>;
65}
66
67impl<A, S> EigVals for ArrayBase<S, Ix2>
68where
69    A: Scalar + Lapack,
70    S: Data<Elem = A>,
71{
72    type EigVal = Array1<A::Complex>;
73
74    fn eigvals(&self) -> Result<Self::EigVal> {
75        let mut a = self.to_owned();
76        let (s, _) = A::eig(false, a.square_layout()?, a.as_allocated_mut()?)?;
77        Ok(ArrayBase::from(s))
78    }
79}
ndarray_linalg/eig.rs

ndarray_linalg/
eig.rs
