diff --git a/src/linalg/bidiagonal.rs b/src/linalg/bidiagonal.rs index 3036e9f15..4faa9c5a0 100644 --- a/src/linalg/bidiagonal.rs +++ b/src/linalg/bidiagonal.rs @@ -8,6 +8,7 @@ use simba::scalar::ComplexField; use crate::geometry::Reflection; use crate::linalg::householder; +use crate::num::Zero; use std::mem::MaybeUninit; /// The bidiagonalization of a general matrix. @@ -227,7 +228,11 @@ where for i in (0..dim - shift).rev() { let axis = self.uv.view_range(i + shift.., i); - // TODO: sometimes, the axis might have a zero magnitude. + + // Sometimes, the axis might have a zero magnitude. + if axis.norm_squared().is_zero() { + continue; + } let refl = Reflection::new(Unit::new_unchecked(axis), T::zero()); let mut res_rows = res.view_range_mut(i + shift.., i..); @@ -263,7 +268,11 @@ where let axis = self.uv.view_range(i, i + shift..); let mut axis_packed = axis_packed.rows_range_mut(i + shift..); axis_packed.tr_copy_from(&axis); - // TODO: sometimes, the axis might have a zero magnitude. + + // Sometimes, the axis might have a zero magnitude. + if axis_packed.norm_squared().is_zero() { + continue; + } let refl = Reflection::new(Unit::new_unchecked(axis_packed), T::zero()); let mut res_rows = res.view_range_mut(i.., i + shift..); diff --git a/src/linalg/givens.rs b/src/linalg/givens.rs index 30e2bed38..dbd254391 100644 --- a/src/linalg/givens.rs +++ b/src/linalg/givens.rs @@ -17,7 +17,7 @@ pub struct GivensRotation { // Matrix = UnitComplex * Matrix impl GivensRotation { - /// The Givents rotation that does nothing. + /// The Givens rotation that does nothing. pub fn identity() -> Self { Self { c: T::RealField::one(), @@ -88,13 +88,13 @@ impl GivensRotation { } } - /// The cos part of this roration. + /// The cos part of this rotation. #[must_use] pub fn c(&self) -> T::RealField { self.c.clone() } - /// The sin part of this roration. + /// The sin part of this rotation. #[must_use] pub fn s(&self) -> T { self.s.clone() diff --git a/tests/linalg/bidiagonal.rs b/tests/linalg/bidiagonal.rs index aaee393fb..800967335 100644 --- a/tests/linalg/bidiagonal.rs +++ b/tests/linalg/bidiagonal.rs @@ -1,6 +1,6 @@ -#![cfg(feature = "proptest-support")] - -macro_rules! gen_tests( +#[cfg(feature = "proptest-support")] +mod proptest_tests { + macro_rules! gen_tests( ($module: ident, $scalar: expr) => { mod $module { #[allow(unused_imports)] @@ -54,8 +54,9 @@ macro_rules! gen_tests( } ); -gen_tests!(complex, complex_f64()); -gen_tests!(f64, PROPTEST_F64); + gen_tests!(complex, complex_f64()); + gen_tests!(f64, PROPTEST_F64); +} #[test] fn bidiagonal_identity() { @@ -74,3 +75,31 @@ fn bidiagonal_identity() { let (u, d, v_t) = bidiagonal.unpack(); assert_eq!(m, &u * d * &v_t); } + +#[test] +fn bidiagonal_regression_issue_1313() { + let s = 6.123234e-16_f32; + let mut m = nalgebra::dmatrix![ + 10.0, 0.0, 0.0, 0.0, -10.0, 0.0, 0.0, 0.0; + s, 10.0, 0.0, 10.0, s, 0.0, 0.0, 0.0; + 20.0, -20.0, 0.0, 20.0, 20.0, 0.0, 0.0, 0.0; + ]; + m.unscale_mut(m.camax()); + let bidiagonal = m.clone().bidiagonalize(); + let (u, d, v_t) = bidiagonal.unpack(); + let m2 = &u * d * &v_t; + assert_relative_eq!(m, m2, epsilon = 1e-6); +} + +#[test] +fn bidiagonal_regression_issue_1313_minimal() { + let s = 6.123234e-17_f32; + let m = nalgebra::dmatrix![ + 1.0, 0.0, -1.0; + s, 1.0, s; + ]; + let bidiagonal = m.clone().bidiagonalize(); + let (u, d, v_t) = bidiagonal.unpack(); + let m2 = &u * &d * &v_t; + assert_relative_eq!(m, m2, epsilon = 1e-6); +} diff --git a/tests/linalg/svd.rs b/tests/linalg/svd.rs index 900901ad3..d8b23d024 100644 --- a/tests/linalg/svd.rs +++ b/tests/linalg/svd.rs @@ -499,3 +499,17 @@ fn svd_regression_issue_1072() { epsilon = 1e-9 ); } + +#[test] +// Exercises bug reported in issue #1313 of nalgebra (https://github.com/dimforge/nalgebra/issues/1313) +fn svd_regression_issue_1313() { + let s = 6.123234e-16_f32; + let m = nalgebra::dmatrix![ + 10.0, 0.0, 0.0, 0.0, -10.0, 0.0, 0.0, 0.0; + s, 10.0, 0.0, 10.0, s, 0.0, 0.0, 0.0; + 20.0, -20.0, 0.0, 20.0, 20.0, 0.0, 0.0, 0.0; + ]; + let svd = m.clone().svd(true, true); + let m2 = svd.recompose().unwrap(); + assert_relative_eq!(&m, &m2, epsilon = 1e-5); +}