use crate::allocator::Allocator;
use crate::base::{DefaultAllocator, OMatrix, OVector, Unit, Vector};
use crate::dimension::Dim;
use crate::storage::StorageMut;
use num::Zero;
use simba::scalar::ComplexField;
use crate::geometry::Reflection;
#[doc(hidden)]
#[inline(always)]
pub fn reflection_axis_mut<T: ComplexField, D: Dim, S: StorageMut<T, D>>(
column: &mut Vector<T, D, S>,
) -> (T, bool) {
let reflection_sq_norm = column.norm_squared();
let reflection_norm = reflection_sq_norm.clone().sqrt();
let factor;
let signed_norm;
unsafe {
let (modulus, sign) = column.vget_unchecked(0).clone().to_exp();
signed_norm = sign.scale(reflection_norm.clone());
factor = (reflection_sq_norm + modulus * reflection_norm) * crate::convert(2.0);
*column.vget_unchecked_mut(0) += signed_norm.clone();
};
if !factor.is_zero() {
column.unscale_mut(factor.sqrt());
let _ = column.normalize_mut();
(-signed_norm, true)
} else {
(signed_norm, false)
}
}
#[doc(hidden)]
#[must_use]
pub fn clear_column_unchecked<T: ComplexField, R: Dim, C: Dim>(
matrix: &mut OMatrix<T, R, C>,
icol: usize,
shift: usize,
bilateral: Option<&mut OVector<T, R>>,
) -> T
where
DefaultAllocator: Allocator<R, C> + Allocator<R>,
{
let (mut left, mut right) = matrix.columns_range_pair_mut(icol, icol + 1..);
let mut axis = left.rows_range_mut(icol + shift..);
let (reflection_norm, not_zero) = reflection_axis_mut(&mut axis);
if not_zero {
let refl = Reflection::new(Unit::new_unchecked(axis), T::zero());
let sign = reflection_norm.clone().signum();
if let Some(work) = bilateral {
refl.reflect_rows_with_sign(&mut right, work, sign.clone());
}
refl.reflect_with_sign(&mut right.rows_range_mut(icol + shift..), sign.conjugate());
}
reflection_norm
}
#[doc(hidden)]
#[must_use]
pub fn clear_row_unchecked<T: ComplexField, R: Dim, C: Dim>(
matrix: &mut OMatrix<T, R, C>,
axis_packed: &mut OVector<T, C>,
work: &mut OVector<T, R>,
irow: usize,
shift: usize,
) -> T
where
DefaultAllocator: Allocator<R, C> + Allocator<R> + Allocator<C>,
{
let (mut top, mut bottom) = matrix.rows_range_pair_mut(irow, irow + 1..);
let mut axis = axis_packed.rows_range_mut(irow + shift..);
axis.tr_copy_from(&top.columns_range(irow + shift..));
let (reflection_norm, not_zero) = reflection_axis_mut(&mut axis);
axis.conjugate_mut(); if not_zero {
let refl = Reflection::new(Unit::new_unchecked(axis), T::zero());
refl.reflect_rows_with_sign(
&mut bottom.columns_range_mut(irow + shift..),
&mut work.rows_range_mut(irow + 1..),
reflection_norm.clone().signum().conjugate(),
);
top.columns_range_mut(irow + shift..)
.tr_copy_from(refl.axis());
} else {
top.columns_range_mut(irow + shift..).tr_copy_from(&axis);
}
reflection_norm
}
#[doc(hidden)]
pub fn assemble_q<T: ComplexField, D: Dim>(m: &OMatrix<T, D, D>, signs: &[T]) -> OMatrix<T, D, D>
where
DefaultAllocator: Allocator<D, D>,
{
assert!(m.is_square());
let dim = m.shape_generic().0;
let mut res = OMatrix::identity_generic(dim, dim);
for i in (0..dim.value() - 1).rev() {
let axis = m.view_range(i + 1.., i);
let refl = Reflection::new(Unit::new_unchecked(axis), T::zero());
let mut res_rows = res.view_range_mut(i + 1.., i..);
refl.reflect_with_sign(&mut res_rows, signs[i].clone().signum());
}
res
}