From 23950400e42200e2a6e4320f3779d69010186eb7 Mon Sep 17 00:00:00 2001 From: metric-space Date: Mon, 24 Jan 2022 23:56:44 -0500 Subject: [PATCH] New wrapper for generalized eigenvalues and associated eigenvectors via LAPACK routines sggev/dggev --- .../src/generalized_eigenvalues.rs | 339 ++++++++++++++++++ nalgebra-lapack/src/lib.rs | 2 + .../tests/linalg/generalized_eigenvalues.rs | 59 +++ nalgebra-lapack/tests/linalg/mod.rs | 1 + 4 files changed, 401 insertions(+) create mode 100644 nalgebra-lapack/src/generalized_eigenvalues.rs create mode 100644 nalgebra-lapack/tests/linalg/generalized_eigenvalues.rs diff --git a/nalgebra-lapack/src/generalized_eigenvalues.rs b/nalgebra-lapack/src/generalized_eigenvalues.rs new file mode 100644 index 00000000..6332f2db --- /dev/null +++ b/nalgebra-lapack/src/generalized_eigenvalues.rs @@ -0,0 +1,339 @@ +#[cfg(feature = "serde-serialize")] +use serde::{Deserialize, Serialize}; + +use num::Zero; +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 lapack; + +/// Generalized eigenvalues and generalized eigenvectors(left and right) of a pair of N*N square matrices. +#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))] +#[cfg_attr( + feature = "serde-serialize", + serde( + bound(serialize = "DefaultAllocator: Allocator + Allocator, + OVector: Serialize, + OMatrix: Serialize") + ) +)] +#[cfg_attr( + feature = "serde-serialize", + serde( + bound(deserialize = "DefaultAllocator: Allocator + Allocator, + OVector: Serialize, + OMatrix: Deserialize<'de>") + ) +)] +#[derive(Clone, Debug)] +pub struct GE +where + DefaultAllocator: Allocator + Allocator, +{ + alphar: OVector, + alphai: OVector, + beta: OVector, + vsl: OMatrix, + vsr: OMatrix, +} + +impl Copy for GE +where + DefaultAllocator: Allocator + Allocator, + OMatrix: Copy, + OVector: Copy, +{ +} + +impl GE +where + DefaultAllocator: Allocator + Allocator, +{ + /// Attempts to compute the generalized eigenvalues (and eigenvectors) via the raw returns from LAPACK's + /// dggev and sggev routines + /// + /// For each e in generalized eigenvalues and the associated eigenvectors e_l and e_r (left andf right) + /// it satisfies e_l*a = e*e_l*b and a*e_r = e*b*e_r + /// + /// Panics if the method did not converge. + pub fn new(a: OMatrix, b: OMatrix) -> Self { + Self::try_new(a, b).expect("Calculation of generalized eigenvalues failed.") + } + + /// Attempts to compute the generalized eigenvalues (and eigenvectors) via the raw returns from LAPACK's + /// dggev and sggev routines + /// + /// For each e in generalized eigenvalues and the associated eigenvectors e_l and e_r (left andf right) + /// it satisfies e_l*a = e*e_l*b and a*e_r = e*b*e_r + /// + /// Returns `None` if the method did not converge. + pub fn try_new(mut a: OMatrix, mut b: OMatrix) -> Option { + assert!( + a.is_square() && b.is_square(), + "Unable to compute the generalized eigenvalues of non-square matrices." + ); + + assert!( + a.shape_generic() == b.shape_generic(), + "Unable to compute the generalized eigenvalues of two square matrices of different dimensions." + ); + + let (nrows, ncols) = a.shape_generic(); + let n = nrows.value(); + + let mut info = 0; + + let mut alphar = Matrix::zeros_generic(nrows, Const::<1>); + let mut alphai = Matrix::zeros_generic(nrows, Const::<1>); + let mut beta = Matrix::zeros_generic(nrows, Const::<1>); + let mut vsl = Matrix::zeros_generic(nrows, ncols); + let mut vsr = Matrix::zeros_generic(nrows, ncols); + + let lwork = T::xggev_work_size( + b'V', + b'V', + n as i32, + a.as_mut_slice(), + n as i32, + b.as_mut_slice(), + n as i32, + alphar.as_mut_slice(), + alphai.as_mut_slice(), + beta.as_mut_slice(), + vsl.as_mut_slice(), + n as i32, + vsr.as_mut_slice(), + n as i32, + &mut info, + ); + lapack_check!(info); + + let mut work = vec![T::zero(); lwork as usize]; + + T::xggev( + b'V', + b'V', + n as i32, + a.as_mut_slice(), + n as i32, + b.as_mut_slice(), + n as i32, + alphar.as_mut_slice(), + alphai.as_mut_slice(), + beta.as_mut_slice(), + vsl.as_mut_slice(), + n as i32, + vsr.as_mut_slice(), + n as i32, + &mut work, + lwork, + &mut info, + ); + lapack_check!(info); + + Some(GE { + alphar, + alphai, + beta, + vsl, + vsr, + }) + } + + /// Calculates the generalized eigenvectors (left and right) associated with the generalized eigenvalues + pub fn eigenvectors(self) -> (OMatrix, D, D>, OMatrix, D, D>) + where + DefaultAllocator: Allocator, D, D> + Allocator, D>, + { + let n = self.vsl.shape().0; + let mut l = self + .vsl + .clone() + .map(|x| Complex::new(x, T::RealField::zero())); + let mut r = self + .vsr + .clone() + .map(|x| Complex::new(x, T::RealField::zero())); + + let eigenvalues = &self.eigenvalues(); + + let mut ll; + let mut c = 0; + while c < n { + if eigenvalues[c].im.abs() > T::RealField::default_epsilon() && c + 1 < n && { + let e_conj = eigenvalues[c].conj(); + let e = eigenvalues[c + 1]; + ((e_conj.re - e.re).abs() < T::RealField::default_epsilon()) + && ((e_conj.im - e.im).abs() < T::RealField::default_epsilon()) + } { + ll = l.column(c + 1).into_owned(); + l.column_mut(c).zip_apply(&ll, |r, i| { + *r = Complex::new(r.re.clone(), i.re); + }); + ll.copy_from(&l.column(c)); + l.column_mut(c + 1).zip_apply(&ll, |r, i| { + *r = i.conj(); + }); + + ll.copy_from(&r.column(c + 1)); + r.column_mut(c).zip_apply(&ll, |r, i| { + *r = Complex::new(r.re, i.re); + }); + ll.copy_from(&r.column(c)); + r.column_mut(c + 1).zip_apply(&ll, |r, i| { + *r = i.conj(); + }); + + c += 2; + } else { + c += 1; + } + } + + (l, r) + } + + /// computes the generalized eigenvalues + #[must_use] + pub fn eigenvalues(&self) -> OVector, D> + where + DefaultAllocator: Allocator, D>, + { + let mut out = Matrix::zeros_generic(self.vsl.shape_generic().0, Const::<1>); + + 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 + } +} + +/* + * + * Lapack functions dispatch. + * + */ +/// Trait implemented by scalars for which Lapack implements the RealField GE decomposition. +pub trait GEScalar: Scalar { + #[allow(missing_docs)] + fn xggev( + jobvsl: u8, + jobvsr: u8, + n: i32, + a: &mut [Self], + lda: i32, + b: &mut [Self], + ldb: i32, + alphar: &mut [Self], + alphai: &mut [Self], + beta: &mut [Self], + vsl: &mut [Self], + ldvsl: i32, + vsr: &mut [Self], + ldvsr: i32, + work: &mut [Self], + lwork: i32, + info: &mut i32, + ); + + #[allow(missing_docs)] + fn xggev_work_size( + jobvsl: u8, + jobvsr: u8, + n: i32, + a: &mut [Self], + lda: i32, + b: &mut [Self], + ldb: i32, + alphar: &mut [Self], + alphai: &mut [Self], + beta: &mut [Self], + vsl: &mut [Self], + ldvsl: i32, + vsr: &mut [Self], + ldvsr: i32, + info: &mut i32, + ) -> i32; +} + +macro_rules! real_eigensystem_scalar_impl ( + ($N: ty, $xggev: path) => ( + impl GEScalar for $N { + #[inline] + fn xggev(jobvsl: u8, + jobvsr: u8, + n: i32, + a: &mut [$N], + lda: i32, + b: &mut [$N], + ldb: i32, + alphar: &mut [$N], + alphai: &mut [$N], + beta : &mut [$N], + vsl: &mut [$N], + ldvsl: i32, + vsr: &mut [$N], + ldvsr: i32, + work: &mut [$N], + lwork: i32, + info: &mut i32) { + unsafe { $xggev(jobvsl, jobvsr, n, a, lda, b, ldb, alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr, work, lwork, info); } + } + + + #[inline] + fn xggev_work_size(jobvsl: u8, + jobvsr: u8, + n: i32, + a: &mut [$N], + lda: i32, + b: &mut [$N], + ldb: i32, + alphar: &mut [$N], + alphai: &mut [$N], + beta : &mut [$N], + vsl: &mut [$N], + ldvsl: i32, + vsr: &mut [$N], + ldvsr: i32, + info: &mut i32) + -> i32 { + let mut work = [ Zero::zero() ]; + let lwork = -1 as i32; + + unsafe { $xggev(jobvsl, jobvsr, n, a, lda, b, ldb, alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr, &mut work, lwork, info); } + ComplexHelper::real_part(work[0]) as i32 + } + } + ) +); + +real_eigensystem_scalar_impl!(f32, lapack::sggev); +real_eigensystem_scalar_impl!(f64, lapack::dggev); diff --git a/nalgebra-lapack/src/lib.rs b/nalgebra-lapack/src/lib.rs index e89ab160..dec4daac 100644 --- a/nalgebra-lapack/src/lib.rs +++ b/nalgebra-lapack/src/lib.rs @@ -87,6 +87,7 @@ mod hessenberg; mod lu; mod qr; mod qz; +mod generalized_eigenvalues; mod schur; mod svd; mod symmetric_eigen; @@ -99,6 +100,7 @@ pub use self::hessenberg::Hessenberg; pub use self::lu::{LUScalar, LU}; pub use self::qr::QR; pub use self::qz::QZ; +pub use self::generalized_eigenvalues::GE; pub use self::schur::Schur; pub use self::svd::SVD; pub use self::symmetric_eigen::SymmetricEigen; diff --git a/nalgebra-lapack/tests/linalg/generalized_eigenvalues.rs b/nalgebra-lapack/tests/linalg/generalized_eigenvalues.rs new file mode 100644 index 00000000..275691c8 --- /dev/null +++ b/nalgebra-lapack/tests/linalg/generalized_eigenvalues.rs @@ -0,0 +1,59 @@ +use na::dimension::{Const, Dynamic}; +use na::{DMatrix, EuclideanNorm, Norm, OMatrix}; +use nl::GE; +use num_complex::Complex; +use simba::scalar::ComplexField; +use std::cmp; + +use crate::proptest::*; +use proptest::{prop_assert, proptest}; + +proptest! { + #[test] + fn ge(n in PROPTEST_MATRIX_DIM) { + let n = cmp::max(1, cmp::min(n, 10)); + let a = DMatrix::::new_random(n, n); + let b = DMatrix::::new_random(n, n); + 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) < 10.0 && b_condition_no.unwrap_or(200000.0) < 10.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)); + + let ge = GE::new(a.clone(), b.clone()); + let (vsl,vsr) = ge.clone().eigenvectors(); + let eigenvalues = ge.clone().eigenvalues(); + + for i in 0..n { + let left_eigenvector = &vsl.column(i); + prop_assert!(relative_eq!((left_eigenvector.transpose()*&a_c - left_eigenvector.transpose()*&b_c*eigenvalues[i]).map(|x| x.modulus()), OMatrix::zeros_generic(Const::<1>,Dynamic::new(n)) ,epsilon = 1.0e-7)); + + let right_eigenvector = &vsr.column(i); + prop_assert!(relative_eq!((&a_c*right_eigenvector - &b_c*right_eigenvector*eigenvalues[i]).map(|x| x.modulus()), OMatrix::zeros_generic(Dynamic::new(n), Const::<1>) ,epsilon = 1.0e-7)); + }; + }; + } + + #[test] + fn ge_static(a in matrix4(), b in matrix4()) { + 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) < 10.0 && b_condition_no.unwrap_or(200000.0) < 10.0{ + let ge = GE::new(a.clone(), b.clone()); + let a_c =a.clone().map(|x| Complex::new(x, 0.0)); + let b_c = b.clone().map(|x| Complex::new(x, 0.0)); + let (vsl,vsr) = ge.eigenvectors(); + let eigenvalues = ge.eigenvalues(); + + for i in 0..4 { + let left_eigenvector = &vsl.column(i); + prop_assert!(relative_eq!((left_eigenvector.transpose()*&a_c - left_eigenvector.transpose()*&b_c*eigenvalues[i]).map(|x| x.modulus()), OMatrix::zeros_generic(Const::<1>,Const::<4>) ,epsilon = 1.0e-7)); + + let right_eigenvector = &vsr.column(i); + prop_assert!(relative_eq!((&a_c*right_eigenvector - &b_c*right_eigenvector*eigenvalues[i]).map(|x| x.modulus()), OMatrix::zeros_generic(Const::<4>, Const::<1>) ,epsilon = 1.0e-7)); + }; + }; + } +} diff --git a/nalgebra-lapack/tests/linalg/mod.rs b/nalgebra-lapack/tests/linalg/mod.rs index a4043469..9fd539c4 100644 --- a/nalgebra-lapack/tests/linalg/mod.rs +++ b/nalgebra-lapack/tests/linalg/mod.rs @@ -2,6 +2,7 @@ mod cholesky; mod lu; mod qr; mod qz; +mod generalized_eigenvalues; mod real_eigensystem; mod schur; mod svd;