working on issue 1106
@@ -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<Vec<Complex<T>>>, Option<Vec<OVector<Complex<T>, D>>>, Option<Vec<OVector<Complex<T>, D>>>) where DefaultAllocator: Allocator<Complex<T>, 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::<Complex<T>>::with_capacity(2*number_of_complex_entries);
+ let mut eigenvectors = match self.eigenvectors.is_some() {
+ true => Some(Vec::<OVector<Complex<T>, D>>::with_capacity(2*number_of_complex_entries)),
+ false => None
+ };
+ let mut left_eigenvectors = match self.left_eigenvectors.is_some() {
+ true => Some(Vec::<OVector<Complex<T>, 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::<T>::new(self.eigenvalues_re[i].clone(), self.eigenvalues_im[i].clone()));
+ eigenvalues.push(Complex::<T>::new(self.eigenvalues_re[i].clone(), -self.eigenvalues_im[i].clone()));
+
+ if eigenvectors.is_some() {
+ let mut r1_vec = OVector::<Complex<T>, D>::zeros(number_of_elements);
+ let mut r1_vec_conj = OVector::<Complex<T>, D>::zeros(number_of_elements);
+
+ for j in 0..number_of_elements {
+ r1_vec[j] = Complex::<T>::new(self.eigenvectors.unwrap()[(i,j)].clone(),self.eigenvectors.unwrap()[(i,j+1)].clone());
+ r1_vec_conj[j] = Complex::<T>::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)
+ }
+ }
+
+
+ }
+
}
/*
