M-Labs
/
nalgebra
3
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

first version

This commit is contained in:
Marc Haubenstock 2022-10-16 11:52:32 +02:00
parent 8ee68afaac
commit ee3f84abba
1 changed files with 50 additions and 48 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

98
nalgebra-lapack/src/eigen.rs
View File

@ -31,7 +31,7 @@ use lapack;
) )
)] )]
#[derive(Clone, Debug)] #[derive(Clone, Debug)]
pub struct Eigen<T: Scalar, D: Dim> pub struct Eigen<T: Scalar, D: DimName>
where where
DefaultAllocator: Allocator<T, D> + Allocator<T, D, D>, DefaultAllocator: Allocator<T, D> + Allocator<T, D, D>,
{ {
@ -45,7 +45,7 @@ where
pub left_eigenvectors: Option<OMatrix<T, D, D>>, pub left_eigenvectors: Option<OMatrix<T, D, D>>,
} }
impl<T: Scalar + Copy, D: Dim> Copy for Eigen<T, D> impl<T: Scalar + Copy, D: DimName> Copy for Eigen<T, D>
where where
DefaultAllocator: Allocator<T, D> + Allocator<T, D, D>, DefaultAllocator: Allocator<T, D> + Allocator<T, D, D>,
OVector<T, D>: Copy, OVector<T, D>: Copy,
@ -53,7 +53,7 @@ where
{ {
} }
impl<T: EigenScalar + RealField, D: Dim> Eigen<T, D> impl<T: EigenScalar + RealField, D: DimName> Eigen<T, D>
where where
DefaultAllocator: Allocator<T, D, D> + Allocator<T, D>, DefaultAllocator: Allocator<T, D, D> + Allocator<T, D>,
{ {
@ -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. /// 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. /// 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> { 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> {
panic!("TODO"); match !self.eigenvalues_are_real() {
// match !self.eigenvalues_are_real() { true => (None, None, None),
// true => (None, None, None), false => {
// false => { let number_of_elements = self.eigenvalues_re.nrows();
// 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 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 eigenvalues = Vec::<Complex<T>>::with_capacity(2*number_of_complex_entries); let mut eigenvectors = match self.eigenvectors.is_some() {
// let mut eigenvectors = match self.eigenvectors.is_some() { true => Some(Vec::<OVector<Complex<T>, D>>::with_capacity(2*number_of_complex_entries)),
// true => Some(Vec::<OVector<Complex<T>, D>>::with_capacity(2*number_of_complex_entries)), false => None
// false => None };
// }; let mut left_eigenvectors = match self.left_eigenvectors.is_some() {
// let mut left_eigenvectors = match self.left_eigenvectors.is_some() { true => Some(Vec::<OVector<Complex<T>, D>>::with_capacity(2*number_of_complex_entries)),
// true => Some(Vec::<OVector<Complex<T>, D>>::with_capacity(2*number_of_complex_entries)), false => None
// false => None };
// };
// let eigenvectors_raw = self.eigenvectors; for mut c in 0..number_of_elements {
// let left_eigenvectors_raw = self.left_eigenvectors; 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::<T>::new(self.eigenvalues_re[c].clone(), self.eigenvalues_im[c].clone()));
eigenvalues.push(Complex::<T>::new(self.eigenvalues_re[c].clone(), -self.eigenvalues_im[c].clone()));
// for mut i in 0..number_of_elements { if eigenvectors.is_some() {
// if self.eigenvalues_im[i] != T::zero() { let mut vec = OVector::<Complex<T>, D>::zeros();
// //Complex conjugate pairs of eigenvalues appear consecutively with the eigenvalue having the positive imaginary part first. let mut vec_conj = OVector::<Complex<T>, D>::zeros();
// 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() { for r in 0..number_of_elements {
// let mut r1_vec = OVector::<Complex<T>, D>::zeros(number_of_elements); vec[r] = Complex::<T>::new((&self.eigenvectors.as_ref()).unwrap()[(r,c)].clone(),(&self.eigenvectors.as_ref()).unwrap()[(r,c+1)].clone());
// let mut r1_vec_conj = OVector::<Complex<T>, D>::zeros(number_of_elements); vec_conj[r] = Complex::<T>::new((&self.eigenvectors.as_ref()).unwrap()[(r,c)].clone(),-(&self.eigenvectors.as_ref()).unwrap()[(r,c+1)].clone());
}
// 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.as_mut().unwrap().push(vec);
// eigenvectors.unwrap().push(r1_vec_conj); eigenvectors.as_mut().unwrap().push(vec_conj);
// } }
if left_eigenvectors.is_some() {
let mut vec = OVector::<Complex<T>, D>::zeros();
let mut vec_conj = OVector::<Complex<T>, D>::zeros();
// if left_eigenvectors.is_some() { for r in 0..number_of_elements {
// //TODO: Do the same for left vec[r] = Complex::<T>::new((&self.left_eigenvectors.as_ref()).unwrap()[(r,c)].clone(),(&self.left_eigenvectors.as_ref()).unwrap()[(r,c+1)].clone());
// } vec_conj[r] = Complex::<T>::new((&self.left_eigenvectors.as_ref()).unwrap()[(r,c)].clone(),-(&self.left_eigenvectors.as_ref()).unwrap()[(r,c+1)].clone());
}
// i += 1; left_eigenvectors.as_mut().unwrap().push(vec);
// } left_eigenvectors.as_mut().unwrap().push(vec_conj);
}
// } //skip next entry
// (Some(eigenvalues), left_eigenvectors, eigenvectors) c += 1;
// } }
// } }
(Some(eigenvalues), left_eigenvectors, eigenvectors)
}
}
} }