nalgebra_lapack/
hessenberg.rs

1#[cfg(feature = "serde-serialize")]
2use serde::{Deserialize, Serialize};
3
4use num::Zero;
5use num_complex::Complex;
6
7use crate::ComplexHelper;
8use na::allocator::Allocator;
9use na::dimension::{Const, DimDiff, DimSub, U1};
10use na::{DefaultAllocator, Matrix, OMatrix, OVector, Scalar};
11
12use lapack;
13
14/// The Hessenberg decomposition of a general matrix.
15#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
16#[cfg_attr(
17    feature = "serde-serialize",
18    serde(bound(serialize = "DefaultAllocator: Allocator<D, D> +
19                           Allocator<DimDiff<D, U1>>,
20         OMatrix<T, D, D>: Serialize,
21         OVector<T, DimDiff<D, U1>>: Serialize"))
22)]
23#[cfg_attr(
24    feature = "serde-serialize",
25    serde(bound(deserialize = "DefaultAllocator: Allocator<D, D> +
26                           Allocator<DimDiff<D, U1>>,
27         OMatrix<T, D, D>: Deserialize<'de>,
28         OVector<T, DimDiff<D, U1>>: Deserialize<'de>"))
29)]
30#[derive(Clone, Debug)]
31pub struct Hessenberg<T: Scalar, D: DimSub<U1>>
32where
33    DefaultAllocator: Allocator<D, D> + Allocator<DimDiff<D, U1>>,
34{
35    h: OMatrix<T, D, D>,
36    tau: OVector<T, DimDiff<D, U1>>,
37}
38
39impl<T: Scalar + Copy, 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: HessenbergScalar + Zero, D: DimSub<U1>> Hessenberg<T, D>
48where
49    DefaultAllocator: Allocator<D, D> + Allocator<DimDiff<D, U1>>,
50{
51    /// Computes the hessenberg decomposition of the matrix `m`.
52    pub fn new(mut m: OMatrix<T, D, D>) -> Self {
53        let nrows = m.shape_generic().0;
54        let n = nrows.value() as i32;
55
56        assert!(
57            m.is_square(),
58            "Unable to compute the hessenberg decomposition of a non-square matrix."
59        );
60        assert!(
61            !m.is_empty(),
62            "Unable to compute the hessenberg decomposition of an empty matrix."
63        );
64
65        let mut tau = Matrix::zeros_generic(nrows.sub(Const::<1>), Const::<1>);
66
67        let mut info = 0;
68        let lwork =
69            T::xgehrd_work_size(n, 1, n, m.as_mut_slice(), n, tau.as_mut_slice(), &mut info);
70        let mut work = vec![T::zero(); lwork as usize];
71
72        lapack_panic!(info);
73
74        T::xgehrd(
75            n,
76            1,
77            n,
78            m.as_mut_slice(),
79            n,
80            tau.as_mut_slice(),
81            &mut work,
82            lwork,
83            &mut info,
84        );
85        lapack_panic!(info);
86
87        Self { h: m, tau }
88    }
89
90    /// Computes the hessenberg matrix of this decomposition.
91    #[inline]
92    #[must_use]
93    pub fn h(&self) -> OMatrix<T, D, D> {
94        let mut h = self.h.clone_owned();
95        h.fill_lower_triangle(T::zero(), 2);
96
97        h
98    }
99}
100
101impl<T: HessenbergReal + Zero, D: DimSub<U1>> Hessenberg<T, D>
102where
103    DefaultAllocator: Allocator<D, D> + Allocator<DimDiff<D, U1>>,
104{
105    /// Computes the matrices `(Q, H)` of this decomposition.
106    #[inline]
107    pub fn unpack(self) -> (OMatrix<T, D, D>, OMatrix<T, D, D>) {
108        (self.q(), self.h())
109    }
110
111    /// Computes the unitary matrix `Q` of this decomposition.
112    #[inline]
113    #[must_use]
114    pub fn q(&self) -> OMatrix<T, D, D> {
115        let n = self.h.nrows() as i32;
116        let mut q = self.h.clone_owned();
117        let mut info = 0;
118
119        let lwork =
120            T::xorghr_work_size(n, 1, n, q.as_mut_slice(), n, self.tau.as_slice(), &mut info);
121        let mut work = vec![T::zero(); lwork as usize];
122
123        T::xorghr(
124            n,
125            1,
126            n,
127            q.as_mut_slice(),
128            n,
129            self.tau.as_slice(),
130            &mut work,
131            lwork,
132            &mut info,
133        );
134
135        q
136    }
137}
138
139/*
140 *
141 * Lapack functions dispatch.
142 *
143 */
144pub trait HessenbergScalar: Scalar + Copy {
145    fn xgehrd(
146        n: i32,
147        ilo: i32,
148        ihi: i32,
149        a: &mut [Self],
150        lda: i32,
151        tau: &mut [Self],
152        work: &mut [Self],
153        lwork: i32,
154        info: &mut i32,
155    );
156    fn xgehrd_work_size(
157        n: i32,
158        ilo: i32,
159        ihi: i32,
160        a: &mut [Self],
161        lda: i32,
162        tau: &mut [Self],
163        info: &mut i32,
164    ) -> i32;
165}
166
167/// Trait implemented by scalars for which Lapack implements the hessenberg decomposition.
168pub trait HessenbergReal: HessenbergScalar {
169    #[allow(missing_docs)]
170    fn xorghr(
171        n: i32,
172        ilo: i32,
173        ihi: i32,
174        a: &mut [Self],
175        lda: i32,
176        tau: &[Self],
177        work: &mut [Self],
178        lwork: i32,
179        info: &mut i32,
180    );
181    #[allow(missing_docs)]
182    fn xorghr_work_size(
183        n: i32,
184        ilo: i32,
185        ihi: i32,
186        a: &mut [Self],
187        lda: i32,
188        tau: &[Self],
189        info: &mut i32,
190    ) -> i32;
191}
192
193macro_rules! hessenberg_scalar_impl(
194    ($N: ty, $xgehrd: path) => (
195        impl HessenbergScalar for $N {
196            #[inline]
197            fn xgehrd(n: i32, ilo: i32, ihi: i32, a: &mut [Self], lda: i32,
198                      tau: &mut [Self], work: &mut [Self], lwork: i32, info: &mut i32) {
199                unsafe { $xgehrd(n, ilo, ihi, a, lda, tau, work, lwork, info) }
200            }
201
202            #[inline]
203            fn xgehrd_work_size(n: i32, ilo: i32, ihi: i32, a: &mut [Self], lda: i32,
204                                tau: &mut [Self], info: &mut i32) -> i32 {
205                let mut work = [ Zero::zero() ];
206                let lwork = -1 as i32;
207
208                unsafe { $xgehrd(n, ilo, ihi, a, lda, tau, &mut work, lwork, info) };
209                ComplexHelper::real_part(work[0]) as i32
210            }
211        }
212    )
213);
214
215macro_rules! hessenberg_real_impl(
216    ($N: ty, $xorghr: path) => (
217        impl HessenbergReal for $N {
218            #[inline]
219            fn xorghr(n: i32, ilo: i32, ihi: i32, a: &mut [Self], lda: i32, tau: &[Self],
220                      work: &mut [Self], lwork: i32, info: &mut i32) {
221                unsafe { $xorghr(n, ilo, ihi, a, lda, tau, work, lwork, info) }
222            }
223
224            #[inline]
225            fn xorghr_work_size(n: i32, ilo: i32, ihi: i32, a: &mut [Self], lda: i32,
226                                tau: &[Self], info: &mut i32) -> i32 {
227                let mut work = [ Zero::zero() ];
228                let lwork = -1 as i32;
229
230                unsafe { $xorghr(n, ilo, ihi, a, lda, tau, &mut work, lwork, info) };
231                ComplexHelper::real_part(work[0]) as i32
232            }
233        }
234    )
235);
236
237hessenberg_scalar_impl!(f32, lapack::sgehrd);
238hessenberg_scalar_impl!(f64, lapack::dgehrd);
239hessenberg_scalar_impl!(Complex<f32>, lapack::cgehrd);
240hessenberg_scalar_impl!(Complex<f64>, lapack::zgehrd);
241
242hessenberg_real_impl!(f32, lapack::sorghr);
243hessenberg_real_impl!(f64, lapack::dorghr);
nalgebra_lapack/hessenberg.rs

nalgebra_lapack/
hessenberg.rs
