From 7230ae1e638893d2aa87b30f96dfd3cfa1e73f46 Mon Sep 17 00:00:00 2001 From: metric-space Date: Tue, 18 Jan 2022 22:35:11 -0500 Subject: [PATCH] First attempt at xgges (qz decomposition), passing tests. Serialization failing across many modules --- nalgebra-lapack/src/lib.rs | 2 + nalgebra-lapack/src/qz.rs | 288 ++++++++++++++++++++++++++++ nalgebra-lapack/tests/linalg/mod.rs | 1 + nalgebra-lapack/tests/linalg/qz.rs | 27 +++ 4 files changed, 318 insertions(+) create mode 100644 nalgebra-lapack/src/qz.rs create mode 100644 nalgebra-lapack/tests/linalg/qz.rs diff --git a/nalgebra-lapack/src/lib.rs b/nalgebra-lapack/src/lib.rs index 9cf0d73d..e89ab160 100644 --- a/nalgebra-lapack/src/lib.rs +++ b/nalgebra-lapack/src/lib.rs @@ -86,6 +86,7 @@ mod eigen; mod hessenberg; mod lu; mod qr; +mod qz; mod schur; mod svd; mod symmetric_eigen; @@ -97,6 +98,7 @@ pub use self::eigen::Eigen; pub use self::hessenberg::Hessenberg; pub use self::lu::{LUScalar, LU}; pub use self::qr::QR; +pub use self::qz::QZ; pub use self::schur::Schur; pub use self::svd::SVD; pub use self::symmetric_eigen::SymmetricEigen; diff --git a/nalgebra-lapack/src/qz.rs b/nalgebra-lapack/src/qz.rs new file mode 100644 index 00000000..47efc08c --- /dev/null +++ b/nalgebra-lapack/src/qz.rs @@ -0,0 +1,288 @@ +#[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; + +/// Eigendecomposition of a real square matrix with complex eigenvalues. +#[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 QZ +where + DefaultAllocator: Allocator + Allocator, +{ + alphar: OVector, + alphai: OVector, + beta: OVector, + vsl: OMatrix, + s: OMatrix, + vsr: OMatrix, + t: OMatrix +} + +impl Copy for QZ +where + DefaultAllocator: Allocator + Allocator, + OMatrix: Copy, + OVector: Copy, +{ +} + +impl QZ +where + DefaultAllocator: Allocator + Allocator, +{ + /// Computes the eigenvalues and real Schur form of the matrix `m`. + /// + /// Panics if the method did not converge. + pub fn new(a: OMatrix, b: OMatrix) -> Self { + Self::try_new(a,b).expect("Schur decomposition: convergence failed.") + } + + /// Computes the eigenvalues and real Schur form of the matrix `m`. + /// + /// 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 qz decomposition of non-square matrices." + ); + + // another assert to compare shape? + + let (nrows, ncols) = a.shape_generic(); + let n = nrows.value(); + + let lda = n as i32; + let ldb = lda.clone(); + + 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); + // Placeholders: + let mut bwork = [0i32]; + let mut unused = 0; + + let lwork = T::xgges_work_size( + b'V', + b'V', + b'N', + n as i32, + a.as_mut_slice(), + n as i32, + b.as_mut_slice(), + n as i32, + &mut unused, + 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 bwork, + &mut info, + ); + lapack_check!(info); + + let mut work = vec![T::zero(); lwork as usize]; + + T::xgges( + b'V', + b'V', + b'N', + n as i32, + a.as_mut_slice(), + n as i32, + b.as_mut_slice(), + n as i32, + &mut unused, + 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 bwork, + &mut info, + ); + lapack_check!(info); + + Some(QZ {alphar, alphai, beta, + vsl, s:a, + vsr, t:b}) + } + + /// Retrieves the unitary matrix `Q` and the upper-quasitriangular matrix `T` such that the + /// decomposed matrix equals `Q * T * Q.transpose()`. + pub fn unpack(self) -> (OMatrix, OMatrix, OMatrix, OMatrix){ + (self.vsl, self.s, self.t, self.vsr) + } + + /// computes the generalized eigenvalues + #[must_use] + pub fn eigenvalues(&self) -> OVector, D> + where + DefaultAllocator: Allocator, D>, + { + let mut out = Matrix::zeros_generic(self.t.shape_generic().0, Const::<1>); + + for i in 0..out.len() { + out[i] = Complex::new(self.alphar[i].clone()/self.beta[i].clone(), + self.alphai[i].clone()/self.beta[i].clone()) + } + + out + } +} + +/* + * + * Lapack functions dispatch. + * + */ +/// Trait implemented by scalars for which Lapack implements the RealField QZ decomposition. +pub trait QZScalar: Scalar { + #[allow(missing_docs)] + fn xgges( + jobvsl: u8, + jobvsr: u8, + sort: u8, + // select: ??? + n: i32, + a: &mut [Self], + lda: i32, + b: &mut [Self], + ldb: i32, + sdim: &mut 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, + bwork: &mut [i32], + info: &mut i32 + ); + + #[allow(missing_docs)] + fn xgges_work_size( + jobvsl: u8, + jobvsr: u8, + sort: u8, + // select: ??? + n: i32, + a: &mut [Self], + lda: i32, + b: &mut [Self], + ldb: i32, + sdim: &mut i32, + alphar: &mut [Self], + alphai: &mut [Self], + beta : &mut [Self], + vsl: &mut [Self], + ldvsl: i32, + vsr: &mut [Self], + ldvsr: i32, + bwork: &mut [i32], + info: &mut i32 + ) -> i32; +} + +macro_rules! real_eigensystem_scalar_impl ( + ($N: ty, $xgges: path) => ( + impl QZScalar for $N { + #[inline] + fn xgges(jobvsl: u8, + jobvsr: u8, + sort: u8, + // select: ??? + n: i32, + a: &mut [$N], + lda: i32, + b: &mut [$N], + ldb: i32, + sdim: &mut 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, + bwork: &mut [i32], + info: &mut i32) { + unsafe { $xgges(jobvsl, jobvsr, sort, None, n, a, lda, b, ldb, sdim, alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr, work, lwork, bwork, info); } + } + + + #[inline] + fn xgges_work_size(jobvsl: u8, + jobvsr: u8, + sort: u8, + // select: ??? + n: i32, + a: &mut [$N], + lda: i32, + b: &mut [$N], + ldb: i32, + sdim: &mut i32, + alphar: &mut [$N], + alphai: &mut [$N], + beta : &mut [$N], + vsl: &mut [$N], + ldvsl: i32, + vsr: &mut [$N], + ldvsr: i32, + bwork: &mut [i32], + info: &mut i32) + -> i32 { + let mut work = [ Zero::zero() ]; + let lwork = -1 as i32; + + unsafe { $xgges(jobvsl, jobvsr, sort, None, n, a, lda, b, ldb, sdim, alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr, &mut work, lwork, bwork, info); } + ComplexHelper::real_part(work[0]) as i32 + } + } + ) +); + +real_eigensystem_scalar_impl!(f32, lapack::sgges); +real_eigensystem_scalar_impl!(f64, lapack::dgges); diff --git a/nalgebra-lapack/tests/linalg/mod.rs b/nalgebra-lapack/tests/linalg/mod.rs index a6742217..a4043469 100644 --- a/nalgebra-lapack/tests/linalg/mod.rs +++ b/nalgebra-lapack/tests/linalg/mod.rs @@ -1,6 +1,7 @@ mod cholesky; mod lu; mod qr; +mod qz; mod real_eigensystem; mod schur; mod svd; diff --git a/nalgebra-lapack/tests/linalg/qz.rs b/nalgebra-lapack/tests/linalg/qz.rs new file mode 100644 index 00000000..6b4efbda --- /dev/null +++ b/nalgebra-lapack/tests/linalg/qz.rs @@ -0,0 +1,27 @@ +use na::DMatrix; +use nl::QZ; +use std::cmp; + +use crate::proptest::*; +use proptest::{prop_assert, proptest}; + +proptest! { + #[test] + fn qz(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 (vsl,s,t,vsr) = QZ::new(a.clone(), b.clone()).unpack(); + + 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)) + } + + #[test] + fn qz_static(a in matrix4(), b in matrix4()) { + let (vsl,s,t,vsr) = QZ::new(a.clone(), b.clone()).unpack(); + 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)) + } +}