From ee3f84abba0fdb027401d0e747f7accf41f2af23 Mon Sep 17 00:00:00 2001 From: Marc Haubenstock Date: Sun, 16 Oct 2022 11:52:32 +0200 Subject: [PATCH] first version --- nalgebra-lapack/src/eigen.rs | 98 ++++++++++++++++++------------------ 1 file changed, 50 insertions(+), 48 deletions(-) diff --git a/nalgebra-lapack/src/eigen.rs b/nalgebra-lapack/src/eigen.rs index d81468ef..c49c6dd3 100644 --- a/nalgebra-lapack/src/eigen.rs +++ b/nalgebra-lapack/src/eigen.rs @@ -31,7 +31,7 @@ use lapack; ) )] #[derive(Clone, Debug)] -pub struct Eigen +pub struct Eigen where DefaultAllocator: Allocator + Allocator, { @@ -45,7 +45,7 @@ where pub left_eigenvectors: Option>, } -impl Copy for Eigen +impl Copy for Eigen where DefaultAllocator: Allocator + Allocator, OVector: Copy, @@ -53,7 +53,7 @@ where { } -impl Eigen +impl Eigen where DefaultAllocator: Allocator + Allocator, { @@ -171,57 +171,59 @@ where /// 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> { - panic!("TODO"); - // 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 - // }; + 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 c in 0..number_of_elements { + if self.eigenvalues_im[c] != 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[c].clone(), self.eigenvalues_im[c].clone())); + eigenvalues.push(Complex::::new(self.eigenvalues_re[c].clone(), -self.eigenvalues_im[c].clone())); - // 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 vec = OVector::, D>::zeros(); + let mut vec_conj = OVector::, D>::zeros(); - // 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()); - // } + for r in 0..number_of_elements { + vec[r] = Complex::::new((&self.eigenvectors.as_ref()).unwrap()[(r,c)].clone(),(&self.eigenvectors.as_ref()).unwrap()[(r,c+1)].clone()); + vec_conj[r] = Complex::::new((&self.eigenvectors.as_ref()).unwrap()[(r,c)].clone(),-(&self.eigenvectors.as_ref()).unwrap()[(r,c+1)].clone()); + } - // eigenvectors.unwrap().push(r1_vec); - // eigenvectors.unwrap().push(r1_vec_conj); - // } + eigenvectors.as_mut().unwrap().push(vec); + eigenvectors.as_mut().unwrap().push(vec_conj); + } + if left_eigenvectors.is_some() { + let mut vec = OVector::, D>::zeros(); + let mut vec_conj = OVector::, D>::zeros(); - // if left_eigenvectors.is_some() { - // //TODO: Do the same for left - // } - - - // i += 1; - // } - - // } - // (Some(eigenvalues), left_eigenvectors, eigenvectors) - // } - // } + for r in 0..number_of_elements { + vec[r] = Complex::::new((&self.left_eigenvectors.as_ref()).unwrap()[(r,c)].clone(),(&self.left_eigenvectors.as_ref()).unwrap()[(r,c+1)].clone()); + vec_conj[r] = Complex::::new((&self.left_eigenvectors.as_ref()).unwrap()[(r,c)].clone(),-(&self.left_eigenvectors.as_ref()).unwrap()[(r,c+1)].clone()); + } + + left_eigenvectors.as_mut().unwrap().push(vec); + left_eigenvectors.as_mut().unwrap().push(vec_conj); + } + //skip next entry + c += 1; + } + } + (Some(eigenvalues), left_eigenvectors, eigenvectors) + } + } }