From a439121641d4eb25c91e3204fbf15cb310bcda4c Mon Sep 17 00:00:00 2001 From: metric-space Date: Mon, 24 Jan 2022 23:58:21 -0500 Subject: [PATCH] Cleanup of QZ module and added GE's calculation of eigenvalues as a test for QZ's calculation of eigenvalues --- nalgebra-lapack/src/qz.rs | 17 +++++++++++----- nalgebra-lapack/tests/linalg/qz.rs | 31 ++++++++++++------------------ 2 files changed, 24 insertions(+), 24 deletions(-) diff --git a/nalgebra-lapack/src/qz.rs b/nalgebra-lapack/src/qz.rs index e3319452..ea775ea6 100644 --- a/nalgebra-lapack/src/qz.rs +++ b/nalgebra-lapack/src/qz.rs @@ -13,7 +13,11 @@ use na::{DefaultAllocator, Matrix, OMatrix, OVector, Scalar}; use lapack; -/// Generalized eigendecomposition of a pair of N*N square matrices. +/// QZ decomposition of a pair of N*N square matrices. +/// Retrieves the left and right matrices of Schur Vectors (VSL and VSR) +/// the upper-quasitriangular matrix `S` and upper triangular matrix `T` such that the +/// decomposed input matrix `a` equals `VSL * S * VSL.transpose()` and +/// decomposed input matrix `b` equals `VSL * T * VSL.transpose()`. #[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))] #[cfg_attr( feature = "serde-serialize", @@ -59,6 +63,11 @@ where { /// Attempts to compute the QZ decomposition of input square matrices `a` and `b`. /// + /// i.e retrieves the left and right matrices of Schur Vectors (VSL and VSR) + /// the upper-quasitriangular matrix `S` and upper triangular matrix `T` such that the + /// decomposed matrix `a` equals `VSL * S * VSL.transpose()` and + /// decomposed matrix `b` equals `VSL * T * VSL.transpose()`. + /// /// Panics if the method did not converge. pub fn new(a: OMatrix, b: OMatrix) -> Self { Self::try_new(a, b).expect("QZ decomposition: convergence failed.") @@ -154,8 +163,8 @@ where /// Retrieves the left and right matrices of Schur Vectors (VSL and VSR) /// the upper-quasitriangular matrix `S` and upper triangular matrix `T` such that the - /// decomposed matrix `A` equals `VSL * S * VSL.transpose()` and - /// decomposed matrix `B` equals `VSL * T * VSL.transpose()`. + /// decomposed input matrix `a` equals `VSL * S * VSL.transpose()` and + /// decomposed input matrix `b` equals `VSL * T * VSL.transpose()`. pub fn unpack( self, ) -> ( @@ -167,8 +176,6 @@ where (self.vsl, self.s, self.t, self.vsr) } - - /// computes the generalized eigenvalues #[must_use] pub fn eigenvalues(&self) -> OVector, D> diff --git a/nalgebra-lapack/tests/linalg/qz.rs b/nalgebra-lapack/tests/linalg/qz.rs index 84a7b030..d7fe4132 100644 --- a/nalgebra-lapack/tests/linalg/qz.rs +++ b/nalgebra-lapack/tests/linalg/qz.rs @@ -1,7 +1,5 @@ -use na::{zero, DMatrix, SMatrix}; -use nl::QZ; -use num_complex::Complex; -use simba::scalar::ComplexField; +use na::DMatrix; +use nl::{GE, QZ}; use std::cmp; use crate::proptest::*; @@ -16,33 +14,28 @@ proptest! { let qz = QZ::new(a.clone(), b.clone()); let (vsl,s,t,vsr) = qz.clone().unpack(); - //let eigenvalues = qz.eigenvalues(); - //let a_c = a.clone().map(|x| Complex::new(x, zero::())); + let eigenvalues = qz.eigenvalues(); + + let ge = GE::new(a.clone(), b.clone()); + let eigenvalues2 = ge.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)); - // spotty test that skips over the first eigenvalue which in some cases is extremely large relative to the other ones - // and fails the condition - //for i in 1..n { - // let b_c = b.clone().map(|x| eigenvalues[i]*Complex::new(x,zero::())); - // prop_assert!(relative_eq!((&a_c - &b_c).determinant().modulus(), 0.0, epsilon = 1.0e-6)); - //} + prop_assert!(eigenvalues == eigenvalues2); } #[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 a_c = a.clone().map(|x| Complex::new(x, zero::())); + let eigenvalues = qz.eigenvalues(); + let eigenvalues2 = ge.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)); - //for i in 0..4 { - // let b_c = b.clone().map(|x| eigenvalues[i]*Complex::new(x,zero::())); - // println!("{}",eigenvalues); - // prop_assert!(relative_eq!((&a_c - &b_c).determinant().modulus(), 0.0, epsilon = 1.0e-4)) - //} + prop_assert!(eigenvalues == eigenvalues2); + } }