diff --git a/nalgebra-lapack/src/eigen.rs b/nalgebra-lapack/src/eigen.rs index 68b34b56..d6dc90d0 100644 --- a/nalgebra-lapack/src/eigen.rs +++ b/nalgebra-lapack/src/eigen.rs @@ -7,9 +7,8 @@ use num_complex::Complex; use simba::scalar::RealField; use crate::ComplexHelper; -use na::allocator::Allocator; -use na::dimension::{Const, Dim}; -use na::{DefaultAllocator, Matrix, OMatrix, OVector, Scalar}; +use na::dimension::{Const, Dim, DimName}; +use na::{DefaultAllocator, Matrix, OMatrix, OVector, Scalar, allocator::Allocator}; use lapack; @@ -148,7 +147,7 @@ where eigenvalues_re: wr, eigenvalues_im: wi, left_eigenvectors: vl, - eigenvectors: vr, + eigenvectors: vr }) } @@ -168,6 +167,64 @@ where det } + + /// Returns a tuple of vectors. The elements of the tuple are the complex eigenvalues, complex left eigenvectors and complex right eigenvectors respectively. + /// The elements appear as conjugate pairs within each vector, with the positive of the pair always being first. + pub fn get_complex_elements(&self) -> (Option>>, Option, D>>>, Option, D>>>) where DefaultAllocator: Allocator, D> { + match !self.eigenvalues_are_real() { + true => (None, None, None), + false => { + let number_of_elements = self.eigenvalues_re.nrows(); + let number_of_complex_entries = self.eigenvalues_im.iter().fold(0, |acc, e| if !e.is_zero() {acc + 1} else {acc}); + let mut eigenvalues = Vec::>::with_capacity(2*number_of_complex_entries); + let mut eigenvectors = match self.eigenvectors.is_some() { + true => Some(Vec::, D>>::with_capacity(2*number_of_complex_entries)), + false => None + }; + let mut left_eigenvectors = match self.left_eigenvectors.is_some() { + true => Some(Vec::, D>>::with_capacity(2*number_of_complex_entries)), + false => None + }; + + let eigenvectors_raw = self.eigenvectors; + let left_eigenvectors_raw = self.left_eigenvectors; + + for mut i in 0..number_of_elements { + if self.eigenvalues_im[i] != T::zero() { + //Complex conjugate pairs of eigenvalues appear consecutively with the eigenvalue having the positive imaginary part first. + eigenvalues.push(Complex::::new(self.eigenvalues_re[i].clone(), self.eigenvalues_im[i].clone())); + eigenvalues.push(Complex::::new(self.eigenvalues_re[i].clone(), -self.eigenvalues_im[i].clone())); + + if eigenvectors.is_some() { + let mut r1_vec = OVector::, D>::zeros(number_of_elements); + let mut r1_vec_conj = OVector::, D>::zeros(number_of_elements); + + for j in 0..number_of_elements { + r1_vec[j] = Complex::::new(self.eigenvectors.unwrap()[(i,j)].clone(),self.eigenvectors.unwrap()[(i,j+1)].clone()); + r1_vec_conj[j] = Complex::::new(self.eigenvectors.unwrap()[(i,j)].clone(),-self.eigenvectors.unwrap()[(i,j+1)].clone()); + } + + eigenvectors.unwrap().push(r1_vec); + eigenvectors.unwrap().push(r1_vec_conj); + } + + + if left_eigenvectors.is_some() { + //TODO: Do the same for left + } + + + i += 1; + } + + } + (Some(eigenvalues), left_eigenvectors, eigenvectors) + } + } + + + } + } /*