use super::*;
use crate::assert;
fn issue81_matrix<T: ComplexField>(make: fn(f64, f64) -> T) -> Mat<T> {
let entries = [
(0, 0, 2.3594724911802776, 0.006177867158102503),
(0, 1, 0.14611279140110442, -0.006177867158102429),
(0, 6, 0.360966787092662, 3.1078018129467733),
(0, 7, 0.01751975679899344, 0.19373242661916215),
(1, 0, -0.0679103720800715, 0.22101999154443816),
(1, 1, -0.009059753243836796, -0.2210199915444376),
(1, 6, -0.28638655061977014, 0.08432691582930801),
(1, 7, 0.29309227638075863, -0.025832938808607833),
(2, 0, 0.01575108897601647, -0.0064054967869236245),
(2, 1, -0.16593049000885574, 0.006405496786924529),
(2, 6, -0.004996754256389675, 0.010278618642714205),
(2, 7, -0.019060840267077314, -0.2201327688699224),
(3, 0, -0.07767698796918937, -0.006873965170532744),
(3, 1, -0.005684608987456119, 0.006873965170532807),
(3, 6, 0.015651423951280847, 0.05736152331917376),
(3, 7, -0.008821719809322738, 0.0022139158468607797),
(4, 0, 7.235274104261472e-08, -0.007137482256639182),
(4, 1, 0.0003427142346573653, 0.0071374822566606395),
(4, 6, -0.005327872093894844, 0.0006106378793937996),
(4, 7, 0.005379612470294339, -0.0001593071066018335),
(5, 0, 0.0003118896463960938, -4.038096293900582e-05),
(5, 1, -0.0055536553672519175, 4.038096293906832e-05),
(5, 6, 2.6332376524858047e-05, -0.0002404425657241616),
(5, 7, 0.00042327870925504816, 0.004162395454531604),
(6, 0, -3.277976868456115e-17, -1.385977603087669e-18),
(6, 2, 2.3594724911802785, -0.006177867158102503),
(6, 3, -0.14611279140110334, -0.0061778671581025815),
(6, 6, -6.0073910004316986e-18, -4.322493044612878e-17),
(6, 8, 0.35170915334452946, 3.108863104318509),
(6, 9, -0.02677739054712593, -0.19267113524742585),
(7, 0, 2.0325162008358368e-18, -4.958487294655612e-17),
(7, 2, -0.06791037208007403, -0.22101999154443824),
(7, 3, 0.009059753243834264, -0.22101999154443888),
(7, 6, 6.535963938086967e-17, -9.235304269143807e-18),
(7, 8, 0.29804028174295966, 0.01732846120549322),
(7, 9, 0.2913345559819713, -0.04116551581520697),
(8, 0, 3.722578309570986e-17, 1.437045319383591e-18),
(8, 2, 0.015751088976016386, 0.006405496786923625),
(8, 3, 0.16593049000885565, 0.006405496786922724),
(8, 6, 7.010985095238214e-18, 4.9072315645412575e-17),
(8, 8, 0.00719322401019955, 0.008881164486720506),
(8, 9, 0.03125081853366654, 0.21873531471392868),
(9, 0, 1.2753172797818566e-18, 1.5421441618840358e-18),
(9, 2, -0.07767698796918937, 0.006873965170532747),
(9, 3, 0.005684608987456119, 0.006873965170532689),
(9, 6, -2.0401709699413828e-18, -3.5913136499417674e-20),
(9, 8, -0.002264169528176444, 0.05941535967699379),
(9, 9, -0.009093873670134558, -0.00016007948904075715),
(10, 0, -7.688644682402164e-20, 1.6012630729900756e-18),
(10, 2, 7.235184428662789e-08, 0.007137482256639182),
(10, 3, -0.0003427142355541213, 0.007137482256617727),
(10, 6, 1.1836697084111227e-18, -2.383132402534741e-19),
(10, 8, 0.005327838778958875, -0.0006109284852386418),
(10, 9, 0.005276098402559383, -0.001062259258030608),
(11, 0, 1.2459384051635648e-18, 9.059293246712989e-21),
(11, 2, 0.00031188964639603114, 4.038096293900582e-05),
(11, 3, 0.005553655367251855, 4.0380962938943335e-05),
(11, 6, -1.1882998220841186e-19, -9.310789395413496e-19),
(11, 8, -8.006258488245128e-05, -0.0002282454905740448),
(11, 9, -0.0005296736706623575, -0.004150198379381488),
];
let mut a = Mat::<T>::zeros(12, 10);
for (i, j, re, im) in entries {
a[(i, j)] = make(re, im);
}
a
}
#[test]
fn test_issue81_complex_qr_stall() {
fn check<T: ComplexField>(a: Mat<T>) {
let expected = [
3.9414676652919867,
3.9414676652919858,
0.48150883185813453,
0.48150883185813415,
0.11330659717393496,
0.11330659717393492,
0.011919830114133985,
0.011919830114133898,
2.8733257864429986e-17,
7.8282169309112482e-18,
];
for scale in [1.0, 2.0, 1e-8, 1e8] {
let a = Mat::<T>::from_fn(12, 10, |i, j| {
a[(i, j)].mul_real(from_f64::<T::Real>(scale))
});
for adjoint in [false, true] {
let a = if adjoint {
a.adjoint().to_owned()
} else {
a.clone()
};
let (m, n) = a.shape();
for threshold in [128, usize::MAX] {
let mut params = Spec::<SvdParams, T>::default();
params.recursion_threshold = threshold;
for (want_u, want_v) in [
(true, true),
(true, false),
(false, true),
(false, false),
] {
let vectors = |want| {
if want {
ComputeSvdVectors::Thin
} else {
ComputeSvdVectors::No
}
};
let mut s = Col::<T>::zeros(10);
let mut u = Mat::<T>::zeros(m, 10);
let mut v = Mat::<T>::zeros(n, 10);
let mut mem = MemBuffer::new(svd_scratch::<T>(
m,
n,
vectors(want_u),
vectors(want_v),
ParShape::Seq,
params,
));
svd(
a.as_ref(),
s.as_mut().as_diagonal_mut(),
if want_u { Some(u.as_mut()) } else { None },
if want_v { Some(v.as_mut()) } else { None },
Par::Seq,
MemStack::new(&mut mem),
params,
)
.unwrap();
assert!(s.is_all_finite());
let tol = from_f64::<T::Real>(64.0 * m.max(n) as f64)
* eps::<T::Real>();
for i in 0..10 {
assert_eq!(s[i].imag(), zero::<T::Real>());
assert!(s[i].real() >= zero::<T::Real>());
if i > 0 {
assert!(s[i - 1].real() >= s[i].real());
}
let value =
s[i].real() / from_f64::<T::Real>(scale);
assert!(
(value - from_f64::<T::Real>(expected[i]))
.abs() <= &tol
* from_f64::<T::Real>(expected[0])
);
}
if want_u {
assert!(
(u.adjoint() * &u - Mat::<T>::identity(10, 10))
.norm_max() <= tol
);
}
if want_v {
assert!(
(v.adjoint() * &v - Mat::<T>::identity(10, 10))
.norm_max() <= tol
);
}
if want_u && want_v {
let back = &u * s.as_diagonal() * v.adjoint();
assert!(
(&back - &a).norm_l2() / a.norm_l2() <= tol
);
}
}
}
}
}
}
check(issue81_matrix::<c64>(c64::new));
check(issue81_matrix::<c32>(|re, im| {
c32::new(re as f32, im as f32)
}));
}