From 9f8cfa207af0b71d141d2eba3f9ce7ecb719d8e3 Mon Sep 17 00:00:00 2001 From: Marc Haubenstock Date: Mon, 2 May 2022 21:55:51 +0200 Subject: [PATCH 01/20] enabled complex eigenvalues for lapack --- nalgebra-lapack/src/eigen.rs | 157 ++++++++++++++++++++++++++++------- 1 file changed, 129 insertions(+), 28 deletions(-) diff --git a/nalgebra-lapack/src/eigen.rs b/nalgebra-lapack/src/eigen.rs index 8eab62d84..594cb0806 100644 --- a/nalgebra-lapack/src/eigen.rs +++ b/nalgebra-lapack/src/eigen.rs @@ -69,8 +69,8 @@ where "Unable to compute the eigenvalue decomposition of a non-square matrix." ); - let ljob = if left_eigenvectors { b'V' } else { b'T' }; - let rjob = if eigenvectors { b'V' } else { b'T' }; + let ljob = if left_eigenvectors { b'V' } else { b'N' }; + let rjob = if eigenvectors { b'V' } else { b'N' }; let (nrows, ncols) = m.shape_generic(); let n = nrows.value(); @@ -232,22 +232,27 @@ where /// The complex eigenvalues of the given matrix. /// /// Panics if the eigenvalue computation does not converge. - pub fn complex_eigenvalues(mut m: OMatrix) -> OVector, D> + pub fn complex_eigenvalues(mut m: OMatrix, left_eigenvectors: bool, eigenvectors: bool) + -> (OVector, D>, Option>, Option>) where - DefaultAllocator: Allocator, D>, + DefaultAllocator: Allocator, D> + Allocator, D, D>, { assert!( m.is_square(), "Unable to compute the eigenvalue decomposition of a non-square matrix." ); - let nrows = m.shape_generic().0; + let ljob = if left_eigenvectors { b'V' } else { b'N' }; + let rjob = if eigenvectors { b'V' } else { b'N' }; + + let (nrows, ncols) = m.shape_generic(); let n = nrows.value(); let lda = n as i32; // TODO: avoid the initialization? let mut wr = Matrix::zeros_generic(nrows, Const::<1>); + // TODO: Tap into the workspace. let mut wi = Matrix::zeros_generic(nrows, Const::<1>); let mut info = 0; @@ -255,8 +260,8 @@ where let mut placeholder2 = [T::zero()]; let lwork = T::xgeev_work_size( - b'T', - b'T', + ljob, + rjob, n as i32, m.as_mut_slice(), lda, @@ -273,31 +278,127 @@ where let mut work = vec![T::zero(); lwork as usize]; - T::xgeev( - b'T', - b'T', - n as i32, - m.as_mut_slice(), - lda, - wr.as_mut_slice(), - wi.as_mut_slice(), - &mut placeholder1, - 1 as i32, - &mut placeholder2, - 1 as i32, - &mut work, - lwork, - &mut info, - ); - lapack_panic!(info); + match (left_eigenvectors, eigenvectors) { + (true, true) => { + // TODO: avoid the initializations? + let mut vl = Matrix::zeros_generic(nrows, ncols); + let mut vr = Matrix::zeros_generic(nrows, ncols); + + T::xgeev( + ljob, + rjob, + n as i32, + m.as_mut_slice(), + lda, + wr.as_mut_slice(), + wi.as_mut_slice(), + &mut vl.as_mut_slice(), + n as i32, + &mut vr.as_mut_slice(), + n as i32, + &mut work, + lwork, + &mut info, + ); + lapack_panic!(info); + + let mut res = Matrix::zeros_generic(nrows, Const::<1>); + + for i in 0..res.len() { + res[i] = Complex::new(wr[i].clone(), wi[i].clone()); + } + + return (res, Some(vl), Some(vr)) + } + (true, false) => { + // TODO: avoid the initialization? + let mut vl = Matrix::zeros_generic(nrows, ncols); + + T::xgeev( + ljob, + rjob, + n as i32, + m.as_mut_slice(), + lda, + wr.as_mut_slice(), + wi.as_mut_slice(), + &mut vl.as_mut_slice(), + n as i32, + &mut placeholder2, + 1 as i32, + &mut work, + lwork, + &mut info, + ); + lapack_panic!(info); + + let mut res = Matrix::zeros_generic(nrows, Const::<1>); + + for i in 0..res.len() { + res[i] = Complex::new(wr[i].clone(), wi[i].clone()); + } + + return (res, Some(vl), None) + } + (false, true) => { + // TODO: avoid the initialization? + let mut vr = Matrix::zeros_generic(nrows, ncols); + + T::xgeev( + ljob, + rjob, + n as i32, + m.as_mut_slice(), + lda, + wr.as_mut_slice(), + wi.as_mut_slice(), + &mut placeholder1, + 1 as i32, + &mut vr.as_mut_slice(), + n as i32, + &mut work, + lwork, + &mut info, + ); + lapack_panic!(info); - let mut res = Matrix::zeros_generic(nrows, Const::<1>); + let mut res = Matrix::zeros_generic(nrows, Const::<1>); - for i in 0..res.len() { - res[i] = Complex::new(wr[i].clone(), wi[i].clone()); + for i in 0..res.len() { + res[i] = Complex::new(wr[i].clone(), wi[i].clone()); + } + + return (res, None, Some(vr)) + } + (false, false) => { + T::xgeev( + ljob, + rjob, + n as i32, + m.as_mut_slice(), + lda, + wr.as_mut_slice(), + wi.as_mut_slice(), + &mut placeholder1, + 1 as i32, + &mut placeholder2, + 1 as i32, + &mut work, + lwork, + &mut info, + ); + lapack_panic!(info); + + let mut res = Matrix::zeros_generic(nrows, Const::<1>); + + for i in 0..res.len() { + res[i] = Complex::new(wr[i].clone(), wi[i].clone()); + } + + return (res, None, None) + } } - res } /// The determinant of the decomposed matrix. From 71ef287fead3449a9914020b016b0207dc68d039 Mon Sep 17 00:00:00 2001 From: Marc Haubenstock Date: Mon, 2 May 2022 22:49:07 +0200 Subject: [PATCH 02/20] removed return statements --- nalgebra-lapack/src/eigen.rs | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/nalgebra-lapack/src/eigen.rs b/nalgebra-lapack/src/eigen.rs index 594cb0806..6c63ee336 100644 --- a/nalgebra-lapack/src/eigen.rs +++ b/nalgebra-lapack/src/eigen.rs @@ -308,7 +308,7 @@ where res[i] = Complex::new(wr[i].clone(), wi[i].clone()); } - return (res, Some(vl), Some(vr)) + (res, Some(vl), Some(vr)) } (true, false) => { // TODO: avoid the initialization? @@ -338,7 +338,7 @@ where res[i] = Complex::new(wr[i].clone(), wi[i].clone()); } - return (res, Some(vl), None) + (res, Some(vl), None) } (false, true) => { // TODO: avoid the initialization? @@ -368,7 +368,7 @@ where res[i] = Complex::new(wr[i].clone(), wi[i].clone()); } - return (res, None, Some(vr)) + (res, None, Some(vr)) } (false, false) => { T::xgeev( @@ -395,7 +395,7 @@ where res[i] = Complex::new(wr[i].clone(), wi[i].clone()); } - return (res, None, None) + (res, None, None) } } From aafb7138489e9cb48072e9d7ecc0f52365f908fa Mon Sep 17 00:00:00 2001 From: Marc Haubenstock Date: Tue, 3 May 2022 19:11:56 +0200 Subject: [PATCH 03/20] removed whitespace --- nalgebra-lapack/src/eigen.rs | 4 ---- 1 file changed, 4 deletions(-) diff --git a/nalgebra-lapack/src/eigen.rs b/nalgebra-lapack/src/eigen.rs index 6c63ee336..0c2106ad8 100644 --- a/nalgebra-lapack/src/eigen.rs +++ b/nalgebra-lapack/src/eigen.rs @@ -307,7 +307,6 @@ where for i in 0..res.len() { res[i] = Complex::new(wr[i].clone(), wi[i].clone()); } - (res, Some(vl), Some(vr)) } (true, false) => { @@ -337,7 +336,6 @@ where for i in 0..res.len() { res[i] = Complex::new(wr[i].clone(), wi[i].clone()); } - (res, Some(vl), None) } (false, true) => { @@ -367,7 +365,6 @@ where for i in 0..res.len() { res[i] = Complex::new(wr[i].clone(), wi[i].clone()); } - (res, None, Some(vr)) } (false, false) => { @@ -394,7 +391,6 @@ where for i in 0..res.len() { res[i] = Complex::new(wr[i].clone(), wi[i].clone()); } - (res, None, None) } } From 743eaa47474b63030a12a5c8cd7aaa7b5d8a890b Mon Sep 17 00:00:00 2001 From: Marc Haubenstock Date: Tue, 13 Sep 2022 19:09:27 +0200 Subject: [PATCH 04/20] renaming method --- nalgebra-lapack/src/eigen.rs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/nalgebra-lapack/src/eigen.rs b/nalgebra-lapack/src/eigen.rs index 0c2106ad8..1d78d036e 100644 --- a/nalgebra-lapack/src/eigen.rs +++ b/nalgebra-lapack/src/eigen.rs @@ -229,10 +229,10 @@ where None } - /// The complex eigenvalues of the given matrix. + /// The complex eigen decomposition of the given matrix. /// /// Panics if the eigenvalue computation does not converge. - pub fn complex_eigenvalues(mut m: OMatrix, left_eigenvectors: bool, eigenvectors: bool) + pub fn complex_eigen_decomposition(mut m: OMatrix, left_eigenvectors: bool, eigenvectors: bool) -> (OVector, D>, Option>, Option>) where DefaultAllocator: Allocator, D> + Allocator, D, D>, From c62801114b3c123ead9c9aba48a1905e9ed373d0 Mon Sep 17 00:00:00 2001 From: Marc Haubenstock Date: Wed, 14 Sep 2022 18:15:16 +0200 Subject: [PATCH 05/20] fixed formatting --- nalgebra-lapack/src/eigen.rs | 1 - 1 file changed, 1 deletion(-) diff --git a/nalgebra-lapack/src/eigen.rs b/nalgebra-lapack/src/eigen.rs index 1d78d036e..77c0f873a 100644 --- a/nalgebra-lapack/src/eigen.rs +++ b/nalgebra-lapack/src/eigen.rs @@ -394,7 +394,6 @@ where (res, None, None) } } - } /// The determinant of the decomposed matrix. From a0412c39f260adf146a52c2b44fc7d212f07be1d Mon Sep 17 00:00:00 2001 From: Marc Haubenstock Date: Wed, 14 Sep 2022 19:07:38 +0200 Subject: [PATCH 06/20] formatting --- nalgebra-lapack/src/eigen.rs | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/nalgebra-lapack/src/eigen.rs b/nalgebra-lapack/src/eigen.rs index 77c0f873a..bab9de83d 100644 --- a/nalgebra-lapack/src/eigen.rs +++ b/nalgebra-lapack/src/eigen.rs @@ -232,8 +232,15 @@ where /// The complex eigen decomposition of the given matrix. /// /// Panics if the eigenvalue computation does not converge. - pub fn complex_eigen_decomposition(mut m: OMatrix, left_eigenvectors: bool, eigenvectors: bool) - -> (OVector, D>, Option>, Option>) + pub fn complex_eigen_decomposition( + mut m: OMatrix, + left_eigenvectors: bool, + eigenvectors: bool, + ) -> ( + OVector, D>, + Option>, + Option>, + ) where DefaultAllocator: Allocator, D> + Allocator, D, D>, { From 3d52327f8298c1769a5b0a2e8e7a340a4d7cd46b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Crozet?= Date: Sun, 9 Oct 2022 16:26:26 +0200 Subject: [PATCH 07/20] nalgebra-lapack: merge both Eigen decomposition function into a single one. --- nalgebra-lapack/src/eigen.rs | 340 +++++------------------------------ 1 file changed, 47 insertions(+), 293 deletions(-) diff --git a/nalgebra-lapack/src/eigen.rs b/nalgebra-lapack/src/eigen.rs index bab9de83d..68b34b56d 100644 --- a/nalgebra-lapack/src/eigen.rs +++ b/nalgebra-lapack/src/eigen.rs @@ -13,7 +13,7 @@ use na::{DefaultAllocator, Matrix, OMatrix, OVector, Scalar}; use lapack; -/// Eigendecomposition of a real square matrix with real eigenvalues. +/// Eigendecomposition of a real square matrix with real or complex eigenvalues. #[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))] #[cfg_attr( feature = "serde-serialize", @@ -36,8 +36,10 @@ pub struct Eigen where DefaultAllocator: Allocator + Allocator, { - /// The eigenvalues of the decomposed matrix. - pub eigenvalues: OVector, + /// The real parts of eigenvalues of the decomposed matrix. + pub eigenvalues_re: OVector, + /// The imaginary parts of the eigenvalues of the decomposed matrix. + pub eigenvalues_im: OVector, /// The (right) eigenvectors of the decomposed matrix. pub eigenvectors: Option>, /// The left eigenvectors of the decomposed matrix. @@ -104,169 +106,27 @@ where lapack_check!(info); let mut work = vec![T::zero(); lwork as usize]; - - match (left_eigenvectors, eigenvectors) { - (true, true) => { - // TODO: avoid the initializations? - let mut vl = Matrix::zeros_generic(nrows, ncols); - let mut vr = Matrix::zeros_generic(nrows, ncols); - - T::xgeev( - ljob, - rjob, - n as i32, - m.as_mut_slice(), - lda, - wr.as_mut_slice(), - wi.as_mut_slice(), - &mut vl.as_mut_slice(), - n as i32, - &mut vr.as_mut_slice(), - n as i32, - &mut work, - lwork, - &mut info, - ); - lapack_check!(info); - - if wi.iter().all(|e| e.is_zero()) { - return Some(Self { - eigenvalues: wr, - left_eigenvectors: Some(vl), - eigenvectors: Some(vr), - }); - } - } - (true, false) => { - // TODO: avoid the initialization? - let mut vl = Matrix::zeros_generic(nrows, ncols); - - T::xgeev( - ljob, - rjob, - n as i32, - m.as_mut_slice(), - lda, - wr.as_mut_slice(), - wi.as_mut_slice(), - &mut vl.as_mut_slice(), - n as i32, - &mut placeholder2, - 1 as i32, - &mut work, - lwork, - &mut info, - ); - lapack_check!(info); - - if wi.iter().all(|e| e.is_zero()) { - return Some(Self { - eigenvalues: wr, - left_eigenvectors: Some(vl), - eigenvectors: None, - }); - } - } - (false, true) => { - // TODO: avoid the initialization? - let mut vr = Matrix::zeros_generic(nrows, ncols); - - T::xgeev( - ljob, - rjob, - n as i32, - m.as_mut_slice(), - lda, - wr.as_mut_slice(), - wi.as_mut_slice(), - &mut placeholder1, - 1 as i32, - &mut vr.as_mut_slice(), - n as i32, - &mut work, - lwork, - &mut info, - ); - lapack_check!(info); - - if wi.iter().all(|e| e.is_zero()) { - return Some(Self { - eigenvalues: wr, - left_eigenvectors: None, - eigenvectors: Some(vr), - }); - } - } - (false, false) => { - T::xgeev( - ljob, - rjob, - n as i32, - m.as_mut_slice(), - lda, - wr.as_mut_slice(), - wi.as_mut_slice(), - &mut placeholder1, - 1 as i32, - &mut placeholder2, - 1 as i32, - &mut work, - lwork, - &mut info, - ); - lapack_check!(info); - - if wi.iter().all(|e| e.is_zero()) { - return Some(Self { - eigenvalues: wr, - left_eigenvectors: None, - eigenvectors: None, - }); - } - } - } - - None - } - - /// The complex eigen decomposition of the given matrix. - /// - /// Panics if the eigenvalue computation does not converge. - pub fn complex_eigen_decomposition( - mut m: OMatrix, - left_eigenvectors: bool, - eigenvectors: bool, - ) -> ( - OVector, D>, - Option>, - Option>, - ) - where - DefaultAllocator: Allocator, D> + Allocator, D, D>, - { - assert!( - m.is_square(), - "Unable to compute the eigenvalue decomposition of a non-square matrix." - ); - - let ljob = if left_eigenvectors { b'V' } else { b'N' }; - let rjob = if eigenvectors { b'V' } else { b'N' }; - - let (nrows, ncols) = m.shape_generic(); - let n = nrows.value(); - - let lda = n as i32; - - // TODO: avoid the initialization? - let mut wr = Matrix::zeros_generic(nrows, Const::<1>); - // TODO: Tap into the workspace. - let mut wi = Matrix::zeros_generic(nrows, Const::<1>); - - let mut info = 0; - let mut placeholder1 = [T::zero()]; - let mut placeholder2 = [T::zero()]; - - let lwork = T::xgeev_work_size( + let mut vl = if left_eigenvectors { + Some(Matrix::zeros_generic(nrows, ncols)) + } else { + None + }; + let mut vr = if eigenvectors { + Some(Matrix::zeros_generic(nrows, ncols)) + } else { + None + }; + + let vl_ref = vl + .as_mut() + .map(|m| m.as_mut_slice()) + .unwrap_or(&mut placeholder1); + let vr_ref = vr + .as_mut() + .map(|m| m.as_mut_slice()) + .unwrap_or(&mut placeholder2); + + T::xgeev( ljob, rjob, n as i32, @@ -274,142 +134,36 @@ where lda, wr.as_mut_slice(), wi.as_mut_slice(), - &mut placeholder1, - n as i32, - &mut placeholder2, - n as i32, + vl_ref, + if left_eigenvectors { n as i32 } else { 1 }, + vr_ref, + if eigenvectors { n as i32 } else { 1 }, + &mut work, + lwork, &mut info, ); + lapack_check!(info); - lapack_panic!(info); - - let mut work = vec![T::zero(); lwork as usize]; - - match (left_eigenvectors, eigenvectors) { - (true, true) => { - // TODO: avoid the initializations? - let mut vl = Matrix::zeros_generic(nrows, ncols); - let mut vr = Matrix::zeros_generic(nrows, ncols); - - T::xgeev( - ljob, - rjob, - n as i32, - m.as_mut_slice(), - lda, - wr.as_mut_slice(), - wi.as_mut_slice(), - &mut vl.as_mut_slice(), - n as i32, - &mut vr.as_mut_slice(), - n as i32, - &mut work, - lwork, - &mut info, - ); - lapack_panic!(info); - - let mut res = Matrix::zeros_generic(nrows, Const::<1>); - - for i in 0..res.len() { - res[i] = Complex::new(wr[i].clone(), wi[i].clone()); - } - (res, Some(vl), Some(vr)) - } - (true, false) => { - // TODO: avoid the initialization? - let mut vl = Matrix::zeros_generic(nrows, ncols); - - T::xgeev( - ljob, - rjob, - n as i32, - m.as_mut_slice(), - lda, - wr.as_mut_slice(), - wi.as_mut_slice(), - &mut vl.as_mut_slice(), - n as i32, - &mut placeholder2, - 1 as i32, - &mut work, - lwork, - &mut info, - ); - lapack_panic!(info); - - let mut res = Matrix::zeros_generic(nrows, Const::<1>); - - for i in 0..res.len() { - res[i] = Complex::new(wr[i].clone(), wi[i].clone()); - } - (res, Some(vl), None) - } - (false, true) => { - // TODO: avoid the initialization? - let mut vr = Matrix::zeros_generic(nrows, ncols); - - T::xgeev( - ljob, - rjob, - n as i32, - m.as_mut_slice(), - lda, - wr.as_mut_slice(), - wi.as_mut_slice(), - &mut placeholder1, - 1 as i32, - &mut vr.as_mut_slice(), - n as i32, - &mut work, - lwork, - &mut info, - ); - lapack_panic!(info); - - let mut res = Matrix::zeros_generic(nrows, Const::<1>); - - for i in 0..res.len() { - res[i] = Complex::new(wr[i].clone(), wi[i].clone()); - } - (res, None, Some(vr)) - } - (false, false) => { - T::xgeev( - ljob, - rjob, - n as i32, - m.as_mut_slice(), - lda, - wr.as_mut_slice(), - wi.as_mut_slice(), - &mut placeholder1, - 1 as i32, - &mut placeholder2, - 1 as i32, - &mut work, - lwork, - &mut info, - ); - lapack_panic!(info); - - let mut res = Matrix::zeros_generic(nrows, Const::<1>); + Some(Self { + eigenvalues_re: wr, + eigenvalues_im: wi, + left_eigenvectors: vl, + eigenvectors: vr, + }) + } - for i in 0..res.len() { - res[i] = Complex::new(wr[i].clone(), wi[i].clone()); - } - (res, None, None) - } - } + /// Returns `true` if all the eigenvalues are real. + pub fn eigenvalues_are_real(&self) -> bool { + self.eigenvalues_im.iter().all(|e| e.is_zero()) } /// The determinant of the decomposed matrix. #[inline] #[must_use] - pub fn determinant(&self) -> T { - let mut det = T::one(); - for e in self.eigenvalues.iter() { - det *= e.clone(); + pub fn determinant(&self) -> Complex { + let mut det: Complex = na::one(); + for (re, im) in self.eigenvalues_re.iter().zip(self.eigenvalues_im.iter()) { + det *= Complex::new(re.clone(), im.clone()); } det From d61f9e29ba1e4dd6a70dc26c9606af2f81350db4 Mon Sep 17 00:00:00 2001 From: Marc Haubenstock Date: Sat, 15 Oct 2022 16:49:13 +0200 Subject: [PATCH 08/20] 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 68b34b56d..d6dc90d0e 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) + } + } + + + } + } /* From b732b75a774b2c51b97b9b0a1df192de8208b809 Mon Sep 17 00:00:00 2001 From: Marc Haubenstock Date: Sat, 15 Oct 2022 19:10:31 +0200 Subject: [PATCH 09/20] added panic so that the code compiles --- nalgebra-lapack/src/eigen.rs | 1 + 1 file changed, 1 insertion(+) diff --git a/nalgebra-lapack/src/eigen.rs b/nalgebra-lapack/src/eigen.rs index d6dc90d0e..0448679c6 100644 --- a/nalgebra-lapack/src/eigen.rs +++ b/nalgebra-lapack/src/eigen.rs @@ -171,6 +171,7 @@ 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 => { From 8ee68afaacfdbe5016727968a1a8e2212e0f4ed6 Mon Sep 17 00:00:00 2001 From: Marc Haubenstock Date: Sat, 15 Oct 2022 19:17:05 +0200 Subject: [PATCH 10/20] commented out so that code compiles --- nalgebra-lapack/src/eigen.rs | 88 ++++++++++++++++++------------------ 1 file changed, 44 insertions(+), 44 deletions(-) diff --git a/nalgebra-lapack/src/eigen.rs b/nalgebra-lapack/src/eigen.rs index 0448679c6..d81468ef5 100644 --- a/nalgebra-lapack/src/eigen.rs +++ b/nalgebra-lapack/src/eigen.rs @@ -172,56 +172,56 @@ where /// 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 - }; - - 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()); - } + // 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); - } + // eigenvectors.unwrap().push(r1_vec); + // eigenvectors.unwrap().push(r1_vec_conj); + // } - if left_eigenvectors.is_some() { - //TODO: Do the same for left - } + // if left_eigenvectors.is_some() { + // //TODO: Do the same for left + // } - i += 1; - } + // i += 1; + // } - } - (Some(eigenvalues), left_eigenvectors, eigenvectors) - } - } + // } + // (Some(eigenvalues), left_eigenvectors, eigenvectors) + // } + // } } From ee3f84abba0fdb027401d0e747f7accf41f2af23 Mon Sep 17 00:00:00 2001 From: Marc Haubenstock Date: Sun, 16 Oct 2022 11:52:32 +0200 Subject: [PATCH 11/20] first version --- nalgebra-lapack/src/eigen.rs | 108 ++++++++++++++++++----------------- 1 file changed, 55 insertions(+), 53 deletions(-) diff --git a/nalgebra-lapack/src/eigen.rs b/nalgebra-lapack/src/eigen.rs index d81468ef5..c49c6dd3f 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 - // }; - - // 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()); - // } + 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 + }; + + 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())); + + if eigenvectors.is_some() { + let mut vec = OVector::, D>::zeros(); + let mut vec_conj = OVector::, D>::zeros(); + + 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); - // } - - - // if left_eigenvectors.is_some() { - // //TODO: Do the same for left - // } - - - // i += 1; - // } - - // } - // (Some(eigenvalues), left_eigenvectors, eigenvectors) - // } - // } + 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(); + + 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) + } + } } From 01c5b9ea564d7b5be6c061423c659df9711b550f Mon Sep 17 00:00:00 2001 From: Marc Haubenstock Date: Sun, 16 Oct 2022 11:55:07 +0200 Subject: [PATCH 12/20] switch back to Dim + zero_generic --- nalgebra-lapack/src/eigen.rs | 25 +++++++++++++------------ 1 file changed, 13 insertions(+), 12 deletions(-) diff --git a/nalgebra-lapack/src/eigen.rs b/nalgebra-lapack/src/eigen.rs index c49c6dd3f..7048ee35a 100644 --- a/nalgebra-lapack/src/eigen.rs +++ b/nalgebra-lapack/src/eigen.rs @@ -7,7 +7,7 @@ use num_complex::Complex; use simba::scalar::RealField; use crate::ComplexHelper; -use na::dimension::{Const, Dim, DimName}; +use na::dimension::{Const, Dim}; use na::{DefaultAllocator, Matrix, OMatrix, OVector, Scalar, allocator::Allocator}; use lapack; @@ -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, { @@ -174,7 +174,8 @@ where match !self.eigenvalues_are_real() { true => (None, None, None), false => { - let number_of_elements = self.eigenvalues_re.nrows(); + let (number_of_elements, _) = self.eigenvalues_re.shape_generic(); + let number_of_elements_value = number_of_elements.value(); 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() { @@ -186,17 +187,17 @@ where false => None }; - for mut c in 0..number_of_elements { + for mut c in 0..number_of_elements_value { 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())); if eigenvectors.is_some() { - let mut vec = OVector::, D>::zeros(); - let mut vec_conj = OVector::, D>::zeros(); + let mut vec = OVector::, D>::zeros_generic(number_of_elements, Const::<1>); + let mut vec_conj = OVector::, D>::zeros_generic(number_of_elements, Const::<1>); - for r in 0..number_of_elements { + for r in 0..number_of_elements_value { 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()); } @@ -206,10 +207,10 @@ where } if left_eigenvectors.is_some() { - let mut vec = OVector::, D>::zeros(); - let mut vec_conj = OVector::, D>::zeros(); + let mut vec = OVector::, D>::zeros_generic(number_of_elements, Const::<1>); + let mut vec_conj = OVector::, D>::zeros_generic(number_of_elements, Const::<1>); - for r in 0..number_of_elements { + for r in 0..number_of_elements_value { 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()); } From 81f1a6d87e350314c5c293a4aab58fbc156f3b8a Mon Sep 17 00:00:00 2001 From: Marc Haubenstock Date: Sun, 16 Oct 2022 12:03:08 +0200 Subject: [PATCH 13/20] review conjugate indexing --- nalgebra-lapack/src/eigen.rs | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/nalgebra-lapack/src/eigen.rs b/nalgebra-lapack/src/eigen.rs index 7048ee35a..995d765e5 100644 --- a/nalgebra-lapack/src/eigen.rs +++ b/nalgebra-lapack/src/eigen.rs @@ -177,13 +177,13 @@ where let (number_of_elements, _) = self.eigenvalues_re.shape_generic(); let number_of_elements_value = number_of_elements.value(); 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 eigenvalues = Vec::>::with_capacity(number_of_complex_entries); let mut eigenvectors = match self.eigenvectors.is_some() { - true => Some(Vec::, D>>::with_capacity(2*number_of_complex_entries)), + true => Some(Vec::, D>>::with_capacity(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)), + true => Some(Vec::, D>>::with_capacity(number_of_complex_entries)), false => None }; @@ -191,7 +191,7 @@ where 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())); + eigenvalues.push(Complex::::new(self.eigenvalues_re[c+1].clone(), self.eigenvalues_im[c+1].clone())); if eigenvectors.is_some() { let mut vec = OVector::, D>::zeros_generic(number_of_elements, Const::<1>); @@ -199,7 +199,7 @@ where for r in 0..number_of_elements_value { 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()); + vec_conj[r] = Complex::::new((&self.eigenvectors.as_ref()).unwrap()[(r,c)].clone(),(&self.eigenvectors.as_ref()).unwrap()[(r,c+1)].clone()); } eigenvectors.as_mut().unwrap().push(vec); @@ -212,7 +212,7 @@ where for r in 0..number_of_elements_value { 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()); + 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); @@ -225,8 +225,6 @@ where (Some(eigenvalues), left_eigenvectors, eigenvectors) } } - - } } From 9acaa35e33964915507a6124fdae318af6503a7f Mon Sep 17 00:00:00 2001 From: Marc Haubenstock Date: Sat, 22 Oct 2022 21:26:09 +0200 Subject: [PATCH 14/20] fixed iteration --- nalgebra-lapack/src/eigen.rs | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/nalgebra-lapack/src/eigen.rs b/nalgebra-lapack/src/eigen.rs index 995d765e5..5f621dbfb 100644 --- a/nalgebra-lapack/src/eigen.rs +++ b/nalgebra-lapack/src/eigen.rs @@ -186,8 +186,9 @@ where true => Some(Vec::, D>>::with_capacity(number_of_complex_entries)), false => None }; - - for mut c in 0..number_of_elements_value { + + let mut c = 0; + while c < number_of_elements_value { 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())); From 4a1fd605e40afe7064cc4b3c306a75868466cdfb Mon Sep 17 00:00:00 2001 From: Marc Haubenstock Date: Sun, 23 Oct 2022 19:23:13 +0200 Subject: [PATCH 15/20] fixed iteration --- nalgebra-lapack/src/eigen.rs | 1 + 1 file changed, 1 insertion(+) diff --git a/nalgebra-lapack/src/eigen.rs b/nalgebra-lapack/src/eigen.rs index 5f621dbfb..3a848d0ab 100644 --- a/nalgebra-lapack/src/eigen.rs +++ b/nalgebra-lapack/src/eigen.rs @@ -222,6 +222,7 @@ where //skip next entry c += 1; } + c+=1; } (Some(eigenvalues), left_eigenvectors, eigenvectors) } From aa89cda0cd9474cf040adbfe55e4c26d75e6d9fd Mon Sep 17 00:00:00 2001 From: Marc Haubenstock Date: Sun, 23 Oct 2022 20:36:06 +0200 Subject: [PATCH 16/20] testing complex eigenvalues --- nalgebra-lapack/Cargo.toml | 4 +++- nalgebra-lapack/examples/complex_eigen.rs | 24 +++++++++++++++++++ nalgebra-lapack/src/eigen.rs | 3 ++- nalgebra-lapack/tests/linalg/complex_eigen.rs | 19 +++++++++++++++ nalgebra-lapack/tests/linalg/mod.rs | 1 + 5 files changed, 49 insertions(+), 2 deletions(-) create mode 100644 nalgebra-lapack/examples/complex_eigen.rs create mode 100644 nalgebra-lapack/tests/linalg/complex_eigen.rs diff --git a/nalgebra-lapack/Cargo.toml b/nalgebra-lapack/Cargo.toml index 91517a8d9..3f165f76d 100644 --- a/nalgebra-lapack/Cargo.toml +++ b/nalgebra-lapack/Cargo.toml @@ -22,7 +22,7 @@ proptest-support = [ "nalgebra/proptest-support" ] arbitrary = [ "nalgebra/arbitrary" ] # For BLAS/LAPACK -default = ["netlib"] +default = ["openblas"] openblas = ["lapack-src/openblas"] netlib = ["lapack-src/netlib"] accelerate = ["lapack-src/accelerate"] @@ -36,6 +36,7 @@ simba = "0.7" serde = { version = "1.0", features = [ "derive" ], optional = true } lapack = { version = "0.19", default-features = false } lapack-src = { version = "0.8", default-features = false } +openblas-src = {version = "*", features = ["static"]} # clippy = "*" [dev-dependencies] @@ -44,3 +45,4 @@ proptest = { version = "1", default-features = false, features = ["std"] } quickcheck = "1" approx = "0.5" rand = "0.8" + diff --git a/nalgebra-lapack/examples/complex_eigen.rs b/nalgebra-lapack/examples/complex_eigen.rs new file mode 100644 index 000000000..8a359d693 --- /dev/null +++ b/nalgebra-lapack/examples/complex_eigen.rs @@ -0,0 +1,24 @@ +extern crate nalgebra as na; +extern crate nalgebra_lapack; +#[macro_use] +extern crate approx; // for assert_relative_eq + +use na::Matrix3; +use nalgebra_lapack::Eigen; +use num_complex::Complex; + + +fn main() { + let m = Matrix3::::new(4.0/5.0, -3.0/5.0, 0.0, 3.0/5.0, 4.0/5.0, 0.0, 1.0, 2.0, 2.0); + let eigen = Eigen::new(m,true,true).expect("Eigen Creation Failed!"); + let (some_eigenvalues, some_left_vec, some_right_vec) = eigen.get_complex_elements(); + let eigenvalues = some_eigenvalues.expect("Eigenvalues Failed"); + let left_eigenvectors = some_left_vec.expect("Left Eigenvectors Failed"); + let eigenvectors = some_right_vec.expect("Right Eigenvectors Failed"); + + assert_relative_eq!(eigenvalues[0].re, Complex::::new(4.0/5.0,3.0/5.0).re); + assert_relative_eq!(eigenvalues[0].im, Complex::::new(4.0/5.0,3.0/5.0).im); + assert_relative_eq!(eigenvalues[1].re, Complex::::new(4.0/5.0,-3.0/5.0).re); + assert_relative_eq!(eigenvalues[1].im, Complex::::new(4.0/5.0,-3.0/5.0).im); + +} \ No newline at end of file diff --git a/nalgebra-lapack/src/eigen.rs b/nalgebra-lapack/src/eigen.rs index 3a848d0ab..d98c500d5 100644 --- a/nalgebra-lapack/src/eigen.rs +++ b/nalgebra-lapack/src/eigen.rs @@ -171,7 +171,8 @@ 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> { - match !self.eigenvalues_are_real() { + + match self.eigenvalues_are_real() { true => (None, None, None), false => { let (number_of_elements, _) = self.eigenvalues_re.shape_generic(); diff --git a/nalgebra-lapack/tests/linalg/complex_eigen.rs b/nalgebra-lapack/tests/linalg/complex_eigen.rs new file mode 100644 index 000000000..108694708 --- /dev/null +++ b/nalgebra-lapack/tests/linalg/complex_eigen.rs @@ -0,0 +1,19 @@ +use std::cmp; + +use na::{Matrix3}; +use nalgebra_lapack::Eigen; + +use crate::proptest::*; +use proptest::{prop_assert, proptest}; + +proptest! { + //#[test] + // fn complex_eigen() { + // let n = cmp::max(1, cmp::min(n, 10)); + // let m = DMatrix::::new_random(n, n); + // let eig = SymmetricEigen::new(m.clone()); + // let recomp = eig.recompose(); + // prop_assert!(relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5)) + // } + +} diff --git a/nalgebra-lapack/tests/linalg/mod.rs b/nalgebra-lapack/tests/linalg/mod.rs index 251bbe7b6..8fc8deebb 100644 --- a/nalgebra-lapack/tests/linalg/mod.rs +++ b/nalgebra-lapack/tests/linalg/mod.rs @@ -7,3 +7,4 @@ mod real_eigensystem; mod schur; mod svd; mod symmetric_eigen; +mod complex_eigen; From 84c37b79ddcfb97f961b76da0cee69d1a07f42b9 Mon Sep 17 00:00:00 2001 From: Marc Haubenstock Date: Sun, 23 Oct 2022 21:43:33 +0200 Subject: [PATCH 17/20] testing eigenvectors --- nalgebra-lapack/examples/complex_eigen.rs | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/nalgebra-lapack/examples/complex_eigen.rs b/nalgebra-lapack/examples/complex_eigen.rs index 8a359d693..ebf5e48c8 100644 --- a/nalgebra-lapack/examples/complex_eigen.rs +++ b/nalgebra-lapack/examples/complex_eigen.rs @@ -7,18 +7,27 @@ use na::Matrix3; use nalgebra_lapack::Eigen; use num_complex::Complex; - +//Matrix taken from https://textbooks.math.gatech.edu/ila/1553/complex-eigenvalues.html fn main() { let m = Matrix3::::new(4.0/5.0, -3.0/5.0, 0.0, 3.0/5.0, 4.0/5.0, 0.0, 1.0, 2.0, 2.0); let eigen = Eigen::new(m,true,true).expect("Eigen Creation Failed!"); let (some_eigenvalues, some_left_vec, some_right_vec) = eigen.get_complex_elements(); let eigenvalues = some_eigenvalues.expect("Eigenvalues Failed"); - let left_eigenvectors = some_left_vec.expect("Left Eigenvectors Failed"); + let _left_eigenvectors = some_left_vec.expect("Left Eigenvectors Failed"); let eigenvectors = some_right_vec.expect("Right Eigenvectors Failed"); assert_relative_eq!(eigenvalues[0].re, Complex::::new(4.0/5.0,3.0/5.0).re); assert_relative_eq!(eigenvalues[0].im, Complex::::new(4.0/5.0,3.0/5.0).im); assert_relative_eq!(eigenvalues[1].re, Complex::::new(4.0/5.0,-3.0/5.0).re); assert_relative_eq!(eigenvalues[1].im, Complex::::new(4.0/5.0,-3.0/5.0).im); + + + assert_relative_eq!(eigenvectors[0][0].re, -12.0/32.7871926215100059134410999); + assert_relative_eq!(eigenvectors[0][0].im, -9.0/32.7871926215100059134410999); + assert_relative_eq!(eigenvectors[0][1].re, -9.0/32.7871926215100059134410999); + assert_relative_eq!(eigenvectors[0][1].im, 12.0/32.7871926215100059134410999); + assert_relative_eq!(eigenvectors[0][2].re, 25.0/32.7871926215100059134410999); + assert_relative_eq!(eigenvectors[0][2].im, 0.0); + } \ No newline at end of file From 9c8b5f0f385a24d8be3e5b1625f511e6957498a8 Mon Sep 17 00:00:00 2001 From: Marc Haubenstock Date: Wed, 26 Oct 2022 16:26:06 +0200 Subject: [PATCH 18/20] added a function to get all the real elements --- nalgebra-lapack/examples/complex_eigen.rs | 4 +-- nalgebra-lapack/src/eigen.rs | 37 ++++++++++++++++++++++- 2 files changed, 37 insertions(+), 4 deletions(-) diff --git a/nalgebra-lapack/examples/complex_eigen.rs b/nalgebra-lapack/examples/complex_eigen.rs index ebf5e48c8..6c452d661 100644 --- a/nalgebra-lapack/examples/complex_eigen.rs +++ b/nalgebra-lapack/examples/complex_eigen.rs @@ -27,7 +27,5 @@ fn main() { assert_relative_eq!(eigenvectors[0][1].re, -9.0/32.7871926215100059134410999); assert_relative_eq!(eigenvectors[0][1].im, 12.0/32.7871926215100059134410999); assert_relative_eq!(eigenvectors[0][2].re, 25.0/32.7871926215100059134410999); - assert_relative_eq!(eigenvectors[0][2].im, 0.0); - - + assert_relative_eq!(eigenvectors[0][2].im, 0.0); } \ No newline at end of file diff --git a/nalgebra-lapack/src/eigen.rs b/nalgebra-lapack/src/eigen.rs index d98c500d5..e76b8ef42 100644 --- a/nalgebra-lapack/src/eigen.rs +++ b/nalgebra-lapack/src/eigen.rs @@ -168,10 +168,45 @@ where det } + /// Returns a tuple of vectors. The elements of the tuple are the real parts of the eigenvalues, left eigenvectors and right eigenvectors respectively. + pub fn get_real_elements(&self) -> (Vec, Option>>, Option>>) where DefaultAllocator: Allocator { + let (number_of_elements, _) = self.eigenvalues_re.shape_generic(); + let number_of_elements_value = number_of_elements.value(); + let mut eigenvalues = Vec::::with_capacity(number_of_elements_value); + let mut eigenvectors = match self.eigenvectors.is_some() { + true => Some(Vec::>::with_capacity(number_of_elements_value)), + false => None + }; + let mut left_eigenvectors = match self.left_eigenvectors.is_some() { + true => Some(Vec::>::with_capacity(number_of_elements_value)), + false => None + }; + + let mut c = 0; + while c < number_of_elements_value { + eigenvalues.push(self.eigenvalues_re[c].clone()); + + if eigenvectors.is_some() { + eigenvectors.as_mut().unwrap().push((&self.eigenvectors.as_ref()).unwrap().column(c).into_owned()); + } + + if left_eigenvectors.is_some() { + left_eigenvectors.as_mut().unwrap().push((&self.left_eigenvectors.as_ref()).unwrap().column(c).into_owned()); + } + if self.eigenvalues_im[c] != T::zero() { + //skip next entry + c += 1; + } + c+=1; + } + (eigenvalues, left_eigenvectors, eigenvectors) + + + } + /// 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 => { From e32f4ee16f402296e4da0fb2a824ac966c90c89d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Crozet?= Date: Sun, 30 Oct 2022 17:22:08 +0100 Subject: [PATCH 19/20] cargo fmt + tests --- nalgebra-lapack/examples/complex_eigen.rs | 31 ---- nalgebra-lapack/src/eigen.rs | 138 +++++++++++++----- nalgebra-lapack/tests/linalg/cholesky.rs | 6 +- nalgebra-lapack/tests/linalg/complex_eigen.rs | 56 +++++-- nalgebra-lapack/tests/linalg/lu.rs | 16 +- nalgebra-lapack/tests/linalg/mod.rs | 2 +- .../tests/linalg/real_eigensystem.rs | 42 +++--- nalgebra-lapack/tests/linalg/schur.rs | 13 +- 8 files changed, 188 insertions(+), 116 deletions(-) delete mode 100644 nalgebra-lapack/examples/complex_eigen.rs diff --git a/nalgebra-lapack/examples/complex_eigen.rs b/nalgebra-lapack/examples/complex_eigen.rs deleted file mode 100644 index 6c452d661..000000000 --- a/nalgebra-lapack/examples/complex_eigen.rs +++ /dev/null @@ -1,31 +0,0 @@ -extern crate nalgebra as na; -extern crate nalgebra_lapack; -#[macro_use] -extern crate approx; // for assert_relative_eq - -use na::Matrix3; -use nalgebra_lapack::Eigen; -use num_complex::Complex; - -//Matrix taken from https://textbooks.math.gatech.edu/ila/1553/complex-eigenvalues.html -fn main() { - let m = Matrix3::::new(4.0/5.0, -3.0/5.0, 0.0, 3.0/5.0, 4.0/5.0, 0.0, 1.0, 2.0, 2.0); - let eigen = Eigen::new(m,true,true).expect("Eigen Creation Failed!"); - let (some_eigenvalues, some_left_vec, some_right_vec) = eigen.get_complex_elements(); - let eigenvalues = some_eigenvalues.expect("Eigenvalues Failed"); - let _left_eigenvectors = some_left_vec.expect("Left Eigenvectors Failed"); - let eigenvectors = some_right_vec.expect("Right Eigenvectors Failed"); - - assert_relative_eq!(eigenvalues[0].re, Complex::::new(4.0/5.0,3.0/5.0).re); - assert_relative_eq!(eigenvalues[0].im, Complex::::new(4.0/5.0,3.0/5.0).im); - assert_relative_eq!(eigenvalues[1].re, Complex::::new(4.0/5.0,-3.0/5.0).re); - assert_relative_eq!(eigenvalues[1].im, Complex::::new(4.0/5.0,-3.0/5.0).im); - - - assert_relative_eq!(eigenvectors[0][0].re, -12.0/32.7871926215100059134410999); - assert_relative_eq!(eigenvectors[0][0].im, -9.0/32.7871926215100059134410999); - assert_relative_eq!(eigenvectors[0][1].re, -9.0/32.7871926215100059134410999); - assert_relative_eq!(eigenvectors[0][1].im, 12.0/32.7871926215100059134410999); - assert_relative_eq!(eigenvectors[0][2].re, 25.0/32.7871926215100059134410999); - assert_relative_eq!(eigenvectors[0][2].im, 0.0); -} \ No newline at end of file diff --git a/nalgebra-lapack/src/eigen.rs b/nalgebra-lapack/src/eigen.rs index e76b8ef42..08f161156 100644 --- a/nalgebra-lapack/src/eigen.rs +++ b/nalgebra-lapack/src/eigen.rs @@ -8,7 +8,7 @@ use simba::scalar::RealField; use crate::ComplexHelper; use na::dimension::{Const, Dim}; -use na::{DefaultAllocator, Matrix, OMatrix, OVector, Scalar, allocator::Allocator}; +use na::{allocator::Allocator, DefaultAllocator, Matrix, OMatrix, OVector, Scalar}; use lapack; @@ -147,7 +147,7 @@ where eigenvalues_re: wr, eigenvalues_im: wi, left_eigenvectors: vl, - eigenvectors: vr + eigenvectors: vr, }) } @@ -168,103 +168,169 @@ where det } - /// Returns a tuple of vectors. The elements of the tuple are the real parts of the eigenvalues, left eigenvectors and right eigenvectors respectively. - pub fn get_real_elements(&self) -> (Vec, Option>>, Option>>) where DefaultAllocator: Allocator { + /// Returns a tuple of vectors. The elements of the tuple are the real parts of the eigenvalues, left eigenvectors and right eigenvectors respectively. + pub fn get_real_elements( + &self, + ) -> ( + Vec, + Option>>, + Option>>, + ) + where + DefaultAllocator: Allocator, + { let (number_of_elements, _) = self.eigenvalues_re.shape_generic(); let number_of_elements_value = number_of_elements.value(); let mut eigenvalues = Vec::::with_capacity(number_of_elements_value); let mut eigenvectors = match self.eigenvectors.is_some() { - true => Some(Vec::>::with_capacity(number_of_elements_value)), - false => None + true => Some(Vec::>::with_capacity( + number_of_elements_value, + )), + false => None, }; let mut left_eigenvectors = match self.left_eigenvectors.is_some() { - true => Some(Vec::>::with_capacity(number_of_elements_value)), - false => None + true => Some(Vec::>::with_capacity( + number_of_elements_value, + )), + false => None, }; - + let mut c = 0; while c < number_of_elements_value { eigenvalues.push(self.eigenvalues_re[c].clone()); if eigenvectors.is_some() { - eigenvectors.as_mut().unwrap().push((&self.eigenvectors.as_ref()).unwrap().column(c).into_owned()); + eigenvectors.as_mut().unwrap().push( + (&self.eigenvectors.as_ref()) + .unwrap() + .column(c) + .into_owned(), + ); } if left_eigenvectors.is_some() { - left_eigenvectors.as_mut().unwrap().push((&self.left_eigenvectors.as_ref()).unwrap().column(c).into_owned()); + left_eigenvectors.as_mut().unwrap().push( + (&self.left_eigenvectors.as_ref()) + .unwrap() + .column(c) + .into_owned(), + ); } if self.eigenvalues_im[c] != T::zero() { //skip next entry c += 1; } - c+=1; + c += 1; } (eigenvalues, left_eigenvectors, eigenvectors) - - } - /// 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. - pub fn get_complex_elements(&self) -> (Option>>, Option, D>>>, Option, D>>>) where DefaultAllocator: Allocator, D> { + 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.shape_generic(); let number_of_elements_value = number_of_elements.value(); - 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::>::with_capacity(number_of_complex_entries); let mut eigenvectors = match self.eigenvectors.is_some() { - true => Some(Vec::, D>>::with_capacity(number_of_complex_entries)), - false => None + true => Some(Vec::, D>>::with_capacity( + number_of_complex_entries, + )), + false => None, }; let mut left_eigenvectors = match self.left_eigenvectors.is_some() { - true => Some(Vec::, D>>::with_capacity(number_of_complex_entries)), - false => None + true => Some(Vec::, D>>::with_capacity( + number_of_complex_entries, + )), + false => None, }; - + let mut c = 0; while c < number_of_elements_value { 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+1].clone(), self.eigenvalues_im[c+1].clone())); + eigenvalues.push(Complex::::new( + self.eigenvalues_re[c].clone(), + self.eigenvalues_im[c].clone(), + )); + eigenvalues.push(Complex::::new( + self.eigenvalues_re[c + 1].clone(), + self.eigenvalues_im[c + 1].clone(), + )); if eigenvectors.is_some() { - let mut vec = OVector::, D>::zeros_generic(number_of_elements, Const::<1>); - let mut vec_conj = OVector::, D>::zeros_generic(number_of_elements, Const::<1>); + let mut vec = OVector::, D>::zeros_generic( + number_of_elements, + Const::<1>, + ); + let mut vec_conj = OVector::, D>::zeros_generic( + number_of_elements, + Const::<1>, + ); for r in 0..number_of_elements_value { - 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()); + 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.as_mut().unwrap().push(vec); eigenvectors.as_mut().unwrap().push(vec_conj); } if left_eigenvectors.is_some() { - let mut vec = OVector::, D>::zeros_generic(number_of_elements, Const::<1>); - let mut vec_conj = OVector::, D>::zeros_generic(number_of_elements, Const::<1>); + let mut vec = OVector::, D>::zeros_generic( + number_of_elements, + Const::<1>, + ); + let mut vec_conj = OVector::, D>::zeros_generic( + number_of_elements, + Const::<1>, + ); for r in 0..number_of_elements_value { - 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()); + 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; } - c+=1; + c += 1; } (Some(eigenvalues), left_eigenvectors, eigenvectors) } } } - } /* diff --git a/nalgebra-lapack/tests/linalg/cholesky.rs b/nalgebra-lapack/tests/linalg/cholesky.rs index 632347b87..0bf74dd43 100644 --- a/nalgebra-lapack/tests/linalg/cholesky.rs +++ b/nalgebra-lapack/tests/linalg/cholesky.rs @@ -58,8 +58,8 @@ proptest! { let sol1 = chol.solve(&b1).unwrap(); let sol2 = chol.solve(&b2).unwrap(); - prop_assert!(relative_eq!(m * sol1, b1, epsilon = 1.0e-7)); - prop_assert!(relative_eq!(m * sol2, b2, epsilon = 1.0e-7)); + prop_assert!(relative_eq!(m * sol1, b1, epsilon = 1.0e-4)); + prop_assert!(relative_eq!(m * sol2, b2, epsilon = 1.0e-4)); } } @@ -84,7 +84,7 @@ proptest! { let id1 = &m * &m1; let id2 = &m1 * &m; - prop_assert!(id1.is_identity(1.0e-5) && id2.is_identity(1.0e-5)) + prop_assert!(id1.is_identity(1.0e-4) && id2.is_identity(1.0e-4)) } } } diff --git a/nalgebra-lapack/tests/linalg/complex_eigen.rs b/nalgebra-lapack/tests/linalg/complex_eigen.rs index 108694708..aa3474b96 100644 --- a/nalgebra-lapack/tests/linalg/complex_eigen.rs +++ b/nalgebra-lapack/tests/linalg/complex_eigen.rs @@ -1,19 +1,47 @@ -use std::cmp; - -use na::{Matrix3}; +use na::Matrix3; use nalgebra_lapack::Eigen; +use num_complex::Complex; -use crate::proptest::*; -use proptest::{prop_assert, proptest}; +#[test] +fn complex_eigen() { + let m = Matrix3::::new( + 4.0 / 5.0, + -3.0 / 5.0, + 0.0, + 3.0 / 5.0, + 4.0 / 5.0, + 0.0, + 1.0, + 2.0, + 2.0, + ); + let eigen = Eigen::new(m, true, true).expect("Eigen Creation Failed!"); + let (some_eigenvalues, some_left_vec, some_right_vec) = eigen.get_complex_elements(); + let eigenvalues = some_eigenvalues.expect("Eigenvalues Failed"); + let _left_eigenvectors = some_left_vec.expect("Left Eigenvectors Failed"); + let eigenvectors = some_right_vec.expect("Right Eigenvectors Failed"); -proptest! { - //#[test] - // fn complex_eigen() { - // let n = cmp::max(1, cmp::min(n, 10)); - // let m = DMatrix::::new_random(n, n); - // let eig = SymmetricEigen::new(m.clone()); - // let recomp = eig.recompose(); - // prop_assert!(relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5)) - // } + assert_relative_eq!( + eigenvalues[0].re, + Complex::::new(4.0 / 5.0, 3.0 / 5.0).re + ); + assert_relative_eq!( + eigenvalues[0].im, + Complex::::new(4.0 / 5.0, 3.0 / 5.0).im + ); + assert_relative_eq!( + eigenvalues[1].re, + Complex::::new(4.0 / 5.0, -3.0 / 5.0).re + ); + assert_relative_eq!( + eigenvalues[1].im, + Complex::::new(4.0 / 5.0, -3.0 / 5.0).im + ); + assert_relative_eq!(eigenvectors[0][0].re, -12.0 / 32.7871926215100059134410999); + assert_relative_eq!(eigenvectors[0][0].im, -9.0 / 32.7871926215100059134410999); + assert_relative_eq!(eigenvectors[0][1].re, -9.0 / 32.7871926215100059134410999); + assert_relative_eq!(eigenvectors[0][1].im, 12.0 / 32.7871926215100059134410999); + assert_relative_eq!(eigenvectors[0][2].re, 25.0 / 32.7871926215100059134410999); + assert_relative_eq!(eigenvectors[0][2].im, 0.0); } diff --git a/nalgebra-lapack/tests/linalg/lu.rs b/nalgebra-lapack/tests/linalg/lu.rs index 9665964e0..b9d452087 100644 --- a/nalgebra-lapack/tests/linalg/lu.rs +++ b/nalgebra-lapack/tests/linalg/lu.rs @@ -51,10 +51,10 @@ proptest! { let tr_sol1 = lup.solve_transpose(&b1).unwrap(); let tr_sol2 = lup.solve_transpose(&b2).unwrap(); - prop_assert!(relative_eq!(&m * sol1, b1, epsilon = 1.0e-7)); - prop_assert!(relative_eq!(&m * sol2, b2, epsilon = 1.0e-7)); - prop_assert!(relative_eq!(m.transpose() * tr_sol1, b1, epsilon = 1.0e-7)); - prop_assert!(relative_eq!(m.transpose() * tr_sol2, b2, epsilon = 1.0e-7)); + prop_assert!(relative_eq!(&m * sol1, b1, epsilon = 1.0e-5)); + prop_assert!(relative_eq!(&m * sol2, b2, epsilon = 1.0e-5)); + prop_assert!(relative_eq!(m.transpose() * tr_sol1, b1, epsilon = 1.0e-5)); + prop_assert!(relative_eq!(m.transpose() * tr_sol2, b2, epsilon = 1.0e-5)); } #[test] @@ -68,10 +68,10 @@ proptest! { let tr_sol1 = lup.solve_transpose(&b1).unwrap(); let tr_sol2 = lup.solve_transpose(&b2).unwrap(); - prop_assert!(relative_eq!(m * sol1, b1, epsilon = 1.0e-7)); - prop_assert!(relative_eq!(m * sol2, b2, epsilon = 1.0e-7)); - prop_assert!(relative_eq!(m.transpose() * tr_sol1, b1, epsilon = 1.0e-7)); - prop_assert!(relative_eq!(m.transpose() * tr_sol2, b2, epsilon = 1.0e-7)); + prop_assert!(relative_eq!(m * sol1, b1, epsilon = 1.0e-5)); + prop_assert!(relative_eq!(m * sol2, b2, epsilon = 1.0e-5)); + prop_assert!(relative_eq!(m.transpose() * tr_sol1, b1, epsilon = 1.0e-5)); + prop_assert!(relative_eq!(m.transpose() * tr_sol2, b2, epsilon = 1.0e-5)); } #[test] diff --git a/nalgebra-lapack/tests/linalg/mod.rs b/nalgebra-lapack/tests/linalg/mod.rs index 8fc8deebb..92425293b 100644 --- a/nalgebra-lapack/tests/linalg/mod.rs +++ b/nalgebra-lapack/tests/linalg/mod.rs @@ -1,4 +1,5 @@ mod cholesky; +mod complex_eigen; mod generalized_eigenvalues; mod lu; mod qr; @@ -7,4 +8,3 @@ mod real_eigensystem; mod schur; mod svd; mod symmetric_eigen; -mod complex_eigen; diff --git a/nalgebra-lapack/tests/linalg/real_eigensystem.rs b/nalgebra-lapack/tests/linalg/real_eigensystem.rs index 3d1c91ebf..599d1b2a2 100644 --- a/nalgebra-lapack/tests/linalg/real_eigensystem.rs +++ b/nalgebra-lapack/tests/linalg/real_eigensystem.rs @@ -13,30 +13,36 @@ proptest! { let m = DMatrix::::new_random(n, n); if let Some(eig) = Eigen::new(m.clone(), true, true) { - let eigvals = DMatrix::from_diagonal(&eig.eigenvalues); - let transformed_eigvectors = &m * eig.eigenvectors.as_ref().unwrap(); - let scaled_eigvectors = eig.eigenvectors.as_ref().unwrap() * &eigvals; - - let transformed_left_eigvectors = m.transpose() * eig.left_eigenvectors.as_ref().unwrap(); - let scaled_left_eigvectors = eig.left_eigenvectors.as_ref().unwrap() * &eigvals; - - prop_assert!(relative_eq!(transformed_eigvectors, scaled_eigvectors, epsilon = 1.0e-7)); - prop_assert!(relative_eq!(transformed_left_eigvectors, scaled_left_eigvectors, epsilon = 1.0e-7)); + // TODO: test the complex case too. + if eig.eigenvalues_are_real() { + let eigvals = DMatrix::from_diagonal(&eig.eigenvalues_re); + let transformed_eigvectors = &m * eig.eigenvectors.as_ref().unwrap(); + let scaled_eigvectors = eig.eigenvectors.as_ref().unwrap() * &eigvals; + + let transformed_left_eigvectors = m.transpose() * eig.left_eigenvectors.as_ref().unwrap(); + let scaled_left_eigvectors = eig.left_eigenvectors.as_ref().unwrap() * &eigvals; + + prop_assert!(relative_eq!(transformed_eigvectors, scaled_eigvectors, epsilon = 1.0e-5)); + prop_assert!(relative_eq!(transformed_left_eigvectors, scaled_left_eigvectors, epsilon = 1.0e-5)); + } } } #[test] fn eigensystem_static(m in matrix4()) { if let Some(eig) = Eigen::new(m, true, true) { - let eigvals = Matrix4::from_diagonal(&eig.eigenvalues); - let transformed_eigvectors = m * eig.eigenvectors.unwrap(); - let scaled_eigvectors = eig.eigenvectors.unwrap() * eigvals; - - let transformed_left_eigvectors = m.transpose() * eig.left_eigenvectors.unwrap(); - let scaled_left_eigvectors = eig.left_eigenvectors.unwrap() * eigvals; - - prop_assert!(relative_eq!(transformed_eigvectors, scaled_eigvectors, epsilon = 1.0e-7)); - prop_assert!(relative_eq!(transformed_left_eigvectors, scaled_left_eigvectors, epsilon = 1.0e-7)); + // TODO: test the complex case too. + if eig.eigenvalues_are_real() { + let eigvals = Matrix4::from_diagonal(&eig.eigenvalues_re); + let transformed_eigvectors = m * eig.eigenvectors.unwrap(); + let scaled_eigvectors = eig.eigenvectors.unwrap() * eigvals; + + let transformed_left_eigvectors = m.transpose() * eig.left_eigenvectors.unwrap(); + let scaled_left_eigvectors = eig.left_eigenvectors.unwrap() * eigvals; + + prop_assert!(relative_eq!(transformed_eigvectors, scaled_eigvectors, epsilon = 1.0e-5)); + prop_assert!(relative_eq!(transformed_left_eigvectors, scaled_left_eigvectors, epsilon = 1.0e-5)); + } } } } diff --git a/nalgebra-lapack/tests/linalg/schur.rs b/nalgebra-lapack/tests/linalg/schur.rs index 0fd1cc33e..ee7bad902 100644 --- a/nalgebra-lapack/tests/linalg/schur.rs +++ b/nalgebra-lapack/tests/linalg/schur.rs @@ -11,14 +11,17 @@ proptest! { let n = cmp::max(1, cmp::min(n, 10)); let m = DMatrix::::new_random(n, n); - let (vecs, vals) = Schur::new(m.clone()).unpack(); - - prop_assert!(relative_eq!(&vecs * vals * vecs.transpose(), m, epsilon = 1.0e-7)) + if let Some(schur) = Schur::try_new(m.clone()) { + let (vecs, vals) = schur.unpack(); + prop_assert!(relative_eq!(&vecs * vals * vecs.transpose(), m, epsilon = 1.0e-5)) + } } #[test] fn schur_static(m in matrix4()) { - let (vecs, vals) = Schur::new(m.clone()).unpack(); - prop_assert!(relative_eq!(vecs * vals * vecs.transpose(), m, epsilon = 1.0e-7)) + if let Some(schur) = Schur::try_new(m.clone()) { + let (vecs, vals) = schur.unpack(); + prop_assert!(relative_eq!(vecs * vals * vecs.transpose(), m, epsilon = 1.0e-5)) + } } } From 24e85b9883e1bd4793dc9be327e5972a4934e0a0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Crozet?= Date: Sun, 30 Oct 2022 17:22:46 +0100 Subject: [PATCH 20/20] Reset nalgebra-lapack cargo.toml to its previous defaults --- nalgebra-lapack/Cargo.toml | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/nalgebra-lapack/Cargo.toml b/nalgebra-lapack/Cargo.toml index 3f165f76d..16baf4b71 100644 --- a/nalgebra-lapack/Cargo.toml +++ b/nalgebra-lapack/Cargo.toml @@ -22,7 +22,7 @@ proptest-support = [ "nalgebra/proptest-support" ] arbitrary = [ "nalgebra/arbitrary" ] # For BLAS/LAPACK -default = ["openblas"] +default = ["netlib"] openblas = ["lapack-src/openblas"] netlib = ["lapack-src/netlib"] accelerate = ["lapack-src/accelerate"] @@ -36,7 +36,6 @@ simba = "0.7" serde = { version = "1.0", features = [ "derive" ], optional = true } lapack = { version = "0.19", default-features = false } lapack-src = { version = "0.8", default-features = false } -openblas-src = {version = "*", features = ["static"]} # clippy = "*" [dev-dependencies]