New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
QZ-decomposition #1067
Merged
sebcrozet
merged 31 commits into
dimforge:dev
from
metric-space:qz-decomposition-lapack
Apr 30, 2022
Merged
QZ-decomposition #1067
Changes from all commits
Commits
Show all changes
31 commits
Select commit
Hold shift + click to select a range
7230ae1
First attempt at xgges (qz decomposition), passing tests. Serializati…
metric-space 6f7ef38
Format file
metric-space 769f20c
Comments more tailored to QZ
metric-space b2c6c6b
Add non-naive way of calculate generalized eigenvalue, write spotty t…
metric-space 6a28306
Commented out failing tests, refactored checks for almost zeroes
metric-space 7e9345d
Correction for not calculating absolurte value
metric-space 7d8fb3d
New wrapper for generalized eigenvalues and associated eigenvectors v…
metric-space 748848f
Cleanup of QZ module and added GE's calculation of eigenvalues as a t…
metric-space ae35d1c
New code and modified tests for generalized_eigenvalues
metric-space 162bb32
New code and modified tests for qz
metric-space d88337e
Formatting
metric-space 55fdd84
Formatting
metric-space 4362f00
Added comment on logic
metric-space 4038e66
Correction to keep naming of left and right eigenvector matrices cons…
metric-space d5069f3
Removed extra memory allocation for buffer (now redundant)
metric-space a4de6a8
Corrected deserialization term in serialization impls
metric-space 497a4d8
Correction in eigenvector matrices build up algorithm
metric-space fb0cb51
Remove condition number, tests pass without. Add proper test generato…
metric-space b7fe6b9
Name change
metric-space 91b7e05
Change name of test generator
metric-space 7510d48
Doc string corrections
metric-space e52f117
Change name of copied macro base
metric-space 5e10ca4
Add another case for when eigenvalues should be mapped to zero. Make …
metric-space c8a920f
Minimal post-processing and fix to documentation
metric-space 3413ab7
Correct typos, move doc portion to comment and fix borrow to clone
metric-space 4413a35
Fix doc
metric-space adf50a6
Fix formatting
metric-space 2743eef
Add in explicit type of matrix element to module overview docs
metric-space 80a844a
Update check for zero
metric-space bc31012
Add newline
metric-space ff2d431
Remove repeated docs
metric-space File filter
Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,350 @@ | ||
#[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 real square matrices. | ||
/// | ||
/// Each generalized eigenvalue (lambda) satisfies determinant(A - lambda*B) = 0 | ||
/// | ||
/// The right eigenvector v(j) corresponding to the eigenvalue lambda(j) | ||
/// of (A,B) satisfies | ||
/// | ||
/// A * v(j) = lambda(j) * B * v(j). | ||
/// | ||
/// 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). | ||
#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))] | ||
#[cfg_attr( | ||
feature = "serde-serialize", | ||
serde( | ||
bound(serialize = "DefaultAllocator: Allocator<T, D, D> + Allocator<T, D>, | ||
OVector<T, D>: Serialize, | ||
OMatrix<T, D, D>: Serialize") | ||
) | ||
)] | ||
#[cfg_attr( | ||
feature = "serde-serialize", | ||
serde( | ||
bound(deserialize = "DefaultAllocator: Allocator<T, D, D> + Allocator<T, D>, | ||
OVector<T, D>: Deserialize<'de>, | ||
OMatrix<T, D, D>: Deserialize<'de>") | ||
) | ||
)] | ||
#[derive(Clone, Debug)] | ||
pub struct GeneralizedEigen<T: Scalar, D: Dim> | ||
where | ||
DefaultAllocator: Allocator<T, D> + Allocator<T, D, D>, | ||
{ | ||
alphar: OVector<T, D>, | ||
alphai: OVector<T, D>, | ||
beta: OVector<T, D>, | ||
vsl: OMatrix<T, D, D>, | ||
vsr: OMatrix<T, D, D>, | ||
} | ||
|
||
impl<T: Scalar + Copy, D: Dim> Copy for GeneralizedEigen<T, D> | ||
where | ||
DefaultAllocator: Allocator<T, D, D> + Allocator<T, D>, | ||
OMatrix<T, D, D>: Copy, | ||
OVector<T, D>: Copy, | ||
{ | ||
} | ||
|
||
impl<T: GeneralizedEigenScalar + RealField + Copy, D: Dim> GeneralizedEigen<T, D> | ||
where | ||
DefaultAllocator: Allocator<T, D, D> + Allocator<T, D>, | ||
{ | ||
/// Attempts to compute the generalized eigenvalues, and left and right associated eigenvectors | ||
/// via the raw returns from LAPACK's dggev and sggev routines | ||
/// | ||
/// Panics if the method did not converge. | ||
pub fn new(a: OMatrix<T, D, D>, b: OMatrix<T, D, D>) -> 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 | ||
/// | ||
/// Returns `None` if the method did not converge. | ||
pub fn try_new(mut a: OMatrix<T, D, D>, mut b: OMatrix<T, D, D>) -> Option<Self> { | ||
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(GeneralizedEigen { | ||
alphar, | ||
alphai, | ||
beta, | ||
vsl, | ||
vsr, | ||
}) | ||
} | ||
|
||
/// Calculates the generalized eigenvectors (left and right) associated with the generalized eigenvalues | ||
/// | ||
/// Outputs two matrices. | ||
/// The first output matrix contains the left eigenvectors of the generalized eigenvalues | ||
/// as columns. | ||
/// The second matrix contains the right eigenvectors of the generalized eigenvalues | ||
/// as columns. | ||
pub fn eigenvectors(&self) -> (OMatrix<Complex<T>, D, D>, OMatrix<Complex<T>, D, D>) | ||
where | ||
DefaultAllocator: | ||
Allocator<Complex<T>, D, D> + Allocator<Complex<T>, D> + Allocator<(Complex<T>, T), D>, | ||
{ | ||
/* | ||
How the eigenvectors are built up: | ||
|
||
Since the input entries are all real, the generalized eigenvalues if complex come in pairs | ||
as a consequence of the [complex conjugate root thorem](https://en.wikipedia.org/wiki/Complex_conjugate_root_theorem) | ||
The Lapack routine output reflects this by expecting the user to unpack the real and 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). | ||
*/ | ||
|
||
let n = self.vsl.shape().0; | ||
|
||
let mut l = self.vsl.map(|x| Complex::new(x, T::RealField::zero())); | ||
|
||
let mut r = self.vsr.map(|x| Complex::new(x, T::RealField::zero())); | ||
|
||
let eigenvalues = self.raw_eigenvalues(); | ||
|
||
let mut c = 0; | ||
|
||
while c < n { | ||
if eigenvalues[c].0.im.abs() != T::RealField::zero() && 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()); | ||
}); | ||
l.column_mut(c + 1).zip_apply(&self.vsl.column(c), |i, r| { | ||
*i = Complex::new(r.clone(), -i.re.clone()); | ||
}); | ||
|
||
// taking care of the right eigenvector matrix | ||
r.column_mut(c).zip_apply(&self.vsr.column(c + 1), |r, i| { | ||
*r = Complex::new(r.re.clone(), i.clone()); | ||
}); | ||
r.column_mut(c + 1).zip_apply(&self.vsr.column(c), |i, r| { | ||
*i = Complex::new(r.clone(), -i.re.clone()); | ||
}); | ||
|
||
c += 2; | ||
} else { | ||
c += 1; | ||
} | ||
} | ||
|
||
(l, r) | ||
} | ||
|
||
/// Outputs the unprocessed (almost) version of generalized eigenvalues ((alphar, alphai), beta) | ||
/// straight from LAPACK | ||
#[must_use] | ||
pub fn raw_eigenvalues(&self) -> OVector<(Complex<T>, T), D> | ||
where | ||
DefaultAllocator: Allocator<(Complex<T>, T), D>, | ||
{ | ||
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] = (Complex::new(self.alphar[i], self.alphai[i]), self.beta[i]) | ||
} | ||
|
||
out | ||
} | ||
} | ||
|
||
/* | ||
* | ||
* Lapack functions dispatch. | ||
* | ||
*/ | ||
/// Trait implemented by scalars for which Lapack implements the RealField GeneralizedEigen decomposition. | ||
pub trait GeneralizedEigenScalar: 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! generalized_eigen_scalar_impl ( | ||
($N: ty, $xggev: path) => ( | ||
impl GeneralizedEigenScalar 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 | ||
} | ||
} | ||
) | ||
); | ||
|
||
generalized_eigen_scalar_impl!(f32, lapack::sggev); | ||
generalized_eigen_scalar_impl!(f64, lapack::dggev); |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Oops, something went wrong.
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Doesn't there need to be an extra newline here to make sure only the first sentence ends up in the module overview?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
re:module over-view: isn't the module overview here https://github.com/dimforge/nalgebra/pull/1067/files#diff-a056a4c6b4fb629735a73fcb65e0d5223386cb47dedb6ab80d5a1b905ebd9d25R16 ?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Sorry, not module overview, but the docs for the struct. Please render them (
cargo doc
) and see :)