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::linalg::householder;
10use crate::Matrix;
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#[derive(Clone, Debug)]
30pub struct Hessenberg<T: ComplexField, D: DimSub<U1>>
31where
32    DefaultAllocator: Allocator<D, D> + Allocator<DimDiff<D, U1>>,
33{
34    hess: OMatrix<T, D, D>,
35    subdiag: OVector<T, DimDiff<D, U1>>,
36}
37
38impl<T: ComplexField, D: DimSub<U1>> Copy for Hessenberg<T, D>
39where
40    DefaultAllocator: Allocator<D, D> + Allocator<DimDiff<D, U1>>,
41    OMatrix<T, D, D>: Copy,
42    OVector<T, DimDiff<D, U1>>: Copy,
43{
44}
45
46impl<T: ComplexField, D: DimSub<U1>> Hessenberg<T, D>
47where
48    DefaultAllocator: Allocator<D, D> + Allocator<D> + Allocator<DimDiff<D, U1>>,
49{
50    /// Computes the Hessenberg decomposition using householder reflections.
51    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.
60    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]
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]
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            .view_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]
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.view_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]
146    pub fn q(&self) -> OMatrix<T, D, D> {
147        householder::assemble_q(&self.hess, self.subdiag.as_slice())
148    }
149
150    #[doc(hidden)]
151    pub fn hess_internal(&self) -> &OMatrix<T, D, D> {
152        &self.hess
153    }
154}
nalgebra/linalg/hessenberg.rs

nalgebra/linalg/
hessenberg.rs
