diff --git a/src/base/matrix_view.rs b/src/base/matrix_view.rs index e45d31f17..51be828c7 100644 --- a/src/base/matrix_view.rs +++ b/src/base/matrix_view.rs @@ -33,14 +33,6 @@ macro_rules! view_storage_impl ( #[deprecated = "Use ViewStorage(Mut) instead."] pub type $legacy_name<'a, T, R, C, RStride, CStride> = $T<'a, T, R, C, RStride, CStride>; - unsafe impl<'a, T: Send, R: Dim, C: Dim, RStride: Dim, CStride: Dim> Send - for $T<'a, T, R, C, RStride, CStride> - {} - - unsafe impl<'a, T: Sync, R: Dim, C: Dim, RStride: Dim, CStride: Dim> Sync - for $T<'a, T, R, C, RStride, CStride> - {} - impl<'a, T, R: Dim, C: Dim, RStride: Dim, CStride: Dim> $T<'a, T, R, C, RStride, CStride> { /// Create a new matrix view without bounds checking and from a raw pointer. /// @@ -136,6 +128,20 @@ impl Copy { } +/// Safety: Equivalent to a shared reference to `T`. All `Dim` type arguments are `Send + Sync`. A +/// shared reference can be sent iff `T: Sync`. +unsafe impl<'a, T: Sync, R: Dim, C: Dim, RStride: Dim, CStride: Dim> Send + for ViewStorage<'a, T, R, C, RStride, CStride> +{ +} + +/// Safety: Equivalent to a shared reference to `T`. All `Dim` type arguments are `Send + Sync`. A +/// shared reference is `Sync` iff `T: Sync`. +unsafe impl<'a, T: Sync, R: Dim, C: Dim, RStride: Dim, CStride: Dim> Sync + for ViewStorage<'a, T, R, C, RStride, CStride> +{ +} + impl Clone for ViewStorage<'_, T, R, C, RStride, CStride> { @@ -145,6 +151,20 @@ impl Clone } } +/// Safety: Equivalent to a unique reference to `T`. All `Dim` type arguments are `Send + Sync`. A +/// unique reference is `Send` iff `T: Send`. +unsafe impl<'a, T: Send, R: Dim, C: Dim, RStride: Dim, CStride: Dim> Send + for ViewStorageMut<'a, T, R, C, RStride, CStride> +{ +} + +/// Safety: Equivalent to a unique reference to `T`. All `Dim` type arguments are `Send + Sync`. A +/// unique reference is `Sync` iff `T: Sync`. +unsafe impl<'a, T: Sync, R: Dim, C: Dim, RStride: Dim, CStride: Dim> Sync + for ViewStorageMut<'a, T, R, C, RStride, CStride> +{ +} + impl<'a, T: Scalar, R: Dim, C: Dim, RStride: Dim, CStride: Dim> ViewStorageMut<'a, T, R, C, RStride, CStride> where diff --git a/src/linalg/givens.rs b/src/linalg/givens.rs index ee8286da8..bf59a20cf 100644 --- a/src/linalg/givens.rs +++ b/src/linalg/givens.rs @@ -58,10 +58,13 @@ impl GivensRotation { /// /// Returns `None` if no rotation is needed (i.e. if `v.y == 0`). Otherwise, this returns the norm /// of `v` and the rotation `r` such that `R * v = [ |v|, 0.0 ]^t` where `|v|` is the norm of `v`. - pub fn cancel_y>(v: &Vector) -> Option<(Self, T)> { - if !v[1].is_zero() { - let (mod0, sign0) = v[0].clone().to_exp(); - let denom = (mod0.clone() * mod0.clone() + v[1].clone().modulus_squared()).sqrt(); + pub fn cancel_y>( + v: &Vector, + eps: T::RealField, + ) -> Option<(Self, T)> { + let (mod0, sign0) = v[0].clone().to_exp(); + let denom = (mod0.clone() * mod0.clone() + v[1].clone().modulus_squared()).sqrt(); + if !v[1].is_zero() && denom > eps { let c = mod0 / denom.clone(); let s = -v[1].clone() / sign0.clone().scale(denom.clone()); let r = sign0.scale(denom); @@ -75,10 +78,13 @@ impl GivensRotation { /// /// Returns `None` if no rotation is needed (i.e. if `v.x == 0`). Otherwise, this returns the norm /// of `v` and the rotation `r` such that `R * v = [ 0.0, |v| ]^t` where `|v|` is the norm of `v`. - pub fn cancel_x>(v: &Vector) -> Option<(Self, T)> { - if !v[0].is_zero() { - let (mod1, sign1) = v[1].clone().to_exp(); - let denom = (mod1.clone() * mod1.clone() + v[0].clone().modulus_squared()).sqrt(); + pub fn cancel_x>( + v: &Vector, + eps: T::RealField, + ) -> Option<(Self, T)> { + let (mod1, sign1) = v[1].clone().to_exp(); + let denom = (mod1.clone() * mod1.clone() + v[0].clone().modulus_squared()).sqrt(); + if !v[0].is_zero() && denom > eps { let c = mod1 / denom.clone(); let s = (v[0].clone().conjugate() * sign1.clone()).unscale(denom.clone()); let r = sign1.scale(denom); diff --git a/src/linalg/svd.rs b/src/linalg/svd.rs index 0509dec00..858525742 100644 --- a/src/linalg/svd.rs +++ b/src/linalg/svd.rs @@ -159,6 +159,26 @@ where let mut diagonal = bi_matrix.diagonal(); let mut off_diagonal = bi_matrix.off_diagonal(); + // Compute the Frobenius norm of the bidiagonal matrix for relative convergence checks, + // so we can use a relative tolerance rather than absolute. + // This avoids stability issues when diagonal elements are very small but not exactly zero. + let mut norm_sqr = T::RealField::zero(); + for i in 0..diagonal.len() { + norm_sqr += diagonal[i].clone() * diagonal[i].clone(); + } + for i in 0..off_diagonal.len() { + norm_sqr += off_diagonal[i].clone() * off_diagonal[i].clone(); + } + let b_norm = norm_sqr.sqrt(); + // Use relative epsilon: eps * dim * ||B|| for zero checks. + // + // Scaling by dim accounts for roundoff accumulation in the Householder + // reflections during bidiagonalization (similar to LAPACK's N*EPS approach + // in DBDSQR). Without this, near-rank-deficient complex matrices can + // have bidiagonal elements that are "almost zero" but above the threshold, + // leading to numerically unstable iterations on a near-zero subproblem. + let eps_rel = eps.clone() * crate::convert::<_, T::RealField>(dim as f64) * b_norm; + let mut niter = 0; let (mut start, mut end) = Self::delimit_subproblem( &mut diagonal, @@ -168,6 +188,7 @@ where bi_matrix.is_upper_diagonal(), dim - 1, eps.clone(), + eps_rel.clone(), ); while end != start { @@ -214,7 +235,7 @@ where m12, ); - match GivensRotation::cancel_y(&vec) { + match GivensRotation::cancel_y(&vec, eps.clone()) { Some((rot1, norm1)) => { rot1.inverse() .rotate_rows(&mut subm.fixed_columns_mut::<2>(0)); @@ -227,8 +248,7 @@ where } let v = Vector2::new(subm[(0, 0)].clone(), subm[(1, 0)].clone()); - // TODO: does the case `v.y == 0` ever happen? - let (rot2, norm2) = GivensRotation::cancel_y(&v) + let (rot2, norm2) = GivensRotation::cancel_y(&v, eps.clone()) .unwrap_or((GivensRotation::identity(), subm[(0, 0)].clone())); rot2.rotate(&mut subm.fixed_columns_mut::<2>(1)); @@ -317,6 +337,7 @@ where bi_matrix.is_upper_diagonal(), end, eps.clone(), + eps_rel.clone(), ); start = sub.0; end = sub.1; @@ -372,6 +393,7 @@ where is_upper_diagonal: bool, end: usize, eps: T::RealField, + eps_rel: T::RealField, ) -> (usize, usize) { let mut n = end; @@ -381,9 +403,10 @@ where if off_diagonal[m].is_zero() || off_diagonal[m].clone().norm1() <= eps.clone() * (diagonal[n].clone().norm1() + diagonal[m].clone().norm1()) + || off_diagonal[m].clone().norm1() <= eps_rel { off_diagonal[m] = T::RealField::zero(); - } else if diagonal[m].clone().norm1() <= eps { + } else if diagonal[m].clone().norm1() <= eps_rel { diagonal[m] = T::RealField::zero(); Self::cancel_horizontal_off_diagonal_elt( diagonal, @@ -405,7 +428,7 @@ where m - 1, ); } - } else if diagonal[n].clone().norm1() <= eps { + } else if diagonal[n].clone().norm1() <= eps_rel { diagonal[n] = T::RealField::zero(); Self::cancel_vertical_off_diagonal_elt( diagonal, @@ -432,12 +455,11 @@ where if off_diagonal[m].clone().norm1() <= eps.clone() * (diagonal[new_start].clone().norm1() + diagonal[m].clone().norm1()) + || off_diagonal[m].clone().norm1() <= eps_rel { off_diagonal[m] = T::RealField::zero(); break; - } - // TODO: write a test that enters this case. - else if diagonal[m].clone().norm1() <= eps { + } else if diagonal[m].clone().norm1() <= eps_rel { diagonal[m] = T::RealField::zero(); Self::cancel_horizontal_off_diagonal_elt( diagonal, @@ -482,7 +504,10 @@ where off_diagonal[i] = T::RealField::zero(); for k in i..end { - match GivensRotation::cancel_x(&v) { + match GivensRotation::cancel_x( + &v, + T::RealField::default_epsilon() * crate::convert(1000.0), + ) { Some((rot, norm)) => { let rot = GivensRotation::new_unchecked(rot.c(), T::from_real(rot.s())); diagonal[k + 1] = norm; @@ -522,7 +547,10 @@ where off_diagonal[i] = T::RealField::zero(); for k in (0..i + 1).rev() { - match GivensRotation::cancel_y(&v) { + match GivensRotation::cancel_y( + &v, + T::RealField::default_epsilon() * crate::convert(1000.0), + ) { Some((rot, norm)) => { let rot = GivensRotation::new_unchecked(rot.c(), T::from_real(rot.s())); diagonal[k] = norm; @@ -846,9 +874,10 @@ where } } -// Explicit formulae inspired from the paper "Computing the Singular Values of 2-by-2 Complex -// Matrices", Sanzheng Qiao and Xiaohong Wang. -// http://www.cas.mcmaster.ca/sqrl/papers/sqrl5.pdf +// Ported from +// https://www.netlib.org/lapack/explore-html/d8/da7/group__lasv2_ga96f9f244300d82921950e2c393b4b20f.html#ga96f9f244300d82921950e2c393b4b20f +// with some clean-up and modification to match our slightly different contract. +// (We want the smaller eigenvalue in the top left.) fn compute_2x2_uptrig_svd( m11: T, m12: T, @@ -862,42 +891,120 @@ fn compute_2x2_uptrig_svd( ) { let two: T::RealField = crate::convert(2.0f64); let half: T::RealField = crate::convert(0.5f64); + let four: T::RealField = crate::convert(4.0f64); - let denom = (m11.clone() + m22.clone()).hypot(m12.clone()) - + (m11.clone() - m22.clone()).hypot(m12.clone()); - - // NOTE: v1 is the singular value that is the closest to m22. - // This prevents cancellation issues when constructing the vector `csv` below. If we chose - // otherwise, we would have v1 ~= m11 when m12 is small. This would cause catastrophic - // cancellation on `v1 * v1 - m11 * m11` below. - let mut v1 = m11.clone() * m22.clone() * two / denom.clone(); - let mut v2 = half * denom; - - let mut u = None; - let mut v_t = None; - - if compute_u || compute_v { - let (csv, sgn_v) = GivensRotation::new( - m11.clone() * m12.clone(), - v1.clone() * v1.clone() - m11.clone() * m11.clone(), - ); - v1 *= sgn_v.clone(); - v2 *= sgn_v; + #[derive(PartialEq)] + enum MatrixEntry { + F, + G, + H, + } - if compute_v { - v_t = Some(csv.clone()); + // Determine the largest entry in the matrix (m11=f, m12=g, m22=h), and + // ensure that |f| >= |h| by swapping if necessary (and correcting the + // result later). + let mut abs_f = m11.clone().abs(); + let abs_g = m12.clone().abs(); + let mut abs_h = m22.clone().abs(); + let swap_f_and_h = abs_h > abs_f; + let largest_entry = if swap_f_and_h { + if abs_g > abs_h { + MatrixEntry::G + } else { + MatrixEntry::H } + } else { + if abs_g > abs_f { + MatrixEntry::G + } else { + MatrixEntry::F + } + }; + let (mut f, g, mut h) = (m11.clone(), m12.clone(), m22.clone()); + if swap_f_and_h { + (f, h) = (h, f); + (abs_f, abs_h) = (abs_h, abs_f); + } - let cu = (m11.scale(csv.c()) + m12 * csv.s()) / v1.clone(); - let su = (m22 * csv.s()) / v1.clone(); - let (csu, sgn_u) = GivensRotation::new(cu, su); - v1 *= sgn_u.clone(); - v2 *= sgn_u; - - if compute_u { - u = Some(csu); + // Solve the problem. + let (mut v1, mut v2, mut su, mut cu, mut sv, mut cv) = if abs_g.is_zero() { + (abs_h, abs_f, -T::one(), T::zero(), T::one(), T::zero()) + } else { + if largest_entry == MatrixEntry::G && abs_f.clone() / abs_g.clone() < T::default_epsilon() { + ( + if abs_h > T::one() { + abs_f / (abs_g.clone() / abs_h) + } else { + (abs_f / abs_g.clone()) * abs_h + }, + abs_g, + T::one(), + h.clone() / g.clone(), + -f.clone() / g.clone(), + T::one(), + ) + } else { + let d = abs_f.clone() - abs_h.clone(); + let l = if d == abs_f { + T::one() + } else { + d.clone() / abs_f.clone() + }; + let m = g.clone() / f.clone(); + let mut t = two.clone() - l.clone(); + let mm = m.clone() * m.clone(); + let s = (t.clone() * t.clone() + mm.clone()).sqrt(); + let r = if l == T::zero() { + m.clone().abs() + } else { + (l.clone() * l.clone() + mm).sqrt() + }; + let a = half * (s.clone() + r.clone()); + t = if m == T::zero() { + if l == T::zero() { + two.clone() * f.clone().signum() * g.clone().signum() + } else { + g.clone() / (d * f.clone().signum()) + m.clone() / t + } + } else { + (m.clone() / (s + t) + m.clone() / (r + l)) * (T::one() + a.clone()) + }; + let b = (t.clone() * t.clone() + four).sqrt(); + let sv = -two / b.clone(); + let cv = t / b; + let su = (sv.clone() - cv.clone() * m) / a.clone(); + let cu = (h.clone() / f.clone()) * cv.clone() / a.clone(); + (abs_h / a.clone(), abs_f * a, su, cu, sv, cv) } + }; + + // Swap if necessary. + if swap_f_and_h { + (su, cu, sv, cv) = (cv, sv, cu, su); } - (u, Vector2::new(v1, v2), v_t) + // Set the signs of the eigenvalues. + let (f_sign, g_sign, h_sign) = (f.signum(), g.signum(), h.signum()); + let tsign = match largest_entry { + MatrixEntry::F => sv.clone().signum() * su.clone().signum() * f_sign.clone(), + MatrixEntry::G => -cv.clone().signum() * su.clone().signum() * g_sign, + MatrixEntry::H => cv.clone().signum() * cu.clone().signum() * h_sign.clone(), + }; + v1 *= tsign.clone() * f_sign * h_sign; + v2 *= tsign; + + // Return the result. + ( + if compute_u { + Some(GivensRotation::new_unchecked(cu.clone(), su.clone())) + } else { + None + }, + Vector2::new(v1, v2), + if compute_v { + Some(GivensRotation::new_unchecked(cv.clone(), sv.clone())) + } else { + None + }, + ) } diff --git a/src/linalg/symmetric_eigen.rs b/src/linalg/symmetric_eigen.rs index 0b4eceec4..7e4f5d0b8 100644 --- a/src/linalg/symmetric_eigen.rs +++ b/src/linalg/symmetric_eigen.rs @@ -152,7 +152,7 @@ where for i in start..n { let j = i + 1; - match GivensRotation::cancel_y(&vec) { + match GivensRotation::cancel_y(&vec, T::RealField::zero()) { Some((rot, norm)) => { if i > start { // Not the first iteration. @@ -205,10 +205,21 @@ where diag[start + 1].clone(), ); let eigvals = m.eigenvalues().unwrap(); - let basis = Vector2::new( - eigvals.x.clone() - diag[start + 1].clone(), - off_diag[start].clone(), - ); + + // Choose the basis least likely to experience cancellation + let basis = if (eigvals.x.clone() - diag[start + 1].clone()).abs() + > (eigvals.x.clone() - diag[start].clone()).abs() + { + Vector2::new( + eigvals.x.clone() - diag[start + 1].clone(), + off_diag[start].clone(), + ) + } else { + Vector2::new( + off_diag[start].clone(), + eigvals.x.clone() - diag[start].clone(), + ) + }; diag[start] = eigvals[0].clone(); diag[start + 1] = eigvals[1].clone(); @@ -348,7 +359,36 @@ where #[cfg(test)] mod test { - use crate::base::Matrix2; + use crate::base::{Matrix2, Matrix4}; + + /// Exercises bug reported in issue #1109 of nalgebra (https://github.com/dimforge/nalgebra/issues/1109) + #[test] + fn symmetric_eigen_regression_issue_1109() { + let m = Matrix4::new( + -19884.07f64, + -10.07188, + 11.277279, + -188560.63, + -10.07188, + 12.518197, + 1.3770627, + -102.97504, + 11.277279, + 1.3770627, + 14.587362, + 113.26099, + -188560.63, + -102.97504, + 113.26099, + -1788112.3, + ); + let eig = m.symmetric_eigen(); + assert!(relative_eq!( + m.lower_triangle(), + eig.recompose().lower_triangle(), + epsilon = 1.0e-5 + )); + } fn expected_shift(m: Matrix2) -> f64 { let vals = m.eigenvalues().unwrap(); diff --git a/tests/linalg/svd.rs b/tests/linalg/svd.rs index d8b23d024..e082d48bf 100644 --- a/tests/linalg/svd.rs +++ b/tests/linalg/svd.rs @@ -513,3 +513,178 @@ fn svd_regression_issue_1313() { let m2 = svd.recompose().unwrap(); assert_relative_eq!(&m, &m2, epsilon = 1e-5); } + +#[test] +// Accuracy bug reported in issue #1172 of nalgebra (https://github.com/dimforge/nalgebra/issues/1172) +fn svd_regression_issue_1172() { + use nalgebra::{Complex, Matrix4}; + type M4C = Matrix4>; + + let m = M4C::new( + Complex::new(0.4846888711394364, -0.0000000000000002529226450498658), + Complex::new(0.4997655143494952, -0.00000000000000001731471891552503), + Complex::new(0.00000000000000001527512211317369, 0.49976551434949495), + Complex::new(0.00000000000000009643636194372324, -0.48468887113943676), + Complex::new(0.4997655143494954, -0.00000000000000023999992468308935), + Complex::new(0.5153111288605629, -0.0000000000000002516108359162637), + Complex::new(-0.000000000000000005044961400919933, 0.5153111288605631), + Complex::new(0.000000000000000042486770760464153, -0.49976551434949495), + Complex::new(-0.000000000000000032383811520876863, -0.49976551434949484), + Complex::new(0.000000000000000014710206963221994, -0.5153111288605633), + Complex::new(0.5153111288605635, 0.00000000000000016459157448729957), + Complex::new(-0.4997655143494946, -0.00000000000000028045339796436483), + Complex::new(0.00000000000000007455816688442185, 0.4846888711394367), + Complex::new(0.0000000000000000641248522548626, 0.49976551434949484), + Complex::new(-0.49976551434949473, -0.00000000000000018487105599880945), + Complex::new(0.48468887113943704, 0.00000000000000015953066497724206), + ); + let m_singular_values = nalgebra::dvector![2.0, 0.0, 0.0, 0.0]; + + let svd = m.svd(true, true); + let sings = svd.singular_values; + let u = svd.u.unwrap(); + let v_t = svd.v_t.unwrap(); + let sigma = M4C::from_diagonal(&sings.cast::>()); + let m1 = u * sigma * v_t; + + // Should be accurate to machine precision + assert_relative_eq!(m, m1, epsilon = 1e-12); + + for (s, expected) in sings.iter().zip(m_singular_values.iter()) { + assert!( + (*s - expected).abs() < 1e-12, + "Singular value {s:e} is not accurate: expected {expected:e}", + ); + } +} + +#[test] +// Accuracy bug reported as feedback to PR #1590 of nalgebra (https://github.com/dimforge/nalgebra/issues/1590) +fn svd_regression_platform_dependent_accuracy() { + use nalgebra::{Complex, Matrix4}; + type M4C = Matrix4>; + + let m = M4C::new( + Complex::new(0.5163130224597328, 0.2640110414676673), + Complex::new(-0.10845827476820835, 0.34220148642893244), + Complex::new(0.14618038920991627, 0.3278663576677311), + Complex::new(0.4834191243671928, -0.32029192524071315), + Complex::new(0.25011023587834597, -0.5693169162970136), + Complex::new(0.3731433798035375, 0.09455746917888716), + Complex::new(0.34176716325705053, -0.17712228258452342), + Complex::new(-0.3732958190779979, -0.49731992995651203), + Complex::new(0.3732958190779976, -0.49731992995651314), + Complex::new(0.34176716325704987, 0.17712228258452062), + Complex::new(0.37314337980353635, -0.09455746917888491), + Complex::new(-0.25011023587834524, -0.5693169162970155), + Complex::new(0.4834191243671925, 0.32029192524071587), + Complex::new(-0.14618038920991572, 0.32786635766773), + Complex::new(0.10845827476820777, 0.3422014864289315), + Complex::new(0.5163130224597317, -0.26401104146766996), + ); + let m_singular_values = nalgebra::dvector![2.0, 0.0, 0.0, 0.0]; + + let svd = m.svd(true, true); + let sings = svd.singular_values; + let u = svd.u.unwrap(); + let v_t = svd.v_t.unwrap(); + let sigma = M4C::from_diagonal(&sings.cast::>()); + let m1 = u * sigma * v_t; + + // Should be accurate to machine precision + assert_relative_eq!(m, m1, epsilon = 1e-12); + + for (s, expected) in sings.iter().zip(m_singular_values.iter()) { + assert!( + (*s - expected).abs() < 1e-12, + "Singular value {s:e} is not accurate: expected {expected:e}", + ); + } +} + +#[test] +// Accuracy bug reported as feedback to PR #1590 of nalgebra (https://github.com/dimforge/nalgebra/issues/1590) +fn svd_regression_1() { + use nalgebra::{Complex, Matrix4}; + type M4C = Matrix4>; + + let m = M4C::new( + Complex::new(0.5507428877419177, -0.22418277708063508), + Complex::new(0.3304118701150042, -0.19300867894787538), + Complex::new(0.38115865977169494, -0.033799854188208585), + Complex::new(-0.5788888093103887, 0.13588006620948034), + Complex::new(-0.5788888093104021, -0.13588006620946957), + Complex::new(-0.3811586597717028, -0.03379985418820041), + Complex::new(-0.33041187011500606, -0.19300867894788695), + Complex::new(0.5507428877419225, 0.22418277708065218), + Complex::new(-0.5507428877419267, 0.22418277708064818), + Complex::new(-0.3304118701150135, 0.1930086789478659), + Complex::new(-0.38115865977169644, 0.03379985418822192), + Complex::new(0.5788888093104035, -0.13588006620947532), + Complex::new(-0.5788888093103961, -0.1358800662094592), + Complex::new(-0.38115865977169044, -0.0337998541882125), + Complex::new(-0.33041187011500334, -0.19300867894787002), + Complex::new(0.5507428877419116, 0.22418277708065648), + ); + let m_singular_values = nalgebra::dvector![2.0, 0.0, 0.0, 0.0]; + + let svd = m.svd(true, true); + let sings = svd.singular_values; + let u = svd.u.unwrap(); + let v_t = svd.v_t.unwrap(); + let sigma = M4C::from_diagonal(&sings.cast::>()); + let m1 = u * sigma * v_t; + + // Should be accurate to machine precision + assert_relative_eq!(m, m1, epsilon = 1e-12); + + for (s, expected) in sings.iter().zip(m_singular_values.iter()) { + assert!( + (*s - expected).abs() < 1e-12, + "Singular value {s:e} is not accurate: expected {expected:e}", + ); + } +} + +#[test] +fn svd_regression_2() { + use nalgebra::{Complex, Matrix4}; + type M4C = Matrix4>; + + let m = M4C::new( + Complex::new(-0.4999999999999969, 0.00000000000000009847209127792333), + Complex::new(-0.5000000000000002, -0.0000000000000005145817176957074), + Complex::new(-0.5000000000000031, -0.00000000000000011326570795854838), + Complex::new(0.4999999999999997, 0.00000000000000028317808990753396), + Complex::new(-0.49999999999999994, 0.0000000000000002324728614724245), + Complex::new(-0.5000000000000029, 0.00000000000000009174755495786714), + Complex::new(-0.49999999999999994, 0.00000000000000037035002208515673), + Complex::new(0.4999999999999969, -0.0000000000000000779181294373365), + Complex::new(-0.5000000000000031, -0.00000000000000010725752877840523), + Complex::new(-0.4999999999999997, -0.0000000000000004179869031467993), + Complex::new(-0.499999999999997, -0.000000000000000027676985470314004), + Complex::new(0.5000000000000002, 0.0000000000000004129804445544165), + Complex::new(0.5000000000000001, -0.000000000000000490409530759043), + Complex::new(0.4999999999999969, -0.00000000000000008015506216445231), + Complex::new(0.5000000000000001, -0.00000000000000032732792387061645), + Complex::new(-0.5000000000000029, 0.00000000000000010807420148690568), + ); + let m_singular_values = nalgebra::dvector![2.0, 0.0, 0.0, 0.0]; + + let svd = m.svd(true, true); + let sings = svd.singular_values; + let u = svd.u.unwrap(); + let v_t = svd.v_t.unwrap(); + let sigma = M4C::from_diagonal(&sings.cast::>()); + let m1 = u * sigma * v_t; + + // Should be accurate to machine precision + assert_relative_eq!(m, m1, epsilon = 1e-12); + + for (s, expected) in sings.iter().zip(m_singular_values.iter()) { + assert!( + (*s - expected).abs() < 1e-12, + "Singular value {s:e} is not accurate: expected {expected:e}", + ); + } +}