1 #[cfg(feature = "serde-serialize-no-std")]
2 use serde::{Deserialize, Serialize};
3 
4 use crate::allocator::Allocator;
5 use crate::base::{DefaultAllocator, OMatrix, OVector};
6 use crate::dimension::{Const, DimDiff, DimSub, U1};
7 use simba::scalar::ComplexField;
8 
9 use crate::linalg::householder;
10 use crate::Matrix;
11 use std::mem::MaybeUninit;
12 
13 /// Hessenberg decomposition of a general matrix.
14 #[cfg_attr(feature = "serde-serialize-no-std", derive(Serialize, Deserialize))]
15 #[cfg_attr(
16     feature = "serde-serialize-no-std",
17     serde(bound(serialize = "DefaultAllocator: Allocator<T, D, D> +
18                            Allocator<T, DimDiff<D, U1>>,
19          OMatrix<T, D, D>: Serialize,
20          OVector<T, DimDiff<D, U1>>: Serialize"))
21 )]
22 #[cfg_attr(
23     feature = "serde-serialize-no-std",
24     serde(bound(deserialize = "DefaultAllocator: Allocator<T, D, D> +
25                            Allocator<T, DimDiff<D, U1>>,
26          OMatrix<T, D, D>: Deserialize<'de>,
27          OVector<T, DimDiff<D, U1>>: Deserialize<'de>"))
28 )]
29 #[derive(Clone, Debug)]
30 pub struct Hessenberg<T: ComplexField, D: DimSub<U1>>
31 where
32     DefaultAllocator: Allocator<T, D, D> + Allocator<T, DimDiff<D, U1>>,
33 {
34     hess: OMatrix<T, D, D>,
35     subdiag: OVector<T, DimDiff<D, U1>>,
36 }
37 
38 impl<T: ComplexField, D: DimSub<U1>> Copy for Hessenberg<T, D>
39 where
40     DefaultAllocator: Allocator<T, D, D> + Allocator<T, DimDiff<D, U1>>,
41     OMatrix<T, D, D>: Copy,
42     OVector<T, DimDiff<D, U1>>: Copy,
43 {
44 }
45 
46 impl<T: ComplexField, D: DimSub<U1>> Hessenberg<T, D>
47 where
48     DefaultAllocator: Allocator<T, D, D> + Allocator<T, D> + Allocator<T, DimDiff<D, U1>>,
49 {
50     /// Computes the Hessenberg decomposition using householder reflections.
new(hess: OMatrix<T, D, D>) -> Self51     pub fn new(hess: OMatrix<T, D, D>) -> Self {
52         let mut work = Matrix::zeros_generic(hess.shape_generic().0, Const::<1>);
53         Self::new_with_workspace(hess, &mut work)
54     }
55 
56     /// Computes the Hessenberg decomposition using householder reflections.
57     ///
58     /// The workspace containing `D` elements must be provided but its content does not have to be
59     /// initialized.
new_with_workspace(mut hess: OMatrix<T, D, D>, work: &mut OVector<T, D>) -> Self60     pub fn new_with_workspace(mut hess: OMatrix<T, D, D>, work: &mut OVector<T, D>) -> Self {
61         assert!(
62             hess.is_square(),
63             "Cannot compute the hessenberg decomposition of a non-square matrix."
64         );
65 
66         let dim = hess.shape_generic().0;
67 
68         assert!(
69             dim.value() != 0,
70             "Cannot compute the hessenberg decomposition of an empty matrix."
71         );
72         assert_eq!(
73             dim.value(),
74             work.len(),
75             "Hessenberg: invalid workspace size."
76         );
77 
78         if dim.value() == 0 {
79             return Hessenberg {
80                 hess,
81                 subdiag: Matrix::zeros_generic(dim.sub(Const::<1>), Const::<1>),
82             };
83         }
84 
85         let mut subdiag = Matrix::uninit(dim.sub(Const::<1>), Const::<1>);
86 
87         for ite in 0..dim.value() - 1 {
88             subdiag[ite] = MaybeUninit::new(householder::clear_column_unchecked(
89                 &mut hess,
90                 ite,
91                 1,
92                 Some(work),
93             ));
94         }
95 
96         // Safety: subdiag is now fully initialized.
97         let subdiag = unsafe { subdiag.assume_init() };
98         Hessenberg { hess, subdiag }
99     }
100 
101     /// Retrieves `(q, h)` with `q` the orthogonal matrix of this decomposition and `h` the
102     /// hessenberg matrix.
103     #[inline]
unpack(self) -> (OMatrix<T, D, D>, OMatrix<T, D, D>)104     pub fn unpack(self) -> (OMatrix<T, D, D>, OMatrix<T, D, D>) {
105         let q = self.q();
106 
107         (q, self.unpack_h())
108     }
109 
110     /// Retrieves the upper trapezoidal submatrix `H` of this decomposition.
111     #[inline]
unpack_h(mut self) -> OMatrix<T, D, D>112     pub fn unpack_h(mut self) -> OMatrix<T, D, D> {
113         let dim = self.hess.nrows();
114         self.hess.fill_lower_triangle(T::zero(), 2);
115         self.hess
116             .slice_mut((1, 0), (dim - 1, dim - 1))
117             .set_partial_diagonal(
118                 self.subdiag
119                     .iter()
120                     .map(|e| T::from_real(e.clone().modulus())),
121             );
122         self.hess
123     }
124 
125     // TODO: add a h that moves out of self.
126     /// Retrieves the upper trapezoidal submatrix `H` of this decomposition.
127     ///
128     /// This is less efficient than `.unpack_h()` as it allocates a new matrix.
129     #[inline]
130     #[must_use]
h(&self) -> OMatrix<T, D, D>131     pub fn h(&self) -> OMatrix<T, D, D> {
132         let dim = self.hess.nrows();
133         let mut res = self.hess.clone();
134         res.fill_lower_triangle(T::zero(), 2);
135         res.slice_mut((1, 0), (dim - 1, dim - 1))
136             .set_partial_diagonal(
137                 self.subdiag
138                     .iter()
139                     .map(|e| T::from_real(e.clone().modulus())),
140             );
141         res
142     }
143 
144     /// Computes the orthogonal matrix `Q` of this decomposition.
145     #[must_use]
q(&self) -> OMatrix<T, D, D>146     pub fn q(&self) -> OMatrix<T, D, D> {
147         householder::assemble_q(&self.hess, self.subdiag.as_slice())
148     }
149 
150     #[doc(hidden)]
hess_internal(&self) -> &OMatrix<T, D, D>151     pub fn hess_internal(&self) -> &OMatrix<T, D, D> {
152         &self.hess
153     }
154 }
155