2017-08-14 01:53:04 +08:00
|
|
|
#[cfg(feature = "serde-serialize")]
|
|
|
|
use serde;
|
|
|
|
|
2017-08-03 01:37:44 +08:00
|
|
|
use alga::general::Real;
|
2018-02-02 19:26:35 +08:00
|
|
|
use core::{DefaultAllocator, MatrixMN, MatrixN, SquareMatrix, VectorN};
|
|
|
|
use dimension::{DimDiff, DimSub, U1};
|
2017-08-03 01:37:44 +08:00
|
|
|
use storage::Storage;
|
|
|
|
use allocator::Allocator;
|
|
|
|
|
|
|
|
use linalg::householder;
|
|
|
|
|
2017-08-14 01:53:04 +08:00
|
|
|
/// Tridiagonalization of a symmetric matrix.
|
|
|
|
#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
|
|
|
|
#[cfg_attr(feature = "serde-serialize",
|
2018-02-02 19:26:35 +08:00
|
|
|
serde(bound(serialize = "DefaultAllocator: Allocator<N, D, D> +
|
2017-08-14 01:53:04 +08:00
|
|
|
Allocator<N, DimDiff<D, U1>>,
|
|
|
|
MatrixN<N, D>: serde::Serialize,
|
|
|
|
VectorN<N, DimDiff<D, U1>>: serde::Serialize")))]
|
|
|
|
#[cfg_attr(feature = "serde-serialize",
|
2018-02-02 19:26:35 +08:00
|
|
|
serde(bound(deserialize = "DefaultAllocator: Allocator<N, D, D> +
|
2017-08-14 01:53:04 +08:00
|
|
|
Allocator<N, DimDiff<D, U1>>,
|
|
|
|
MatrixN<N, D>: serde::Deserialize<'de>,
|
|
|
|
VectorN<N, DimDiff<D, U1>>: serde::Deserialize<'de>")))]
|
2017-08-14 01:53:00 +08:00
|
|
|
#[derive(Clone, Debug)]
|
2017-08-03 01:37:44 +08:00
|
|
|
pub struct SymmetricTridiagonal<N: Real, D: DimSub<U1>>
|
2018-02-02 19:26:35 +08:00
|
|
|
where
|
|
|
|
DefaultAllocator: Allocator<N, D, D> + Allocator<N, DimDiff<D, U1>>,
|
|
|
|
{
|
|
|
|
tri: MatrixN<N, D>,
|
|
|
|
off_diagonal: VectorN<N, DimDiff<D, U1>>,
|
2017-08-03 01:37:44 +08:00
|
|
|
}
|
|
|
|
|
2017-08-14 01:53:04 +08:00
|
|
|
impl<N: Real, D: DimSub<U1>> Copy for SymmetricTridiagonal<N, D>
|
2018-02-02 19:26:35 +08:00
|
|
|
where
|
|
|
|
DefaultAllocator: Allocator<N, D, D> + Allocator<N, DimDiff<D, U1>>,
|
|
|
|
MatrixN<N, D>: Copy,
|
|
|
|
VectorN<N, DimDiff<D, U1>>: Copy,
|
|
|
|
{
|
|
|
|
}
|
2017-08-14 01:53:00 +08:00
|
|
|
|
2017-08-03 01:37:44 +08:00
|
|
|
impl<N: Real, D: DimSub<U1>> SymmetricTridiagonal<N, D>
|
2018-02-02 19:26:35 +08:00
|
|
|
where
|
|
|
|
DefaultAllocator: Allocator<N, D, D> + Allocator<N, DimDiff<D, U1>>,
|
|
|
|
{
|
2017-08-03 01:37:44 +08:00
|
|
|
/// Computes the tridiagonalization of the symmetric matrix `m`.
|
|
|
|
///
|
2017-08-14 01:53:04 +08:00
|
|
|
/// Only the lower-triangular part (including the diagonal) of `m` is read.
|
2017-08-03 01:37:44 +08:00
|
|
|
pub fn new(mut m: MatrixN<N, D>) -> Self {
|
|
|
|
let dim = m.data.shape().0;
|
|
|
|
|
2018-02-02 19:26:35 +08:00
|
|
|
assert!(
|
|
|
|
m.is_square(),
|
|
|
|
"Unable to compute the symmetric tridiagonal decomposition of a non-square matrix."
|
|
|
|
);
|
|
|
|
assert!(
|
|
|
|
dim.value() != 0,
|
|
|
|
"Unable to compute the symmetric tridiagonal decomposition of an empty matrix."
|
|
|
|
);
|
2017-08-03 01:37:44 +08:00
|
|
|
|
|
|
|
let mut off_diagonal = unsafe { MatrixMN::new_uninitialized_generic(dim.sub(U1), U1) };
|
2018-02-02 19:26:35 +08:00
|
|
|
let mut p = unsafe { MatrixMN::new_uninitialized_generic(dim.sub(U1), U1) };
|
2017-08-03 01:37:44 +08:00
|
|
|
|
2018-02-02 19:26:35 +08:00
|
|
|
for i in 0..dim.value() - 1 {
|
|
|
|
let mut m = m.rows_range_mut(i + 1..);
|
|
|
|
let (mut axis, mut m) = m.columns_range_pair_mut(i, i + 1..);
|
2017-08-03 01:37:44 +08:00
|
|
|
|
|
|
|
let (norm, not_zero) = householder::reflection_axis_mut(&mut axis);
|
|
|
|
off_diagonal[i] = norm;
|
|
|
|
|
|
|
|
if not_zero {
|
2018-02-02 19:26:35 +08:00
|
|
|
let mut p = p.rows_range_mut(i..);
|
|
|
|
|
2017-08-03 01:37:44 +08:00
|
|
|
p.gemv_symm(::convert(2.0), &m, &axis, N::zero());
|
|
|
|
let dot = axis.dot(&p);
|
|
|
|
p.axpy(-dot, &axis, N::one());
|
|
|
|
m.ger_symm(-N::one(), &p, &axis, N::one());
|
|
|
|
m.ger_symm(-N::one(), &axis, &p, N::one());
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
SymmetricTridiagonal {
|
2018-02-02 19:26:35 +08:00
|
|
|
tri: m,
|
|
|
|
off_diagonal: off_diagonal,
|
2017-08-03 01:37:44 +08:00
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
#[doc(hidden)]
|
|
|
|
// For debugging.
|
|
|
|
pub fn internal_tri(&self) -> &MatrixN<N, D> {
|
|
|
|
&self.tri
|
|
|
|
}
|
|
|
|
|
|
|
|
/// Retrieve the orthogonal transformation, diagonal, and off diagonal elements of this
|
|
|
|
/// decomposition.
|
|
|
|
pub fn unpack(self) -> (MatrixN<N, D>, VectorN<N, D>, VectorN<N, DimDiff<D, U1>>)
|
2018-02-02 19:26:35 +08:00
|
|
|
where
|
|
|
|
DefaultAllocator: Allocator<N, D>,
|
|
|
|
{
|
2017-08-03 01:37:44 +08:00
|
|
|
let diag = self.diagonal();
|
2018-02-02 19:26:35 +08:00
|
|
|
let q = self.q();
|
2017-08-03 01:37:44 +08:00
|
|
|
|
|
|
|
(q, diag, self.off_diagonal)
|
2017-08-06 23:04:40 +08:00
|
|
|
}
|
|
|
|
|
|
|
|
/// Retrieve the diagonal, and off diagonal elements of this decomposition.
|
|
|
|
pub fn unpack_tridiagonal(self) -> (VectorN<N, D>, VectorN<N, DimDiff<D, U1>>)
|
2018-02-02 19:26:35 +08:00
|
|
|
where
|
|
|
|
DefaultAllocator: Allocator<N, D>,
|
|
|
|
{
|
2017-08-06 23:04:40 +08:00
|
|
|
let diag = self.diagonal();
|
|
|
|
|
|
|
|
(diag, self.off_diagonal)
|
2017-08-03 01:37:44 +08:00
|
|
|
}
|
|
|
|
|
|
|
|
/// The diagonal components of this decomposition.
|
|
|
|
pub fn diagonal(&self) -> VectorN<N, D>
|
2018-02-02 19:26:35 +08:00
|
|
|
where
|
|
|
|
DefaultAllocator: Allocator<N, D>,
|
|
|
|
{
|
2017-08-03 01:37:44 +08:00
|
|
|
self.tri.diagonal()
|
|
|
|
}
|
|
|
|
|
|
|
|
/// The off-diagonal components of this decomposition.
|
|
|
|
pub fn off_diagonal(&self) -> &VectorN<N, DimDiff<D, U1>>
|
2018-02-02 19:26:35 +08:00
|
|
|
where
|
|
|
|
DefaultAllocator: Allocator<N, D>,
|
|
|
|
{
|
2017-08-03 01:37:44 +08:00
|
|
|
&self.off_diagonal
|
|
|
|
}
|
|
|
|
|
|
|
|
/// Computes the orthogonal matrix `Q` of this decomposition.
|
|
|
|
pub fn q(&self) -> MatrixN<N, D> {
|
|
|
|
householder::assemble_q(&self.tri)
|
|
|
|
}
|
|
|
|
|
|
|
|
/// Recomputes the original symmetric matrix.
|
|
|
|
pub fn recompose(mut self) -> MatrixN<N, D> {
|
|
|
|
let q = self.q();
|
|
|
|
self.tri.fill_lower_triangle(N::zero(), 2);
|
|
|
|
self.tri.fill_upper_triangle(N::zero(), 2);
|
|
|
|
|
2018-02-02 19:26:35 +08:00
|
|
|
for i in 0..self.off_diagonal.len() {
|
2017-08-03 01:37:44 +08:00
|
|
|
self.tri[(i + 1, i)] = self.off_diagonal[i];
|
|
|
|
self.tri[(i, i + 1)] = self.off_diagonal[i];
|
|
|
|
}
|
|
|
|
|
|
|
|
&q * self.tri * q.transpose()
|
|
|
|
}
|
|
|
|
}
|
2017-08-14 01:52:46 +08:00
|
|
|
|
|
|
|
impl<N: Real, D: DimSub<U1>, S: Storage<N, D, D>> SquareMatrix<N, D, S>
|
2018-02-02 19:26:35 +08:00
|
|
|
where
|
|
|
|
DefaultAllocator: Allocator<N, D, D> + Allocator<N, DimDiff<D, U1>>,
|
|
|
|
{
|
2017-08-14 01:52:46 +08:00
|
|
|
/// Computes the tridiagonalization of this symmetric matrix.
|
|
|
|
///
|
2017-08-14 01:53:04 +08:00
|
|
|
/// Only the lower-triangular part (including the diagonal) of `m` is read.
|
2017-08-14 01:52:46 +08:00
|
|
|
pub fn symmetric_tridiagonalize(self) -> SymmetricTridiagonal<N, D> {
|
|
|
|
SymmetricTridiagonal::new(self.into_owned())
|
|
|
|
}
|
|
|
|
}
|