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
38        // Normalize again, making sure the vector is unit-sized.
39        // If `factor` had a very small value, the first normalization
40        // (dividing by `factor.sqrt()`) might end up with a slightly
41        // non-unit vector (especially when using 32-bits float).
42        // Decompositions strongly rely on that unit-vector property,
43        // so we run a second normalization (that is much more numerically
44        // stable since the norm is close to 1) to ensure it has a unit
45        // size.
46        let _ = column.normalize_mut();
47
48        (-signed_norm, true)
49    } else {
50        // TODO: not sure why we don't have a - sign here.
51        (signed_norm, false)
52    }
53}
54
55/// Uses an householder reflection to zero out the `icol`-th column, starting with the `shift + 1`-th
56/// subdiagonal element.
57///
58/// Returns the signed norm of the column.
59#[doc(hidden)]
60#[must_use]
61pub fn clear_column_unchecked<T: ComplexField, R: Dim, C: Dim>(
62    matrix: &mut OMatrix<T, R, C>,
63    icol: usize,
64    shift: usize,
65    bilateral: Option<&mut OVector<T, R>>,
66) -> T
67where
68    DefaultAllocator: Allocator<R, C> + Allocator<R>,
69{
70    let (mut left, mut right) = matrix.columns_range_pair_mut(icol, icol + 1..);
71    let mut axis = left.rows_range_mut(icol + shift..);
72
73    let (reflection_norm, not_zero) = reflection_axis_mut(&mut axis);
74
75    if not_zero {
76        let refl = Reflection::new(Unit::new_unchecked(axis), T::zero());
77        let sign = reflection_norm.clone().signum();
78        if let Some(work) = bilateral {
79            refl.reflect_rows_with_sign(&mut right, work, sign.clone());
80        }
81        refl.reflect_with_sign(&mut right.rows_range_mut(icol + shift..), sign.conjugate());
82    }
83
84    reflection_norm
85}
86
87/// Uses an householder reflection to zero out the `irow`-th row, ending before the `shift + 1`-th
88/// superdiagonal element.
89///
90/// Returns the signed norm of the column.
91#[doc(hidden)]
92#[must_use]
93pub fn clear_row_unchecked<T: ComplexField, R: Dim, C: Dim>(
94    matrix: &mut OMatrix<T, R, C>,
95    axis_packed: &mut OVector<T, C>,
96    work: &mut OVector<T, R>,
97    irow: usize,
98    shift: usize,
99) -> T
100where
101    DefaultAllocator: Allocator<R, C> + Allocator<R> + Allocator<C>,
102{
103    let (mut top, mut bottom) = matrix.rows_range_pair_mut(irow, irow + 1..);
104    let mut axis = axis_packed.rows_range_mut(irow + shift..);
105    axis.tr_copy_from(&top.columns_range(irow + shift..));
106
107    let (reflection_norm, not_zero) = reflection_axis_mut(&mut axis);
108    axis.conjugate_mut(); // So that reflect_rows actually cancels the first row.
109
110    if not_zero {
111        let refl = Reflection::new(Unit::new_unchecked(axis), T::zero());
112        refl.reflect_rows_with_sign(
113            &mut bottom.columns_range_mut(irow + shift..),
114            &mut work.rows_range_mut(irow + 1..),
115            reflection_norm.clone().signum().conjugate(),
116        );
117        top.columns_range_mut(irow + shift..)
118            .tr_copy_from(refl.axis());
119    } else {
120        top.columns_range_mut(irow + shift..).tr_copy_from(&axis);
121    }
122
123    reflection_norm
124}
125
126/// Computes the orthogonal transformation described by the elementary reflector axii stored on
127/// the lower-diagonal element of the given matrix.
128/// matrices.
129#[doc(hidden)]
130pub fn assemble_q<T: ComplexField, D: Dim>(m: &OMatrix<T, D, D>, signs: &[T]) -> OMatrix<T, D, D>
131where
132    DefaultAllocator: Allocator<D, D>,
133{
134    assert!(m.is_square());
135    let dim = m.shape_generic().0;
136
137    // NOTE: we could build the identity matrix and call p_mult on it.
138    // Instead we don't so that we take in account the matrix sparseness.
139    let mut res = OMatrix::identity_generic(dim, dim);
140
141    for i in (0..dim.value() - 1).rev() {
142        let axis = m.view_range(i + 1.., i);
143        let refl = Reflection::new(Unit::new_unchecked(axis), T::zero());
144
145        let mut res_rows = res.view_range_mut(i + 1.., i..);
146        refl.reflect_with_sign(&mut res_rows, signs[i].clone().signum());
147    }
148
149    res
150}
nalgebra/linalg/householder.rs

nalgebra/linalg/
householder.rs
