From c8a920ff2c0a0b433e04e039502eb9d951de742e Mon Sep 17 00:00:00 2001 From: metric-space Date: Sun, 27 Feb 2022 17:17:31 -0500 Subject: [PATCH] Minimal post-processing and fix to documentation --- .../src/generalized_eigenvalues.rs | 82 ++++++++----------- 1 file changed, 36 insertions(+), 46 deletions(-) diff --git a/nalgebra-lapack/src/generalized_eigenvalues.rs b/nalgebra-lapack/src/generalized_eigenvalues.rs index dac8004c..132be1b7 100644 --- a/nalgebra-lapack/src/generalized_eigenvalues.rs +++ b/nalgebra-lapack/src/generalized_eigenvalues.rs @@ -180,25 +180,44 @@ where } /// Calculates the generalized eigenvectors (left and right) associated with the generalized eigenvalues - /// Outputs two matrices, the first one containing the left eigenvectors of the generalized eigenvalues - /// as columns and the second matrix contains the right eigenvectors of the generalized eigenvalues - /// as columns + /// Outputs two matrices. + /// The first output matix contains the left eigenvectors of the generalized eigenvalues + /// as columns. + /// The second matrix contains the right eigenvectors of the generalized eigenvalues + /// as columns. /// - /// The right eigenvector v(j) corresponding to the eigenvalue lambda(j) - /// of (A,B) satisfies + /// The right eigenvector v(j) corresponding to the eigenvalue lambda(j) + /// of (A,B) satisfies /// - /// A * v(j) = lambda(j) * B * v(j). + /// A * v(j) = lambda(j) * B * v(j) /// - /// The left eigenvector u(j) corresponding to the eigenvalue lambda(j) - /// of (A,B) satisfies + /// The left eigenvector u(j) corresponding to the eigenvalue lambda(j) + /// of (A,B) satisfies /// - /// u(j)**H * A = lambda(j) * u(j)**H * B . - /// where u(j)**H is the conjugate-transpose of u(j). + /// u(j)**H * A = lambda(j) * u(j)**H * B + /// where u(j)**H is the conjugate-transpose of u(j). + /// + /// How the eigenvectors are build up: + /// + /// Since the input entries are all real, the generalized eigenvalues if complex come in pairs + /// as a consequence of + /// The Lapack routine output reflects this by expecting the user to unpack the complex eigenvalues associated + /// eigenvectors from the real matrix output via the following procedure + /// + /// (Note: VL stands for the lapack real matrix output containing the left eigenvectors as columns, + /// VR stands for the lapack real matrix output containing the right eigenvectors as columns) + /// + /// If the j-th and (j+1)-th eigenvalues form a complex conjugate pair, + /// then + /// + /// u(j) = VL(:,j)+i*VL(:,j+1) + /// u(j+1) = VL(:,j)-i*VL(:,j+1) + /// + /// and + /// + /// u(j) = VR(:,j)+i*VR(:,j+1) + /// v(j+1) = VR(:,j)-i*VR(:,j+1). /// - /// What is going on below? - /// If the j-th and (j+1)-th eigenvalues form a complex conjugate pair, - /// then u(j) = VSL(:,j)+i*VSL(:,j+1) and u(j+1) = VSL(:,j)-i*VSL(:,j+1). - /// and then v(j) = VSR(:,j)+i*VSR(:,j+1) and v(j+1) = VSR(:,j)-i*VSR(:,j+1). pub fn eigenvectors(self) -> (OMatrix, D, D>, OMatrix, D, D>) where DefaultAllocator: @@ -216,18 +235,14 @@ where .clone() .map(|x| Complex::new(x, T::RealField::zero())); - let eigenvalues = &self.eigenvalues(); + let eigenvalues = &self.raw_eigenvalues(); let mut c = 0; let epsilon = T::RealField::default_epsilon(); while c < n { - if eigenvalues[c].im.abs() > epsilon && c + 1 < n && { - let e_conj = eigenvalues[c].conj(); - let e = eigenvalues[c + 1]; - (&e_conj.re).ulps_eq(&e.re, epsilon, 6) && (&e_conj.im).ulps_eq(&e.im, epsilon, 6) - } { + if eigenvalues[c].0.im.abs() > epsilon && c + 1 < n { // taking care of the left eigenvector matrix l.column_mut(c).zip_apply(&self.vsl.column(c + 1), |r, i| { *r = Complex::new(r.re.clone(), i.clone()); @@ -253,32 +268,7 @@ where (l, r) } - // only used for internal calculation for assembling eigenvectors based on realness of - // eigenvalues and complex-conjugate checks of subsequent non-real eigenvalues - fn eigenvalues(&self) -> OVector, D> - where - DefaultAllocator: Allocator, D>, - { - let mut out = Matrix::zeros_generic(self.vsl.shape_generic().0, Const::<1>); - - let epsilon = T::RealField::default_epsilon(); - - for i in 0..out.len() { - out[i] = if self.beta[i].clone().abs() < epsilon - || (self.alphai[i].clone().abs() < epsilon - && self.alphar[i].clone().abs() < epsilon) - { - Complex::zero() - } else { - Complex::new(self.alphar[i].clone(), self.alphai[i].clone()) - * (Complex::new(self.beta[i].clone(), T::RealField::zero()).inv()) - } - } - - out - } - - /// outputs the unprocessed (almost) version of generalized eigenvalues ((alphar, alpai), beta) + /// outputs the unprocessed (almost) version of generalized eigenvalues ((alphar, alphai), beta) /// straight from LAPACK #[must_use] pub fn raw_eigenvalues(&self) -> OVector<(Complex, T), D>