diff --git a/nalgebra-lapack/src/lib.rs b/nalgebra-lapack/src/lib.rs index 41e6e8d5..5afc3399 100644 --- a/nalgebra-lapack/src/lib.rs +++ b/nalgebra-lapack/src/lib.rs @@ -21,16 +21,18 @@ mod cholesky; mod lu; mod qr; mod hessenberg; +mod schur; use num_complex::Complex; pub use self::svd::SVD; pub use self::cholesky::{Cholesky, CholeskyScalar}; pub use self::lu::{LU, LUScalar}; -pub use self::eigen::RealEigensystem; +pub use self::eigen::Eigen; pub use self::symmetric_eigen::SymmetricEigen; pub use self::qr::QR; pub use self::hessenberg::Hessenberg; +pub use self::schur::RealSchur; trait ComplexHelper { diff --git a/nalgebra-lapack/src/schur.rs b/nalgebra-lapack/src/schur.rs new file mode 100644 index 00000000..82156803 --- /dev/null +++ b/nalgebra-lapack/src/schur.rs @@ -0,0 +1,191 @@ +use num::Zero; +use num_complex::Complex; + +use alga::general::Real; + +use ::ComplexHelper; +use na::{Scalar, DefaultAllocator, Matrix, VectorN, MatrixN}; +use na::dimension::{Dim, U1}; +use na::storage::Storage; +use na::allocator::Allocator; + +use lapack::fortran as interface; + +/// Eigendecomposition of a real square matrix with real eigenvalues. +pub struct RealSchur + where DefaultAllocator: Allocator + + Allocator { + + re: VectorN, + im: VectorN, + t: MatrixN, + q: MatrixN +} + + +impl RealSchur + where DefaultAllocator: Allocator + + Allocator { + /// Computes the eigenvalues and real Schur foorm of the matrix `m`. + /// + /// Panics if the method did not converge. + pub fn new(m: MatrixN) -> Self { + Self::try_new(m).expect("RealSchur decomposition: convergence failed.") + } + + /// Computes the eigenvalues and real Schur foorm of the matrix `m`. + /// + /// Returns `None` if the method did not converge. + pub fn try_new(mut m: MatrixN) -> Option { + 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 info = 0; + + let mut wr = unsafe { Matrix::new_uninitialized_generic(nrows, U1) }; + let mut wi = unsafe { Matrix::new_uninitialized_generic(nrows, U1) }; + let mut q = unsafe { Matrix::new_uninitialized_generic(nrows, ncols) }; + // Placeholders: + let mut bwork = [ 0i32 ]; + let mut unused = 0; + + let lwork = N::xgees_work_size(b'V', b'N', n as i32, m.as_mut_slice(), lda, &mut unused, + wr.as_mut_slice(), wi.as_mut_slice(), q.as_mut_slice(), n as i32, + &mut bwork, &mut info); + lapack_check!(info); + + let mut work = unsafe { ::uninitialized_vec(lwork as usize) }; + + N::xgees(b'V', b'N', n as i32, m.as_mut_slice(), lda, &mut unused, + wr.as_mut_slice(), wi.as_mut_slice(), q.as_mut_slice(), + n as i32, &mut work, lwork, &mut bwork, &mut info); + lapack_check!(info); + + Some(RealSchur { re: wr, im: wi, t: m, q: q }) + } + + /// 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) -> (MatrixN, MatrixN) { + (self.q, self.t) + } + + /// Computes the real eigenvalues of the decomposed matrix. + /// + /// Return `None` if some eigenvalues are complex. + pub fn eigenvalues(&self) -> Option> { + if self.im.iter().all(|e| e.is_zero()) { + Some(self.re.clone()) + } + else { + None + } + } + + /// Computes the complex eigenvalues of the decomposed matrix. + pub fn complex_eigenvalues(&self) -> VectorN, D> + where DefaultAllocator: Allocator, D> { + + let mut out = unsafe { VectorN::new_uninitialized_generic(self.t.data.shape().0, U1) }; + + for i in 0 .. out.len() { + out[i] = Complex::new(self.re[i], self.im[i]) + } + + out + } +} + + +/* + * + * Lapack functions dispatch. + * + */ +pub trait EigenScalar: Scalar { + fn xgees(jobvs: u8, + sort: u8, + // select: ??? + n: i32, + a: &mut [Self], + lda: i32, + sdim: &mut i32, + wr: &mut [Self], + wi: &mut [Self], + vs: &mut [Self], + ldvs: i32, + work: &mut [Self], + lwork: i32, + bwork: &mut [i32], + info: &mut i32); + + fn xgees_work_size(jobvs: u8, + sort: u8, + // select: ??? + n: i32, + a: &mut [Self], + lda: i32, + sdim: &mut i32, + wr: &mut [Self], + wi: &mut [Self], + vs: &mut [Self], + ldvs: i32, + bwork: &mut [i32], + info: &mut i32) + -> i32; +} + +macro_rules! real_eigensystem_scalar_impl ( + ($N: ty, $xgees: path) => ( + impl EigenScalar for $N { + #[inline] + fn xgees(jobvs: u8, + sort: u8, + // select: ??? + n: i32, + a: &mut [$N], + lda: i32, + sdim: &mut i32, + wr: &mut [$N], + wi: &mut [$N], + vs: &mut [$N], + ldvs: i32, + work: &mut [$N], + lwork: i32, + bwork: &mut [i32], + info: &mut i32) { + $xgees(jobvs, sort, None, n, a, lda, sdim, wr, wi, vs, ldvs, work, lwork, bwork, info); + } + + + #[inline] + fn xgees_work_size(jobvs: u8, + sort: u8, + // select: ??? + n: i32, + a: &mut [$N], + lda: i32, + sdim: &mut i32, + wr: &mut [$N], + wi: &mut [$N], + vs: &mut [$N], + ldvs: i32, + bwork: &mut [i32], + info: &mut i32) + -> i32 { + let mut work = [ Zero::zero() ]; + let lwork = -1 as i32; + + $xgees(jobvs, sort, None, n, a, lda, sdim, wr, wi, vs, ldvs, &mut work, lwork, bwork, info); + ComplexHelper::real_part(work[0]) as i32 + } + } + ) +); + +real_eigensystem_scalar_impl!(f32, interface::sgees); +real_eigensystem_scalar_impl!(f64, interface::dgees); diff --git a/nalgebra-lapack/tests/linalg/mod.rs b/nalgebra-lapack/tests/linalg/mod.rs index 0e2afd52..f692fa88 100644 --- a/nalgebra-lapack/tests/linalg/mod.rs +++ b/nalgebra-lapack/tests/linalg/mod.rs @@ -4,3 +4,4 @@ mod cholesky; mod lu; mod qr; mod svd; +mod real_schur; diff --git a/nalgebra-lapack/tests/linalg/real_schur.rs b/nalgebra-lapack/tests/linalg/real_schur.rs new file mode 100644 index 00000000..4511d925 --- /dev/null +++ b/nalgebra-lapack/tests/linalg/real_schur.rs @@ -0,0 +1,21 @@ +use std::cmp; +use nl::RealSchur; +use na::{DMatrix, Matrix4}; + +quickcheck! { + fn schur(n: usize) -> bool { + let n = cmp::max(1, cmp::min(n, 10)); + let m = DMatrix::::new_random(n, n); + + let (vecs, vals) = RealSchur::new(m.clone()).unpack(); + + relative_eq!(&vecs * vals * vecs.transpose(), m, epsilon = 1.0e-7) + } + + fn schur_static(m: Matrix4) -> bool { + let (vecs, vals) = RealSchur::new(m.clone()).unpack(); + + relative_eq!(vecs * vals * vecs.transpose(), m, epsilon = 1.0e-7) + } +} +