From 6eb0d8a78679a47f755640d178d9a3d9ed9c1e3f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Crozet?= Date: Sun, 6 Aug 2017 16:12:05 +0200 Subject: [PATCH] nalgebra-lapack: add Symmetric eigensystems. --- nalgebra-lapack/src/lib.rs | 2 + nalgebra-lapack/src/symmetric_eigen.rs | 146 ++++++++++++++++++ nalgebra-lapack/tests/linalg/mod.rs | 1 + .../tests/linalg/symmetric_eigen.rs | 20 +++ 4 files changed, 169 insertions(+) create mode 100644 nalgebra-lapack/src/symmetric_eigen.rs create mode 100644 nalgebra-lapack/tests/linalg/symmetric_eigen.rs diff --git a/nalgebra-lapack/src/lib.rs b/nalgebra-lapack/src/lib.rs index 9ad8f41e..41e6e8d5 100644 --- a/nalgebra-lapack/src/lib.rs +++ b/nalgebra-lapack/src/lib.rs @@ -16,6 +16,7 @@ extern crate nalgebra as na; mod lapack_check; mod svd; mod eigen; +mod symmetric_eigen; mod cholesky; mod lu; mod qr; @@ -27,6 +28,7 @@ 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::symmetric_eigen::SymmetricEigen; pub use self::qr::QR; pub use self::hessenberg::Hessenberg; diff --git a/nalgebra-lapack/src/symmetric_eigen.rs b/nalgebra-lapack/src/symmetric_eigen.rs new file mode 100644 index 00000000..e0911c55 --- /dev/null +++ b/nalgebra-lapack/src/symmetric_eigen.rs @@ -0,0 +1,146 @@ +use num::Zero; +use std::ops::MulAssign; + +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 SymmetricEigen + where DefaultAllocator: Allocator + + Allocator { + pub eigenvalues: VectorN, + pub eigenvectors: MatrixN, +} + + +impl SymmetricEigen + where DefaultAllocator: Allocator + + Allocator { + + /// Computes the eigenvalues and eigenvectors of the symmetric matrix `m`. + /// + /// Only the lower-triangular part of `m` is read. If `eigenvectors` is `false` then, the + /// eigenvectors are not computed explicitly. Panics if the method did not converge. + pub fn new(m: MatrixN) -> Self { + let (vals, vecs) = Self::do_decompose(m, true).expect("SymmetricEigen: convergence failure."); + SymmetricEigen { eigenvalues: vals, eigenvectors: vecs.unwrap() } + } + + /// Computes the eigenvalues and eigenvectors of the symmetric matrix `m`. + /// + /// Only the lower-triangular part of `m` is read. If `eigenvectors` is `false` then, the + /// eigenvectors are not computed explicitly. Returns `None` if the method did not converge. + pub fn try_new(m: MatrixN) -> Option { + Self::do_decompose(m, true).map(|(vals, vecs)| { + SymmetricEigen { eigenvalues: vals, eigenvectors: vecs.unwrap() } + }) + } + + fn do_decompose(mut m: MatrixN, eigenvectors: bool) -> Option<(VectorN, Option>)> { + assert!(m.is_square(), "Unable to compute the eigenvalue decomposition of a non-square matrix."); + + let jobz = if eigenvectors { b'V' } else { b'N' }; + + let (nrows, ncols) = m.data.shape(); + let n = nrows.value(); + + let lda = n as i32; + + let mut values = unsafe { Matrix::new_uninitialized_generic(nrows, U1) }; + let mut info = 0; + + let lwork = N::xsyev_work_size(jobz, b'L', n as i32, m.as_mut_slice(), lda, &mut info); + lapack_check!(info); + + let mut work = unsafe { ::uninitialized_vec(lwork as usize) }; + + N::xsyev(jobz, b'L', n as i32, m.as_mut_slice(), lda, values.as_mut_slice(), &mut work, lwork, &mut info); + lapack_check!(info); + + let vectors = if eigenvectors { Some(m) } else { None }; + Some((values, vectors)) + } + + /// Computes only the eigenvalues of the input matrix. + /// + /// Panics if the method does not converge. + pub fn eigenvalues(mut m: MatrixN) -> VectorN { + Self::do_decompose(m, false).expect("SymmetricEigen eigenvalues: convergence failure.").0 + } + + /// Computes only the eigenvalues of the input matrix. + /// + /// Returns `None` if the method does not converge. + pub fn try_eigenvalues(mut m: MatrixN) -> Option> { + Self::do_decompose(m, false).map(|res| res.0) + } + + /// The determinant of the decomposed matrix. + #[inline] + pub fn determinant(&self) -> N { + let mut det = N::one(); + for e in self.eigenvalues.iter() { + det *= *e; + } + + det + } + + /// Rebuild the original matrix. + /// + /// This is useful if some of the eigenvalues have been manually modified. + pub fn recompose(&self) -> MatrixN { + let mut u_t = self.eigenvectors.clone(); + for i in 0 .. self.eigenvalues.len() { + let val = self.eigenvalues[i]; + u_t.column_mut(i).mul_assign(val); + } + u_t.transpose_mut(); + &self.eigenvectors * u_t + } +} + + +/* + * + * Lapack functions dispatch. + * + */ +pub trait RealEigensystemScalar: 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; +} + +macro_rules! real_eigensystem_scalar_impl ( + ($N: ty, $xsyev: path) => ( + impl RealEigensystemScalar 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) { + $xsyev(jobz, uplo, n, a, lda, w, work, lwork, info) + } + + + #[inline] + fn xsyev_work_size(jobz: u8, uplo: u8, n: i32, a: &mut [Self], lda: i32, info: &mut i32) -> i32 { + let mut work = [ Zero::zero() ]; + let mut w = [ Zero::zero() ]; + let lwork = -1 as i32; + + $xsyev(jobz, uplo, n, a, lda, &mut w, &mut work, lwork, info); + ComplexHelper::real_part(work[0]) as i32 + } + } + ) +); + +real_eigensystem_scalar_impl!(f32, interface::ssyev); +real_eigensystem_scalar_impl!(f64, interface::dsyev); diff --git a/nalgebra-lapack/tests/linalg/mod.rs b/nalgebra-lapack/tests/linalg/mod.rs index 8caadde1..0e2afd52 100644 --- a/nalgebra-lapack/tests/linalg/mod.rs +++ b/nalgebra-lapack/tests/linalg/mod.rs @@ -1,4 +1,5 @@ mod real_eigensystem; +mod symmetric_eigen; mod cholesky; mod lu; mod qr; diff --git a/nalgebra-lapack/tests/linalg/symmetric_eigen.rs b/nalgebra-lapack/tests/linalg/symmetric_eigen.rs new file mode 100644 index 00000000..a1ebdfa2 --- /dev/null +++ b/nalgebra-lapack/tests/linalg/symmetric_eigen.rs @@ -0,0 +1,20 @@ +use std::cmp; + +use nl::SymmetricEigen; +use na::{DMatrix, Matrix4}; + +quickcheck!{ + fn symmetric_eigen(n: usize) -> bool { + let n = cmp::max(1, cmp::min(n, 10)); + let m = DMatrix::::new_random(n, n); + let eig = SymmetricEigen::new(m.clone()); + let recomp = eig.recompose(); + relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5) + } + + fn symmetric_eigen_static(m: Matrix4) -> bool { + let eig = SymmetricEigen::new(m); + let recomp = eig.recompose(); + relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5) + } +}