diff --git a/nalgebra-lapack/src/qz.rs b/nalgebra-lapack/src/qz.rs index ea775ea6..ee0e6208 100644 --- a/nalgebra-lapack/src/qz.rs +++ b/nalgebra-lapack/src/qz.rs @@ -42,11 +42,11 @@ where { alphar: OVector, alphai: OVector, - beta: OVector, - vsl: OMatrix, - s: OMatrix, - vsr: OMatrix, - t: OMatrix, + beta: OVector, + vsl: OMatrix, + s: OMatrix, + vsr: OMatrix, + t: OMatrix, } impl Copy for QZ @@ -176,36 +176,19 @@ where (self.vsl, self.s, self.t, self.vsr) } - /// computes the generalized eigenvalues + /// outputs the unprocessed (almost) version of generalized eigenvalues ((alphar, alpai), beta) + /// straight from LAPACK #[must_use] - pub fn eigenvalues(&self) -> OVector, D> + pub fn raw_eigenvalues(&self) -> OVector<(Complex, T), D> where - DefaultAllocator: Allocator, D>, + DefaultAllocator: Allocator<(Complex, T), D>, { - let mut out = Matrix::zeros_generic(self.t.shape_generic().0, Const::<1>); + let mut out = Matrix::from_element_generic(self.vsl.shape_generic().0, Const::<1>, (Complex::zero(), T::RealField::zero())); for i in 0..out.len() { - out[i] = if self.beta[i].clone().abs() < T::RealField::default_epsilon() { - Complex::zero() - } else { - let mut cr = self.alphar[i].clone(); - let mut ci = self.alphai[i].clone(); - let b = self.beta[i].clone(); - - if cr.clone().abs() < T::RealField::default_epsilon() { - cr = T::RealField::zero() - } else { - cr = cr / b.clone() - }; - - if ci.clone().abs() < T::RealField::default_epsilon() { - ci = T::RealField::zero() - } else { - ci = ci / b - }; - - Complex::new(cr, ci) - } + out[i] = (Complex::new(self.alphar[i].clone(), + self.alphai[i].clone()), + self.beta[i].clone()) } out diff --git a/nalgebra-lapack/tests/linalg/qz.rs b/nalgebra-lapack/tests/linalg/qz.rs index d7fe4132..6f9cf7f8 100644 --- a/nalgebra-lapack/tests/linalg/qz.rs +++ b/nalgebra-lapack/tests/linalg/qz.rs @@ -1,5 +1,7 @@ -use na::DMatrix; -use nl::{GE, QZ}; +use na::{DMatrix, EuclideanNorm, Norm}; +use nl::QZ; +use num_complex::Complex; +use simba::scalar::ComplexField; use std::cmp; use crate::proptest::*; @@ -14,28 +16,59 @@ proptest! { let qz = QZ::new(a.clone(), b.clone()); let (vsl,s,t,vsr) = qz.clone().unpack(); - let eigenvalues = qz.eigenvalues(); - - let ge = GE::new(a.clone(), b.clone()); - let eigenvalues2 = ge.eigenvalues(); + let eigenvalues = qz.raw_eigenvalues(); prop_assert!(relative_eq!(&vsl * s * vsr.transpose(), a.clone(), epsilon = 1.0e-7)); prop_assert!(relative_eq!(vsl * t * vsr.transpose(), b.clone(), epsilon = 1.0e-7)); - prop_assert!(eigenvalues == eigenvalues2); + + let a_condition_no = a.clone().try_inverse().and_then(|x| Some(EuclideanNorm.norm(&x)* EuclideanNorm.norm(&a))); + let b_condition_no = b.clone().try_inverse().and_then(|x| Some(EuclideanNorm.norm(&x)* EuclideanNorm.norm(&b))); + + if a_condition_no.unwrap_or(200000.0) < 5.0 && b_condition_no.unwrap_or(200000.0) < 5.0 { + let a_c = a.clone().map(|x| Complex::new(x, 0.0)); + let b_c = b.clone().map(|x| Complex::new(x, 0.0)); + + + for (alpha,beta) in eigenvalues.iter() { + let l_a = a_c.clone() * Complex::new(*beta, 0.0); + let l_b = b_c.clone() * *alpha; + + prop_assert!( + relative_eq!( + (&l_a - &l_b).determinant().modulus(), + 0.0, + epsilon = 1.0e-7)); + + }; + }; } #[test] fn qz_static(a in matrix4(), b in matrix4()) { let qz = QZ::new(a.clone(), b.clone()); - let ge = GE::new(a.clone(), b.clone()); let (vsl,s,t,vsr) = qz.unpack(); - let eigenvalues = qz.eigenvalues(); - let eigenvalues2 = ge.eigenvalues(); + let eigenvalues = qz.raw_eigenvalues(); prop_assert!(relative_eq!(&vsl * s * vsr.transpose(), a, epsilon = 1.0e-7)); prop_assert!(relative_eq!(vsl * t * vsr.transpose(), b, epsilon = 1.0e-7)); - prop_assert!(eigenvalues == eigenvalues2); + let a_condition_no = a.clone().try_inverse().and_then(|x| Some(EuclideanNorm.norm(&x)* EuclideanNorm.norm(&a))); + let b_condition_no = b.clone().try_inverse().and_then(|x| Some(EuclideanNorm.norm(&x)* EuclideanNorm.norm(&b))); + if a_condition_no.unwrap_or(200000.0) < 5.0 && b_condition_no.unwrap_or(200000.0) < 5.0 { + let a_c =a.clone().map(|x| Complex::new(x, 0.0)); + let b_c = b.clone().map(|x| Complex::new(x, 0.0)); + + for (alpha,beta) in eigenvalues.iter() { + let l_a = a_c.clone() * Complex::new(*beta, 0.0); + let l_b = b_c.clone() * *alpha; + + prop_assert!( + relative_eq!( + (&l_a - &l_b).determinant().modulus(), + 0.0, + epsilon = 1.0e-7)); + } + }; } }