nalgebra/linalg/
hessenberg.rs

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

nalgebra/linalg/
hessenberg.rs
