From d61f9e29ba1e4dd6a70dc26c9606af2f81350db4 Mon Sep 17 00:00:00 2001
From: Marc Haubenstock <marc.dhaubenstock@gmail.com>
Date: Sat, 15 Oct 2022 16:49:13 +0200
Subject: [PATCH] working on issue 1106

---
 nalgebra-lapack/src/eigen.rs | 65 +++++++++++++++++++++++++++++++++---
 1 file changed, 61 insertions(+), 4 deletions(-)

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<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)
+            }
+        }
+
+
+    }
+
 }
 
 /*