2017-08-14 01:53:04 +08:00
|
|
|
|
#[cfg(feature = "serde-serialize")]
|
2018-10-22 13:00:10 +08:00
|
|
|
|
use serde::{Deserialize, Serialize};
|
2017-08-14 01:53:04 +08:00
|
|
|
|
|
2019-03-03 02:33:49 +08:00
|
|
|
|
use alga::general::Complex;
|
2019-03-23 21:29:07 +08:00
|
|
|
|
use crate::allocator::Allocator;
|
|
|
|
|
use crate::base::{DefaultAllocator, MatrixMN, MatrixN, SquareMatrix, VectorN};
|
|
|
|
|
use crate::dimension::{DimDiff, DimSub, U1};
|
|
|
|
|
use crate::storage::Storage;
|
2017-08-03 01:37:44 +08:00
|
|
|
|
|
2019-03-23 21:29:07 +08:00
|
|
|
|
use crate::linalg::householder;
|
2017-08-03 01:37:44 +08:00
|
|
|
|
|
2017-08-14 01:53:04 +08:00
|
|
|
|
/// Hessenberg decomposition of a general matrix.
|
|
|
|
|
#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
|
2018-05-19 23:15:15 +08:00
|
|
|
|
#[cfg_attr(
|
|
|
|
|
feature = "serde-serialize",
|
2018-10-22 13:00:10 +08:00
|
|
|
|
serde(bound(
|
|
|
|
|
serialize = "DefaultAllocator: Allocator<N, D, D> +
|
2017-08-14 01:53:04 +08:00
|
|
|
|
Allocator<N, DimDiff<D, U1>>,
|
2018-09-13 12:55:58 +08:00
|
|
|
|
MatrixN<N, D>: Serialize,
|
|
|
|
|
VectorN<N, DimDiff<D, U1>>: Serialize"
|
2018-10-22 13:00:10 +08:00
|
|
|
|
))
|
2018-05-19 23:15:15 +08:00
|
|
|
|
)]
|
|
|
|
|
#[cfg_attr(
|
|
|
|
|
feature = "serde-serialize",
|
2018-10-22 13:00:10 +08:00
|
|
|
|
serde(bound(
|
|
|
|
|
deserialize = "DefaultAllocator: Allocator<N, D, D> +
|
2017-08-14 01:53:04 +08:00
|
|
|
|
Allocator<N, DimDiff<D, U1>>,
|
2018-09-13 12:55:58 +08:00
|
|
|
|
MatrixN<N, D>: Deserialize<'de>,
|
|
|
|
|
VectorN<N, DimDiff<D, U1>>: Deserialize<'de>"
|
2018-10-22 13:00:10 +08:00
|
|
|
|
))
|
2018-05-19 23:15:15 +08:00
|
|
|
|
)]
|
2017-08-14 01:53:00 +08:00
|
|
|
|
#[derive(Clone, Debug)]
|
2019-03-03 02:33:49 +08:00
|
|
|
|
pub struct Hessenberg<N: Complex, D: DimSub<U1>>
|
2018-10-22 13:00:10 +08:00
|
|
|
|
where DefaultAllocator: Allocator<N, D, D> + Allocator<N, DimDiff<D, U1>>
|
2018-02-02 19:26:35 +08:00
|
|
|
|
{
|
|
|
|
|
hess: MatrixN<N, D>,
|
|
|
|
|
subdiag: VectorN<N, DimDiff<D, U1>>,
|
2017-08-03 01:37:44 +08:00
|
|
|
|
}
|
|
|
|
|
|
2019-03-03 02:33:49 +08:00
|
|
|
|
impl<N: Complex, D: DimSub<U1>> Copy for Hessenberg<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,
|
2018-10-22 13:00:10 +08:00
|
|
|
|
{}
|
2017-08-14 01:53:00 +08:00
|
|
|
|
|
2019-03-03 02:33:49 +08:00
|
|
|
|
impl<N: Complex, D: DimSub<U1>> Hessenberg<N, D>
|
2018-10-22 13:00:10 +08:00
|
|
|
|
where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D> + Allocator<N, DimDiff<D, U1>>
|
2018-02-02 19:26:35 +08:00
|
|
|
|
{
|
2017-08-03 01:37:44 +08:00
|
|
|
|
/// Computes the Hessenberg decomposition using householder reflections.
|
|
|
|
|
pub fn new(hess: MatrixN<N, D>) -> Self {
|
|
|
|
|
let mut work = unsafe { MatrixMN::new_uninitialized_generic(hess.data.shape().0, U1) };
|
|
|
|
|
Self::new_with_workspace(hess, &mut work)
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
/// Computes the Hessenberg decomposition using householder reflections.
|
|
|
|
|
///
|
2017-08-14 01:53:04 +08:00
|
|
|
|
/// The workspace containing `D` elements must be provided but its content does not have to be
|
|
|
|
|
/// initialized.
|
2017-08-03 01:37:44 +08:00
|
|
|
|
pub fn new_with_workspace(mut hess: MatrixN<N, D>, work: &mut VectorN<N, D>) -> Self {
|
2018-02-02 19:26:35 +08:00
|
|
|
|
assert!(
|
|
|
|
|
hess.is_square(),
|
|
|
|
|
"Cannot compute the hessenberg decomposition of a non-square matrix."
|
|
|
|
|
);
|
2017-08-03 01:37:44 +08:00
|
|
|
|
|
|
|
|
|
let dim = hess.data.shape().0;
|
|
|
|
|
|
2018-02-02 19:26:35 +08:00
|
|
|
|
assert!(
|
|
|
|
|
dim.value() != 0,
|
|
|
|
|
"Cannot compute the hessenberg decomposition of an empty matrix."
|
|
|
|
|
);
|
|
|
|
|
assert_eq!(
|
|
|
|
|
dim.value(),
|
|
|
|
|
work.len(),
|
|
|
|
|
"Hessenberg: invalid workspace size."
|
|
|
|
|
);
|
2017-08-03 01:37:44 +08:00
|
|
|
|
|
|
|
|
|
let mut subdiag = unsafe { MatrixMN::new_uninitialized_generic(dim.sub(U1), U1) };
|
|
|
|
|
|
|
|
|
|
if dim.value() == 0 {
|
|
|
|
|
return Hessenberg { hess, subdiag };
|
|
|
|
|
}
|
|
|
|
|
|
2018-02-02 19:26:35 +08:00
|
|
|
|
for ite in 0..dim.value() - 1 {
|
2017-08-03 01:37:44 +08:00
|
|
|
|
householder::clear_column_unchecked(&mut hess, &mut subdiag[ite], ite, 1, Some(work));
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
Hessenberg { hess, subdiag }
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
/// Retrieves `(q, h)` with `q` the orthogonal matrix of this decomposition and `h` the
|
|
|
|
|
/// hessenberg matrix.
|
|
|
|
|
#[inline]
|
2019-03-23 21:13:00 +08:00
|
|
|
|
pub fn unpack(self) -> (MatrixN<N, D>, MatrixN<N, D>) {
|
2017-08-03 01:37:44 +08:00
|
|
|
|
let q = self.q();
|
|
|
|
|
|
|
|
|
|
(q, self.unpack_h())
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
/// Retrieves the upper trapezoidal submatrix `H` of this decomposition.
|
|
|
|
|
#[inline]
|
2019-03-23 21:13:00 +08:00
|
|
|
|
pub fn unpack_h(mut self) -> MatrixN<N, D> {
|
2017-08-03 01:37:44 +08:00
|
|
|
|
let dim = self.hess.nrows();
|
|
|
|
|
self.hess.fill_lower_triangle(N::zero(), 2);
|
2018-02-02 19:26:35 +08:00
|
|
|
|
self.hess
|
|
|
|
|
.slice_mut((1, 0), (dim - 1, dim - 1))
|
2019-03-23 21:13:00 +08:00
|
|
|
|
.set_partial_diagonal(self.subdiag.iter().map(|e| N::from_real(e.modulus())));
|
2017-08-03 01:37:44 +08:00
|
|
|
|
self.hess
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
// FIXME: add a h that moves out of self.
|
|
|
|
|
/// Retrieves the upper trapezoidal submatrix `H` of this decomposition.
|
|
|
|
|
///
|
2017-08-14 01:53:04 +08:00
|
|
|
|
/// This is less efficient than `.unpack_h()` as it allocates a new matrix.
|
2017-08-03 01:37:44 +08:00
|
|
|
|
#[inline]
|
2019-03-23 21:13:00 +08:00
|
|
|
|
pub fn h(&self) -> MatrixN<N, D> {
|
2017-08-03 01:37:44 +08:00
|
|
|
|
let dim = self.hess.nrows();
|
|
|
|
|
let mut res = self.hess.clone();
|
|
|
|
|
res.fill_lower_triangle(N::zero(), 2);
|
2018-02-02 19:26:35 +08:00
|
|
|
|
res.slice_mut((1, 0), (dim - 1, dim - 1))
|
2019-03-23 21:13:00 +08:00
|
|
|
|
.set_partial_diagonal(self.subdiag.iter().map(|e| N::from_real(e.modulus())));
|
2017-08-03 01:37:44 +08:00
|
|
|
|
res
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
/// Computes the orthogonal matrix `Q` of this decomposition.
|
|
|
|
|
pub fn q(&self) -> MatrixN<N, D> {
|
2019-03-19 19:00:10 +08:00
|
|
|
|
householder::assemble_q(&self.hess, self.subdiag.as_slice())
|
2017-08-03 01:37:44 +08:00
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
#[doc(hidden)]
|
|
|
|
|
pub fn hess_internal(&self) -> &MatrixN<N, D> {
|
|
|
|
|
&self.hess
|
|
|
|
|
}
|
|
|
|
|
}
|
2017-08-14 01:52:46 +08:00
|
|
|
|
|
2019-03-03 02:33:49 +08:00
|
|
|
|
impl<N: Complex, D: DimSub<U1>, S: Storage<N, D, D>> SquareMatrix<N, D, S>
|
2018-10-22 13:00:10 +08:00
|
|
|
|
where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D> + Allocator<N, DimDiff<D, U1>>
|
2018-02-02 19:26:35 +08:00
|
|
|
|
{
|
2017-08-14 01:52:46 +08:00
|
|
|
|
/// Computes the Hessenberg decomposition of this matrix using householder reflections.
|
|
|
|
|
pub fn hessenberg(self) -> Hessenberg<N, D> {
|
|
|
|
|
Hessenberg::new(self.into_owned())
|
|
|
|
|
}
|
|
|
|
|
}
|