diff --git a/src/linalg/symmetric_eigen.rs b/src/linalg/symmetric_eigen.rs index 748f94f8..424533f0 100644 --- a/src/linalg/symmetric_eigen.rs +++ b/src/linalg/symmetric_eigen.rs @@ -2,8 +2,9 @@ use num_complex::Complex; use std::ops::MulAssign; use alga::general::Real; -use core::{MatrixN, VectorN, DefaultAllocator, Matrix2, Vector2}; +use core::{MatrixN, VectorN, DefaultAllocator, Matrix2, Vector2, SquareMatrix}; use dimension::{Dim, DimSub, DimDiff, U1, U2}; +use storage::Storage; use allocator::Allocator; use linalg::givens; @@ -46,7 +47,19 @@ impl SymmetricEigen /// * `max_niter` − maximum total number of iterations performed by the algorithm. If this /// number of iteration is exceeded, `None` is returned. If `niter == 0`, then the algorithm /// continues indefinitely until convergence. - pub fn try_new(mut m: MatrixN, eps: N, max_niter: usize) -> Option + pub fn try_new(m: MatrixN, eps: N, max_niter: usize) -> Option + where D: DimSub, + DefaultAllocator: Allocator> { + Self::do_decompose(m, true, eps, max_niter).map(|(vals, vecs)| { + SymmetricEigen { + eigenvectors: vecs.unwrap(), + eigenvalues: vals + } + }) + } + + fn do_decompose(mut m: MatrixN, eigenvectors: bool, eps: N, max_niter: usize) + -> Option<(VectorN, Option>)> where D: DimSub, DefaultAllocator: Allocator> { @@ -59,15 +72,24 @@ impl SymmetricEigen m /= m_amax; } - let (mut q, mut diag, mut off_diag) = SymmetricTridiagonal::new(m).unpack(); + let (mut q, mut diag, mut off_diag); + + if eigenvectors { + let res = SymmetricTridiagonal::new(m).unpack(); + q = Some(res.0); + diag = res.1; + off_diag = res.2; + } + else { + let res = SymmetricTridiagonal::new(m).unpack_tridiagonal(); + q = None; + diag = res.0; + off_diag = res.1; + } if dim == 1 { diag *= m_amax; - - return Some(SymmetricEigen { - eigenvectors: q, - eigenvalues: diag - }); + return Some((diag, q)); } let mut niter = 0; @@ -114,7 +136,9 @@ impl SymmetricEigen off_diag[i + 1] *= rot.cos_angle(); } - rot.inverse().rotate_rows(&mut q.fixed_columns_mut::(i)); + if let Some(ref mut q) = q { + rot.inverse().rotate_rows(&mut q.fixed_columns_mut::(i)); + } } else { break; @@ -134,9 +158,11 @@ impl SymmetricEigen diag[start + 0] = eigvals[0]; diag[start + 1] = eigvals[1]; - if let Some(basis) = basis.try_normalize(eps) { - let rot = UnitComplex::new_unchecked(Complex::new(basis.x, basis.y)); - rot.rotate_rows(&mut q.fixed_columns_mut::(start)); + if let Some(ref mut q) = q { + if let Some(basis) = basis.try_normalize(eps) { + let rot = UnitComplex::new_unchecked(Complex::new(basis.x, basis.y)); + rot.rotate_rows(&mut q.fixed_columns_mut::(start)); + } } end -= 1; @@ -156,12 +182,7 @@ impl SymmetricEigen diag *= m_amax; - // Solve the remaining 2x2 subproblem. - - Some(SymmetricEigen { - eigenvectors: q, - eigenvalues: diag - }) + Some((diag, q)) } fn delimit_subproblem(diag: &VectorN, @@ -236,6 +257,28 @@ pub fn wilkinson_shift(tmm: N, tnn: N, tmn: N) -> N { } } + +/* + * + * Computations of eigenvalues for symmetric matrices. + * + */ +impl, S: Storage> SquareMatrix + where DefaultAllocator: Allocator + + Allocator + + Allocator> { + /// Computes the eigenvalues of this symmetric matrix. + /// + /// Only the lower-triangular part of the matrix is read. + pub fn symmetric_eigenvalues(&self) -> VectorN { + SymmetricEigen::do_decompose(self.clone_owned(), false, N::default_epsilon(), 0).unwrap().0 + } +} + + + + + #[cfg(test)] mod test { use core::Matrix2; diff --git a/src/linalg/symmetric_tridiagonal.rs b/src/linalg/symmetric_tridiagonal.rs index 39327c41..5ab4cca0 100644 --- a/src/linalg/symmetric_tridiagonal.rs +++ b/src/linalg/symmetric_tridiagonal.rs @@ -71,6 +71,14 @@ impl> SymmetricTridiagonal (q, diag, self.off_diagonal) } + /// Retrieve the diagonal, and off diagonal elements of this decomposition. + pub fn unpack_tridiagonal(self) -> (VectorN, VectorN>) + where DefaultAllocator: Allocator { + let diag = self.diagonal(); + + (diag, self.off_diagonal) + } + /// The diagonal components of this decomposition. pub fn diagonal(&self) -> VectorN where DefaultAllocator: Allocator {