nalgebra/linalg/
householder.rs

1//! Construction of householder elementary reflections.
2
3use crate::allocator::Allocator;
4use crate::base::{DefaultAllocator, OMatrix, OVector, Unit, Vector};
5use crate::dimension::Dim;
6use crate::storage::StorageMut;
7use num::Zero;
8use simba::scalar::ComplexField;
9
10use crate::geometry::Reflection;
11
12/// Replaces `column` by the axis of the householder reflection that transforms `column` into
13/// `(+/-|column|, 0, ..., 0)`.
14///
15/// The unit-length axis is output to `column`. Returns what would be the first component of
16/// `column` after reflection and `false` if no reflection was necessary.
17#[doc(hidden)]
18#[inline(always)]
19pub fn reflection_axis_mut<T: ComplexField, D: Dim, S: StorageMut<T, D>>(
20    column: &mut Vector<T, D, S>,
21) -> (T, bool) {
22    let reflection_sq_norm = column.norm_squared();
23    let reflection_norm = reflection_sq_norm.clone().sqrt();
24
25    let factor;
26    let signed_norm;
27
28    unsafe {
29        let (modulus, sign) = column.vget_unchecked(0).clone().to_exp();
30        signed_norm = sign.scale(reflection_norm.clone());
31        factor = (reflection_sq_norm + modulus * reflection_norm) * crate::convert(2.0);
32        *column.vget_unchecked_mut(0) += signed_norm.clone();
33    };
34
35    if !factor.is_zero() {
36        column.unscale_mut(factor.sqrt());
37        (-signed_norm, true)
38    } else {
39        // TODO: not sure why we don't have a - sign here.
40        (signed_norm, false)
41    }
42}
43
44/// Uses an householder reflection to zero out the `icol`-th column, starting with the `shift + 1`-th
45/// subdiagonal element.
46///
47/// Returns the signed norm of the column.
48#[doc(hidden)]
49#[must_use]
50pub fn clear_column_unchecked<T: ComplexField, R: Dim, C: Dim>(
51    matrix: &mut OMatrix<T, R, C>,
52    icol: usize,
53    shift: usize,
54    bilateral: Option<&mut OVector<T, R>>,
55) -> T
56where
57    DefaultAllocator: Allocator<T, R, C> + Allocator<T, R>,
58{
59    let (mut left, mut right) = matrix.columns_range_pair_mut(icol, icol + 1..);
60    let mut axis = left.rows_range_mut(icol + shift..);
61
62    let (reflection_norm, not_zero) = reflection_axis_mut(&mut axis);
63
64    if not_zero {
65        let refl = Reflection::new(Unit::new_unchecked(axis), T::zero());
66        let sign = reflection_norm.clone().signum();
67        if let Some(mut work) = bilateral {
68            refl.reflect_rows_with_sign(&mut right, &mut work, sign.clone());
69        }
70        refl.reflect_with_sign(&mut right.rows_range_mut(icol + shift..), sign.conjugate());
71    }
72
73    reflection_norm
74}
75
76/// Uses an householder reflection to zero out the `irow`-th row, ending before the `shift + 1`-th
77/// superdiagonal element.
78///
79/// Returns the signed norm of the column.
80#[doc(hidden)]
81#[must_use]
82pub fn clear_row_unchecked<T: ComplexField, R: Dim, C: Dim>(
83    matrix: &mut OMatrix<T, R, C>,
84    axis_packed: &mut OVector<T, C>,
85    work: &mut OVector<T, R>,
86    irow: usize,
87    shift: usize,
88) -> T
89where
90    DefaultAllocator: Allocator<T, R, C> + Allocator<T, R> + Allocator<T, C>,
91{
92    let (mut top, mut bottom) = matrix.rows_range_pair_mut(irow, irow + 1..);
93    let mut axis = axis_packed.rows_range_mut(irow + shift..);
94    axis.tr_copy_from(&top.columns_range(irow + shift..));
95
96    let (reflection_norm, not_zero) = reflection_axis_mut(&mut axis);
97    axis.conjugate_mut(); // So that reflect_rows actually cancels the first row.
98
99    if not_zero {
100        let refl = Reflection::new(Unit::new_unchecked(axis), T::zero());
101        refl.reflect_rows_with_sign(
102            &mut bottom.columns_range_mut(irow + shift..),
103            &mut work.rows_range_mut(irow + 1..),
104            reflection_norm.clone().signum().conjugate(),
105        );
106        top.columns_range_mut(irow + shift..)
107            .tr_copy_from(refl.axis());
108    } else {
109        top.columns_range_mut(irow + shift..).tr_copy_from(&axis);
110    }
111
112    reflection_norm
113}
114
115/// Computes the orthogonal transformation described by the elementary reflector axii stored on
116/// the lower-diagonal element of the given matrix.
117/// matrices.
118#[doc(hidden)]
119pub fn assemble_q<T: ComplexField, D: Dim>(m: &OMatrix<T, D, D>, signs: &[T]) -> OMatrix<T, D, D>
120where
121    DefaultAllocator: Allocator<T, D, D>,
122{
123    assert!(m.is_square());
124    let dim = m.shape_generic().0;
125
126    // NOTE: we could build the identity matrix and call p_mult on it.
127    // Instead we don't so that we take in account the matrix sparseness.
128    let mut res = OMatrix::identity_generic(dim, dim);
129
130    for i in (0..dim.value() - 1).rev() {
131        let axis = m.slice_range(i + 1.., i);
132        let refl = Reflection::new(Unit::new_unchecked(axis), T::zero());
133
134        let mut res_rows = res.slice_range_mut(i + 1.., i..);
135        refl.reflect_with_sign(&mut res_rows, signs[i].clone().signum());
136    }
137
138    res
139}
nalgebra/linalg/householder.rs

nalgebra/linalg/
householder.rs
