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

ndarray_linalg/
eig.rs
