From 89e8b4975932fab4a4a5adc8932a9022657a4cc8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Crozet?= Date: Sun, 6 Aug 2017 19:41:33 +0200 Subject: [PATCH] nalgebra-lapack: add computation of complex eigenvalues. Also renames RealEigensystem -> Eigen --- nalgebra-lapack/src/eigen.rs | 61 ++++++++++++++++--- nalgebra-lapack/src/symmetric_eigen.rs | 8 +-- .../tests/linalg/real_eigensystem.rs | 6 +- 3 files changed, 59 insertions(+), 16 deletions(-) diff --git a/nalgebra-lapack/src/eigen.rs b/nalgebra-lapack/src/eigen.rs index 648fd88a..51812fb5 100644 --- a/nalgebra-lapack/src/eigen.rs +++ b/nalgebra-lapack/src/eigen.rs @@ -1,4 +1,5 @@ use num::Zero; +use num_complex::Complex; use alga::general::Real; @@ -11,7 +12,7 @@ use na::allocator::Allocator; use lapack::fortran as interface; /// Eigendecomposition of a real square matrix with real eigenvalues. -pub struct RealEigensystem +pub struct Eigen where DefaultAllocator: Allocator + Allocator { pub eigenvalues: VectorN, @@ -20,14 +21,14 @@ pub struct RealEigensystem } -impl RealEigensystem +impl Eigen where DefaultAllocator: Allocator + Allocator { /// Computes the eigenvalues and eigenvectors of the square matrix `m`. /// /// If `eigenvectors` is `false` then, the eigenvectors are not computed explicitly. pub fn new(mut m: MatrixN, left_eigenvectors: bool, eigenvectors: bool) - -> Option> { + -> Option> { assert!(m.is_square(), "Unable to compute the eigenvalue decomposition of a non-square matrix."); @@ -67,7 +68,7 @@ impl RealEigensystem lapack_check!(info); if wi.iter().all(|e| e.is_zero()) { - return Some(RealEigensystem { + return Some(Eigen { eigenvalues: wr, left_eigenvectors: Some(vl), eigenvectors: Some(vr) }) } @@ -81,7 +82,7 @@ impl RealEigensystem lapack_check!(info); if wi.iter().all(|e| e.is_zero()) { - return Some(RealEigensystem { + return Some(Eigen { eigenvalues: wr, left_eigenvectors: Some(vl), eigenvectors: None }); } @@ -95,7 +96,7 @@ impl RealEigensystem lapack_check!(info); if wi.iter().all(|e| e.is_zero()) { - return Some(RealEigensystem { + return Some(Eigen { eigenvalues: wr, left_eigenvectors: None, eigenvectors: Some(vr) }); } @@ -107,7 +108,7 @@ impl RealEigensystem lapack_check!(info); if wi.iter().all(|e| e.is_zero()) { - return Some(RealEigensystem { + return Some(Eigen { eigenvalues: wr, left_eigenvectors: None, eigenvectors: None }); } @@ -117,6 +118,48 @@ impl RealEigensystem None } + /// The complex eigenvalues of the given matrix. + /// + /// Panics if the eigenvalue computation does not converge. + pub fn complex_eigenvalues(mut m: MatrixN) -> VectorN, D> + where DefaultAllocator: Allocator, D> { + assert!(m.is_square(), "Unable to compute the eigenvalue decomposition of a non-square matrix."); + + let (nrows, ncols) = m.data.shape(); + let n = nrows.value(); + + let lda = n as i32; + + let mut wr = unsafe { Matrix::new_uninitialized_generic(nrows, U1) }; + let mut wi = unsafe { Matrix::new_uninitialized_generic(nrows, U1) }; + + + let mut info = 0; + let mut placeholder1 = [ N::zero() ]; + let mut placeholder2 = [ N::zero() ]; + + let lwork = N::xgeev_work_size(b'N', b'N', n as i32, m.as_mut_slice(), lda, + wr.as_mut_slice(), wi.as_mut_slice(), &mut placeholder1, + n as i32, &mut placeholder2, n as i32, &mut info); + + lapack_panic!(info); + + let mut work = unsafe { ::uninitialized_vec(lwork as usize) }; + + N::xgeev(b'N', b'N', n as i32, m.as_mut_slice(), lda, wr.as_mut_slice(), + wi.as_mut_slice(), &mut placeholder1, 1 as i32, &mut placeholder2, + 1 as i32, &mut work, lwork, &mut info); + lapack_panic!(info); + + let mut res = unsafe { Matrix::new_uninitialized_generic(nrows, U1) }; + + for i in 0 .. res.len() { + res[i] = Complex::new(wr[i], wi[i]); + } + + res + } + /// The determinant of the decomposed matrix. #[inline] pub fn determinant(&self) -> N { @@ -138,7 +181,7 @@ impl RealEigensystem * Lapack functions dispatch. * */ -pub trait RealEigensystemScalar: Scalar { +pub trait EigenScalar: Scalar { fn xgeev(jobvl: u8, jobvr: u8, n: i32, a: &mut [Self], lda: i32, wr: &mut [Self], wi: &mut [Self], vl: &mut [Self], ldvl: i32, vr: &mut [Self], ldvr: i32, @@ -150,7 +193,7 @@ pub trait RealEigensystemScalar: Scalar { macro_rules! real_eigensystem_scalar_impl ( ($N: ty, $xgeev: path) => ( - impl RealEigensystemScalar for $N { + impl EigenScalar for $N { #[inline] fn xgeev(jobvl: u8, jobvr: u8, n: i32, a: &mut [Self], lda: i32, wr: &mut [Self], wi: &mut [Self], diff --git a/nalgebra-lapack/src/symmetric_eigen.rs b/nalgebra-lapack/src/symmetric_eigen.rs index e0911c55..4338118f 100644 --- a/nalgebra-lapack/src/symmetric_eigen.rs +++ b/nalgebra-lapack/src/symmetric_eigen.rs @@ -11,7 +11,7 @@ use na::allocator::Allocator; use lapack::fortran as interface; -/// Eigendecomposition of a real square matrix with real eigenvalues. +/// SymmetricEigendecomposition of a real square matrix with real eigenvalues. pub struct SymmetricEigen where DefaultAllocator: Allocator + Allocator { @@ -20,7 +20,7 @@ pub struct SymmetricEigen } -impl SymmetricEigen +impl SymmetricEigen where DefaultAllocator: Allocator + Allocator { @@ -113,7 +113,7 @@ impl SymmetricEigen * Lapack functions dispatch. * */ -pub trait RealEigensystemScalar: Scalar { +pub trait SymmetricEigenScalar: Scalar { fn xsyev(jobz: u8, uplo: u8, n: i32, a: &mut [Self], lda: i32, w: &mut [Self], work: &mut [Self], lwork: i32, info: &mut i32); fn xsyev_work_size(jobz: u8, uplo: u8, n: i32, a: &mut [Self], lda: i32, info: &mut i32) -> i32; @@ -121,7 +121,7 @@ pub trait RealEigensystemScalar: Scalar { macro_rules! real_eigensystem_scalar_impl ( ($N: ty, $xsyev: path) => ( - impl RealEigensystemScalar for $N { + impl SymmetricEigenScalar for $N { #[inline] fn xsyev(jobz: u8, uplo: u8, n: i32, a: &mut [Self], lda: i32, w: &mut [Self], work: &mut [Self], lwork: i32, info: &mut i32) { diff --git a/nalgebra-lapack/tests/linalg/real_eigensystem.rs b/nalgebra-lapack/tests/linalg/real_eigensystem.rs index a616bbfd..70a3d787 100644 --- a/nalgebra-lapack/tests/linalg/real_eigensystem.rs +++ b/nalgebra-lapack/tests/linalg/real_eigensystem.rs @@ -1,6 +1,6 @@ use std::cmp; -use nl::RealEigensystem; +use nl::Eigen; use na::{DMatrix, Matrix4}; quickcheck!{ @@ -9,7 +9,7 @@ quickcheck!{ let n = cmp::min(n, 25); let m = DMatrix::::new_random(n, n); - match RealEigensystem::new(m.clone(), true, true) { + match Eigen::new(m.clone(), true, true) { Some(eig) => { let eigvals = DMatrix::from_diagonal(&eig.eigenvalues); let transformed_eigvectors = &m * eig.eigenvectors.as_ref().unwrap(); @@ -30,7 +30,7 @@ quickcheck!{ } fn eigensystem_static(m: Matrix4) -> bool { - match RealEigensystem::new(m, true, true) { + match Eigen::new(m, true, true) { Some(eig) => { let eigvals = Matrix4::from_diagonal(&eig.eigenvalues); let transformed_eigvectors = m * eig.eigenvectors.unwrap();