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