From ae6879769366cce4fc5c721227cfd59a7938f1bc Mon Sep 17 00:00:00 2001 From: Hiroshi Shinaoka Date: Mon, 9 Nov 2026 13:53:26 +0900 Subject: [PATCH 1/5] fix(gevd): rotate the columns of Q when the complex QZ chases a zero down B (#324) hessenberg_to_qz_unblocked (the complex QZ path for n < 15) chases a negligible T[(j, j)] at j == ilo, or below a negligible H[(j, j - 1)], to the bottom with rotations on rows jch, jch + 1 of H and T. Q has to take the adjoint rotation on its columns, as the QZ sweep, the other zero-chase branch and the blocked path already do (LAPACK zhgeqz rotates Q(1,JCH) and Q(1,JCH+1) with C, DCONJG(S)). It rotated rows jch, jch + 1 of Q with s. Q stayed unitary but A = Q S Z^H no longer held, so every left eigenvector of gevd_cplx was wrong, while the eigenvalues and the right eigenvectors, which only use S, T and Z, were right. Whether a pencil reaches the branch depends on rounding: the issue's 12 x 12 c32 pencil does on x86-64 with AVX2/FMA (left residual 1.962e-1 -> 2.079e-7) and does not on aarch64. Tests: a unit test hands the unblocked iteration an exact zero on B's diagonal (j == ilo and below a zero subdiagonal, c32/c64) and checks the reconstruction; an integration test embeds the issue's bit patterns and checks left and right residuals at c32 and c64. The same lines are in upstream faer (Codeberg main 5833f4e). Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 5 + PATCHSET.md | 1 + faer/src/linalg/gevd/qz_cplx/mod.rs | 102 ++++++++++- faer/tests/gevd_infinite_eigenvalue.rs | 226 +++++++++++++++++++++++++ 4 files changed, 329 insertions(+), 5 deletions(-) create mode 100644 faer/tests/gevd_infinite_eigenvalue.rs diff --git a/CHANGELOG.md b/CHANGELOG.md index 8c2c26c8..9a1f539a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,11 @@ Changes since `t4a-faer` 0.32.0 (tag `faer-v0.32.0`). Numbers are tensor4all/fae ## Correctness +- `gevd_cplx` returned wrong left eigenvectors for some pencils with an infinite eigenvalue (a zero on the diagonal of `B` during the QZ + iteration, `n < 15`), while the eigenvalues and the right eigenvectors were right. The rotations that chase that zero to the bottom were + applied to the rows of `Q` instead of its columns. Whether a pencil reaches that branch depends on rounding, so the same input could fail on + x86-64 with AVX2 and pass on aarch64 (#324). + - `stats::col_mean` and `stats::col_varm` wrote their results in reverse order when the input had a negative row stride (for example `mat.reverse_rows()`), and `row_mean`/`row_varm` did the same for a negative column stride (`mat.reverse_cols()`). Both `NanHandling` modes were affected, and `varm` also subtracted the mean of the mirrored row. The input rows were reversed to get a positive stride, but `out` and the diff --git a/PATCHSET.md b/PATCHSET.md index 367f7309..5a1529e3 100644 --- a/PATCHSET.md +++ b/PATCHSET.md @@ -193,6 +193,7 @@ gate has changed and two kinds of adoption are now in the ledger. | `norm_max` propagates NaN (#316) | `linalg::reductions::norm_max` (and `MatRef`/`ColRef`/`RowRef::norm_max`) skipped NaN entries: the SIMD kernel accumulates with `simd.max`, a `select(acc > x, acc, x)` that keeps the non-NaN operand, and the pairwise halves and the scalar loop combine with `fmax` (IEEE `maxNum`). It returned the largest finite magnitude, and 0 for an all-NaN input, while `norm_l1`/`norm_l2` returned NaN. It now returns NaN when any entry, or either part of a complex entry, is NaN, as LAPACK's `xLANGE('M')` and Julia's `norm(x, Inf)` do. The SIMD kernel keeps its max accumulators and adds a second set that sums the magnitudes: every term is `>= 0`, `+inf` or NaN, so the sum is NaN exactly when an entry is NaN, and overflow gives `inf`, not NaN. The pairwise combine propagates NaN, and the scalar loop returns NaN at the first NaN entry. An empty matrix still returns 0. This changes behaviour: internal callers that take `norm_max` of a matrix holding NaN (QZ, tridiagonal EVD, bidiagonal SVD, Bunch-Kaufman scaling) now see NaN instead of a finite norm. Both files are byte-identical to the upstream sync base `0ecd846`, so upstream faer carries the same behaviour there; whether later upstream changed it is not checked, and feedback is deferred. | EXPERIMENT | tensor4all/faer#316 | Correctness only, no timing claim. The SIMD inner loop does one extra vector add per batch. New tests in `linalg::reductions::norm_max::tests`: f32/f64/c32/c64 with NaN in a real entry and in the real or imaginary part of a complex entry, at the first, middle and last position, and all-NaN, across shapes on both sides of the SIMD tail and the pairwise split. Each case goes through the contiguous (SIMD) path and a row-stride-2 view (scalar loop), and through the `MatRef`/`ColRef`/`RowRef` methods. Finite and `inf` inputs give the same result on both paths, and `norm_l1`/`norm_l2` agree on NaN. Removing the SIMD sum line or the scalar NaN return fails `test_norm_max_propagates_nan`. | | `stats` column/row reductions with a reversed output axis (#317) | `col_mean`/`col_varm`, and through them `row_mean`/`row_varm`, normalize a negative row stride with `mat.reverse_rows()` before the kernel runs. Rows are the output axis of a column reduction, and `out` (and, for `varm`, `col_mean`) were not reversed with them. Entry `i` was written to `nrows - 1 - i`, and `varm` subtracted the mean of the mirrored row. All four internal functions (`col_{mean,varm}_{ignore,propagate}`) now reverse `out`, and `col_mean` for `varm`, together with `mat`. Reversing the reduced axis (`reverse_cols` in `col_*`) was already harmless. Public API unchanged. `meanvar.rs` is byte-identical to the upstream sync base `0ecd846`, so upstream faer carries the same defect there; later upstream is not checked, and feedback is deferred. | EXPERIMENT | tensor4all/faer#317 | Correctness only, no timing claim. New tests in `stats::meanvar::tests`: `col_*` on `reverse_rows()` and `row_*` on `reverse_cols()`, mean and varm, `Propagate` and `Ignore`, f32/f64/c32/c64, on column-major (fallback kernel) and row-major (row-major SIMD kernel) storage with one NaN. Each is compared with the same matrix copied to an owned one. The issue's reproducer is pinned as well. Removing the output reversal from any one of the four functions fails both tests. | | Complex `abs_max` drops a NaN in the real part (#320) | `ComplexField::simd_abs_max` takes the larger of the two part magnitudes with an unordered `re > im` select, and that select returns `im` when `re` is NaN: a NaN real part read as the finite imaginary magnitude, while a NaN imaginary part did return NaN. The imaginary side already lands on `im`, so the real side takes an explicit `!(re == re)` mask in the generic `Complex` SIMD path, and the scalar fallbacks of `ComplexImpl` / `ComplexImpl` test both parts before the comparison. No in-tree type reaches either path (native `c32` / `c64` go through the real coercion in `norm_max`), so no dispatch, allocation or numerical result on an existing input changes. | KEEP | tensor4all/faer#320 | Correctness; the `faer-traits` unit test `complex_abs_max_propagates_nan` drives `SimdCtx::abs_max` on `Complex` and fails before / passes after (real-part NaN, imaginary-part NaN, both, and three finite controls). No performance claim: two mask ops per SIMD chunk on a path no in-tree type reaches. | +| Complex QZ: a zero chased down `B`'s diagonal updates the columns of `Q` (#324) | When `hessenberg_to_qz_unblocked` (the complex QZ path for `n < 15`) finds `T[(j, j)]` negligible with `j == ilo` or a negligible `H[(j, j - 1)]`, it chases the zero to the bottom with rotations on rows `jch, jch + 1` of `H` and `T`. `Q` must take the adjoint rotation on its columns, as the QZ sweep, the other zero-chase branch and the blocked path do (LAPACK `zhgeqz`: `ZROT(N, Q(1,JCH), 1, Q(1,JCH+1), 1, C, DCONJG(S))`); it rotated rows `jch, jch + 1` of `Q` with `s` instead. `Q` stayed unitary but `A = Q S Z^H` failed (`O(1)`), so every left eigenvector of `gevd_cplx` was wrong while the eigenvalues and the right eigenvectors, which use only `S`, `T` and `Z`, were right. It now rotates columns `jch, jch + 1` with `conj(s)`. | KEEP | tensor4all/faer#324 | Correctness only, no timing and no performance claim: the same number of rotations on a column pair instead of a row pair. Whether the branch runs depends on rounding: the issue's 12 x 12 `c32` pencil takes it on x86-64 with AVX2/FMA (left residual `1.962e-1` before, `2.079e-7` after; right `1.889e-7` both) and not on aarch64 or on x86-64 without AVX (Rosetta), where it passes before and after; `c64` does not take it. Tests: `linalg::gevd::qz_cplx::tests::unblocked_zero_chase_updates_q_columns` gives the unblocked iteration an exact zero at `j == ilo` and below a zero subdiagonal at `c32`/`c64`, `n = 2, 3, 6, 12`, and checks `Q`, `Z` unitary and `A = Q S Z^H`, `B = Q T Z^H` to `10 n eps` (fails before on aarch64: `2.51` against `7.0e-15`); `faer/tests/gevd_infinite_eigenvalue.rs` embeds the issue's bit patterns and checks the left and right residuals at `c32` and `c64` (fails before on x86-64 AVX2 only). The same code is in upstream faer `main` `5833f4e` at `faer/src/linalg/gevd/qz_cplx/mod.rs:123-127`, from `57e25ac` (generalized evd initial impl); verified by reading that source, not by running it, and no `upstream-reports/` entry is added (the deferred-feedback decision of `AGENTS.md`). | ## Ground rules diff --git a/faer/src/linalg/gevd/qz_cplx/mod.rs b/faer/src/linalg/gevd/qz_cplx/mod.rs index a4a56c41..b53cd803 100644 --- a/faer/src/linalg/gevd/qz_cplx/mod.rs +++ b/faer/src/linalg/gevd/qz_cplx/mod.rs @@ -120,11 +120,18 @@ fn hessenberg_to_qz_unblocked( .get_mut(.., jch + 1..ilastm + 1) .two_rows_mut(jch, jch + 1); rot(c.copy(), s.copy(), x, y); - if let Some(mut Q) = Q.rb_mut() { - let (x, y) = Q - .rb_mut() - .two_rows_mut(jch, jch + 1); - rot(c.copy(), s.copy(), x, y); + // the rotation acts on rows `jch, jch + 1` of `H` + // and `T`, so `Q` takes its adjoint on the right: + // columns `jch, jch + 1`, with `conj(s)` (faer#324) + if let Some(Q) = Q.rb_mut() { + let (x, y) = + Q.two_cols_mut(jch, jch + 1); + rot( + c.copy(), + s.conj(), + x.transpose_mut(), + y.transpose_mut(), + ); } if ilazr2 { H[(jch, jch - 1)] = @@ -1574,6 +1581,91 @@ mod tests { } } } + /// faer#324: a zero on the diagonal of `B` at `j`, with `j == ilo` or a + /// negligible `A[(j, j - 1)]`, is chased to the bottom by rotations applied + /// from the left. those rotations must update the columns of `Q`; they used + /// to rotate its rows, so `Q` stayed unitary but no longer satisfied + /// `A = Q S Z^H`, and every left eigenvector built from it was wrong. + fn check_unblocked_zero_chase( + A: Mat, + B: Mat, + zero_at: usize, + ) { + let n = A.nrows(); + let mut A = A; + let mut B = B; + for j in 0..n { + for i in j + 2..n { + A[(i, j)] = zero(); + } + for i in j + 1..n { + B[(i, j)] = zero(); + } + } + B[(zero_at, zero_at)] = zero(); + if zero_at > 0 { + A[(zero_at, zero_at - 1)] = zero(); + } + let mut Q = Mat::::identity(n, n); + let mut Z = Mat::::identity(n, n); + let mut S = A.clone(); + let mut P = B.clone(); + let mut alpha = Col::::zeros(n); + let mut beta = Col::::zeros(n); + hessenberg_to_qz_unblocked( + 0, + n - 1, + S.as_mut(), + P.as_mut(), + Some(Q.as_mut()), + Some(Z.as_mut()), + alpha.as_mut(), + beta.as_mut(), + false, + ); + // LAPACK's QZ tests divide by `n * ulp` and accept ratios up to 30 + let ref tol = from_f64::(10.0 * n as f64) * eps::(); + let I = Mat::::identity(n, n); + assert!((Q.adjoint() * &Q - &I).norm_max() < *tol); + assert!((Z.adjoint() * &Z - &I).norm_max() < *tol); + assert!((&Q * &S * Z.adjoint() - &A).norm_max() < tol * A.norm_max()); + assert!((&Q * &P * Z.adjoint() - &B).norm_max() < tol * B.norm_max()); + for j in 0..n { + for i in j + 1..n { + ::core::assert!( + S[(i, j)] == zero::() && P[(i, j)] == zero::() + ); + } + } + ::core::assert_eq!( + beta.iter().filter(|b| **b == zero::()).count(), + 1 + ); + } + #[test] + fn unblocked_zero_chase_updates_q_columns() { + let rng = &mut StdRng::seed_from_u64(324); + for n in [2, 3, 6, 12] { + let rand = stats::CwiseMatDistribution { + nrows: n, + ncols: n, + dist: ComplexDistribution::new(StandardNormal, StandardNormal), + }; + let A: Mat = rand.sample(rng); + let B: Mat = rand.sample(rng); + let to_c32 = |M: &Mat| { + Mat::from_fn(n, n, |i, j| { + c32::new(M[(i, j)].re as f32, M[(i, j)].im as f32) + }) + }; + // `zero_at == 0` takes the `j == ilo` path, `zero_at == n / 2` the + // negligible subdiagonal path + for zero_at in [0, n / 2] { + check_unblocked_zero_chase(A.clone(), B.clone(), zero_at); + check_unblocked_zero_chase(to_c32(&A), to_c32(&B), zero_at); + } + } + } #[test] fn test_qz_cplx_blocked() { let rng = &mut StdRng::seed_from_u64(0); diff --git a/faer/tests/gevd_infinite_eigenvalue.rs b/faer/tests/gevd_infinite_eigenvalue.rs new file mode 100644 index 00000000..9086e95a --- /dev/null +++ b/faer/tests/gevd_infinite_eigenvalue.rs @@ -0,0 +1,226 @@ +//! faer#324: `gevd_cplx` returned wrong left eigenvectors for a pencil with an +//! infinite eigenvalue, while its eigenvalues and right eigenvectors were right. +//! +//! When the complex QZ iteration (`hessenberg_to_qz_unblocked`, the path for +//! `n < 15`) finds a zero on the diagonal of `B` at `j == ilo` or below a +//! negligible subdiagonal of `A`, it chases that zero to the bottom with +//! rotations applied to rows `jch, jch + 1` of `A` and `B`. `Q` must then +//! take the adjoint rotation on columns `jch, jch + 1`; it rotated rows +//! instead. `Q` stayed unitary but `A = Q S Z^H` no longer held, so every +//! left eigenvector `Q * y` was wrong. `Z` was untouched, so the right +//! eigenvectors were fine. +//! +//! Whether the iteration takes that branch depends on rounding: the issue's +//! `c32` pencil takes it on x86-64 with AVX2/FMA and not on aarch64 or on +//! x86-64 without AVX, so this test guards the x86-64 CI path. The branch +//! itself is covered on every target by +//! `linalg::gevd::qz_cplx::tests::unblocked_zero_chase_updates_q_columns`, +//! which hands the unblocked iteration an exact zero on the diagonal of `B`. +//! +//! The test checks, for every column, the relative residuals +//! `||beta u^H A - alpha u^H B|| / (||u|| (|beta| ||A|| + |alpha| ||B||))` and +//! `||beta A v - alpha B v|| / (||v|| (|beta| ||A|| + |alpha| ||B||))`, +//! computed in `c64`, against `10 * n * eps` of the scalar type. +use faer::diag::Diag; +use faer::dyn_stack::{MemBuffer, MemStack}; +use faer::linalg::evd::ComputeEigenvectors; +use faer::linalg::gevd::{gevd_cplx, gevd_scratch}; +use faer::traits::ComplexField; +use faer::{Mat, Par, Scale, c32, c64}; + +/// `A` and `B` of faer#324, 12 x 12 column-major `c32` as `(re, im)` f32 bit +/// patterns. Made in Julia by `rng = Xoshiro(23)`, `U` and `V` the Q factors of +/// `randn(rng, ComplexF32, 12, 12)`, `A = randn(rng, ComplexF32, 12, 12)` and +/// `B = U * Diagonal([ones(Float32, 11); 0]) * V'`. +const A_BITS: [u32; 288] = [ + 0xbf4f6163, 0xbfaf8d0f, 0xbf284eb8, 0x3e5450f6, 0x3f792599, 0x3f468b81, + 0xbfd6e280, 0x3e5ac4dd, 0x3e801a4f, 0x3ed3389f, 0x3f3b7e67, 0x3ee2cce8, + 0xbf18d76b, 0xbfec8d86, 0x3dfade19, 0x3ea50f5a, 0xbe51783e, 0x3ebe98a9, + 0xbec60d56, 0xbfcd7c07, 0xbf391480, 0xbcdd4260, 0x3f2eb1f3, 0x3edba827, + 0x3eb9a75f, 0x3eabe882, 0x3ef40f6c, 0x3f170c99, 0xbf469536, 0xbd56f76f, + 0x3ea18fa1, 0xbf3f9766, 0x3f39edcf, 0x3f2b1e16, 0xbf8cc624, 0xbf00995e, + 0xbf848bb3, 0x3e452b55, 0x3f07d692, 0xbf15e47f, 0xbf78c96b, 0x3dac89b0, + 0x3f7b5738, 0xbdbf1b5d, 0xbee59f63, 0x3f126a87, 0xbf10c015, 0x40010505, + 0xbf3eb471, 0xbf0e8290, 0xbedcd6f8, 0x3f8fd58d, 0xbf9e60eb, 0xbec11184, + 0xbf84f3ec, 0xbf86c7a2, 0x3f32ee92, 0xbe982796, 0x3ebdfbfc, 0xbdd556d5, + 0x3f813d6e, 0xbe680d0d, 0xbf8b6f2f, 0x3ed8a27d, 0xbef1f445, 0xbf228e89, + 0xbe0dbcda, 0xbe5c7b1e, 0x3dcd9628, 0x3ed9b4eb, 0xbd8ba4cf, 0xbe5272d0, + 0xbe6f596a, 0xbe0c7a03, 0x3eede3c1, 0x3e97e81b, 0xbef31572, 0xbf089f76, + 0xbed23891, 0xbeb51c5a, 0xbeb03be6, 0x3f98a2cf, 0xbea89963, 0xbed6e037, + 0xbed59ec8, 0x3ea35097, 0xbf20ced1, 0xbfb19680, 0x3d6d6559, 0xbed0f1e4, + 0xbeade97f, 0x3ddf8a4e, 0xbfe078b6, 0x3d60fda8, 0xbf8e0d51, 0x3f363eb8, + 0xbe77a522, 0xbf6625ee, 0xbf4d0900, 0xbdde8b22, 0xbe6923d1, 0x3ddbb595, + 0x3c9e5539, 0xbecf00a1, 0xbc918357, 0xbf9d3f59, 0xbf5e30c5, 0x3f9699f5, + 0x3d27c04e, 0xbee2252a, 0x3f85e4e5, 0xbfc15574, 0x3e4ea842, 0x3f2f37c4, + 0xbf8bce39, 0xbf9317e6, 0xbf3531f5, 0x3f7276ad, 0xbe91378b, 0x3f32558b, + 0x3f2246fb, 0x3e2b5f34, 0xbee5ef32, 0x3c876061, 0xbedd50c2, 0x3f16c641, + 0xbf8ee74f, 0x3df84626, 0x3f4285c0, 0x3f8c59c3, 0xbe96eada, 0x3e116f7f, + 0x3f558f6c, 0xbbc2eb96, 0x3e9dbe69, 0xbf847a15, 0x3ee472ae, 0x3fad3734, + 0x3f51c643, 0x3e886971, 0xbedb1ddd, 0xbefa6b43, 0x3f091bbd, 0xbf73555a, + 0x3f406c28, 0xbe96084c, 0x3e841a5b, 0x4002c851, 0xbf60d371, 0xbf879afb, + 0x3f0316a2, 0x3d39df0b, 0x3eb84e69, 0x3e8d8780, 0x3e1fc026, 0xbe67b5f9, + 0xbee8ca57, 0xbf523235, 0xbe983fad, 0x3fa56928, 0xbfef51dd, 0x3ef021bc, + 0x3bf816ac, 0xbf083b99, 0x3eadcc26, 0xbdd7c5d8, 0xbe9fbcc4, 0x3f1576cd, + 0x3e867594, 0x3f43aea7, 0xbf6e2796, 0x3e1d360e, 0xbf6e38c8, 0x3f2046ab, + 0x3fbf11c5, 0xbfbfc438, 0x3f99b1ee, 0xbc612ab1, 0x3ec37e6f, 0x3e7754d3, + 0xbb6cb962, 0xbf73023d, 0x3d91c1f4, 0x3e817cf4, 0x3ce5ebcb, 0xbf2105dc, + 0x3db5639c, 0xbf8052b5, 0xbf440e73, 0xbec72e16, 0xbf3e555d, 0xbf677c26, + 0x3edf571d, 0x3f157aad, 0xbfa69ded, 0xbf30c8b2, 0xbf359983, 0xbdedffdb, + 0x3f5b4fed, 0xbe725d63, 0x3e4a4fd9, 0xbf3a1e8e, 0x3eced515, 0x3d2fe6e3, + 0x3dae1562, 0x3eccd06e, 0xbf4f71bc, 0x3dee43b8, 0x3eaec0da, 0xbe238bad, + 0x3e1921ae, 0x3f8d09a0, 0xbf3a55d5, 0xbf047854, 0x3f8e4e32, 0x3fae52b6, + 0x3f28457c, 0xbc905c85, 0x3f139872, 0xbf5e587a, 0xbf2fecf0, 0x3e07d7a8, + 0x3e197c0e, 0x3f12897a, 0x3fa47f64, 0x3ddf12b2, 0x3e9abf88, 0xbe4869d7, + 0x3f89204a, 0x3e9de574, 0x3da829bd, 0xba92e9c3, 0x3eb443ae, 0x3f893eef, + 0x3ea84a00, 0xbdd90d75, 0x3f3de106, 0xbe28e6ae, 0x3e74f4b7, 0x3f361e14, + 0xbf6f0f53, 0xbf77ef24, 0xbf85828d, 0xbf40f1a9, 0x3ea16a7e, 0xbfa51ce5, + 0x3beb5de9, 0x3f45e6dc, 0xbf3ba6a4, 0x3f5134b5, 0xbc3c375b, 0xbf3ef52a, + 0x3c8b5053, 0xbeb80981, 0x3f40ac6a, 0x3f1b1eb7, 0x3f54b29d, 0x3f23070c, + 0xc02b5709, 0xbe80d831, 0xbf830608, 0x3d6bbfbf, 0x3e877dd7, 0xbe190d19, + 0x3f54c252, 0x3f3c7e74, 0xbeb1582b, 0xbf0cda46, 0x3f262f3b, 0x3ea8b9a9, + 0x3eacc3db, 0xbeb5fb1b, 0x3e1f302a, 0xbde174ef, 0x3eeea84d, 0xbf6eec14, + 0xbf3444f5, 0xbf1df5eb, 0x3f22c2c6, 0x3f8c34c9, 0x3efdf863, 0x3d1f5bcd, + 0x3f2a0e32, 0xbf2f3100, 0x3ea0e062, 0x3d961833, 0xbf885b14, 0xbe826741, +]; +const B_BITS: [u32; 288] = [ + 0xbe214c0a, 0x3c180ab2, 0x3e53942f, 0xbe2d70a4, 0x3d7eace2, 0x3ee59f96, + 0xbdaa50ce, 0xbe2bd95e, 0xbe3f58ca, 0xbe05a0c8, 0x3ef358df, 0xbe83a756, + 0xbe8fdf06, 0x3d48267f, 0x3e149e40, 0x3e096784, 0xbdaadec9, 0xbdfa972b, + 0x3e984722, 0xbb3e4340, 0x3def0a2c, 0xbc41c210, 0x3e276ff4, 0x3be8f226, + 0x3de57fbc, 0x3ebbd2e2, 0x3e4dbeb1, 0x3da2dab4, 0xbdf5d5b0, 0x3e05ec2d, + 0x3e049296, 0x3e104455, 0xbe8b4ac4, 0xbe4a33f6, 0xbd837fcc, 0xbdecdcc3, + 0x3efaf2b9, 0xbda655c6, 0x3e6281e7, 0xbc4e524c, 0x3e766e16, 0xbe4565eb, + 0xbdd643d4, 0xbe014954, 0xbe126970, 0xbdc8f2b7, 0x3eaf43a6, 0x3e3f3c1d, + 0x3de8c144, 0xbc3926b6, 0xbeafd6ac, 0x3d744101, 0xbdd2b237, 0x3e2221fe, + 0xbd8ff132, 0xbe86f6c0, 0x3ea7e3b2, 0x3ed46ff9, 0xbd2550bc, 0xbea942b3, + 0x3e33eca9, 0xbe269600, 0xbe036129, 0x3e070a86, 0x3e87b637, 0xbd8cc3a4, + 0x3e34b6e9, 0xbcabfb6c, 0xbd2d2511, 0xbdbd2efd, 0x3cc7281c, 0x3e4b2214, + 0x3df2d99f, 0x3e40d4df, 0x3ec5fdac, 0x3e93ac9f, 0x3ec75b99, 0xbe1eaa3e, + 0x3da2549b, 0xbe490acb, 0x3d66abd0, 0x3c3f094b, 0x3de20da6, 0xbd9efde0, + 0x3c12aed3, 0xbe52fbb5, 0xbd9b8bbb, 0xbde4e7ae, 0x3eaa4ae8, 0x3ebd3359, + 0x3e67831e, 0x3e4f015c, 0xbd98e34b, 0xba928700, 0xbdbe9ac8, 0xbd1166ed, + 0xbdc5900c, 0x3e5f9b7e, 0xbe362eee, 0xbc9b602e, 0xbdf7c275, 0xbe28665d, + 0x3ea67435, 0x3e304614, 0x3e90fa5e, 0xbd103f64, 0x3eb1bed0, 0xbe474954, + 0xbe437fa1, 0xbe22a1e1, 0xbd953260, 0xbe05e26c, 0x3d4d570f, 0xbe8d674e, + 0xbe358232, 0x3d27dc00, 0xbcc7c07f, 0xbed99233, 0x3c4f2690, 0xbe44443d, + 0xbd665cb8, 0x3ec32396, 0x3c08c567, 0xbe142bc2, 0xbeb5bfe9, 0x3d687635, + 0x3d87f331, 0xbdfe645e, 0xbdb0537c, 0x3e3a031d, 0xbe80bde0, 0x3d7e0321, + 0xbd60dde7, 0xbdcf81ae, 0x3e6fc8d6, 0x3e4e0e0e, 0xbd119fd0, 0x3b566150, + 0xbcd944ca, 0x3db25f92, 0x3e14174c, 0x3e446615, 0xbe9c0a74, 0x3db603cd, + 0x3d960d0c, 0x3e952cdf, 0xbea1bbd4, 0xbe36be2f, 0x3e46a5f1, 0xbd3fc041, + 0x3e8a9a2c, 0x3e8750e5, 0x3dc74849, 0xbe6cba24, 0xbd049520, 0x3e1df621, + 0xbe239c0c, 0xbd44ce9b, 0x3e16b6d3, 0xbe39eb22, 0x3d7421f0, 0xbd81a260, + 0x3ee11557, 0xbd2ba59d, 0xbe06ec40, 0x3e95ccae, 0xbe1455c7, 0x3dad765a, + 0x3dfc3786, 0xbe0a1c67, 0x3e91699f, 0xbe7f44ab, 0x3e35997a, 0x3e56074a, + 0x3e48ba1f, 0xbda86bb0, 0x3e619686, 0x3d6b51ee, 0x3de74721, 0x3e70e215, + 0x3e1e871e, 0x3e8480af, 0x3c71663b, 0xbcfdecf8, 0x3c790499, 0xbdda71ec, + 0xbe69086f, 0xbe666ca7, 0xbeae5070, 0xbe32facd, 0xbebf909a, 0x3e0acc7a, + 0xbe09fa07, 0xbd188d12, 0xbe586622, 0x3e07db2b, 0xbe0d881f, 0x3d8510b4, + 0xbe12139c, 0xbe607175, 0xbd4383ca, 0x3c6848f5, 0x3e7fedf8, 0x3d73fd22, + 0xbd4cd5a5, 0xbe0075a8, 0x3ebe613e, 0xbd4b5a62, 0x3e928ef4, 0x3e13ef02, + 0xbe0fbbdb, 0xbec3be67, 0xbea68554, 0x3e69adc4, 0xbe0308c2, 0xbeaf1abc, + 0x3e4ff81a, 0xbdf6996b, 0xbcf3e916, 0x3da93af8, 0x3dd91840, 0x3dffefc5, + 0x3d8f6344, 0x3e506fe1, 0x3db32152, 0xbd7a5993, 0x3e443a4b, 0xbdbcd106, + 0x3ef0897e, 0xbd28b80c, 0xbd8e2c30, 0x3eab0c19, 0xbda095fc, 0xbe54b551, + 0x3e59276c, 0xbd096141, 0x3e7a8623, 0x3e5fab98, 0xbe9d6296, 0xbea3a79e, + 0x3ec1b7cb, 0xbeaf0bd4, 0x3d6ba9d5, 0x3daa1207, 0xbe3755bc, 0xbe18176b, + 0x3de716ca, 0x3ecb1289, 0x3d346a84, 0x3e170f90, 0x3e1075aa, 0xbe8bcacc, + 0xbe0f039f, 0x3e19172a, 0x3ed7aead, 0x3db7bf46, 0x3cba2610, 0x3e448594, + 0xbd80054b, 0x3e3de570, 0xbe09613f, 0x3e0b7ea0, 0x3dc99b48, 0x3dd54d7d, + 0xbe85bba1, 0xbdb9fe3e, 0x3c123b6f, 0x3e4e185d, 0x3df3054a, 0x3e39ebc0, + 0x3e42fed0, 0x3ea190a8, 0xbe6dd6ee, 0x3ebb508b, 0xbc42d4bd, 0x3caa4f60, + 0xbe40490b, 0xbdc9682a, 0xbe9a10be, 0x3e13860e, 0x3ccfbf18, 0xbdce51b2, + 0x3dbc018d, 0xbe434aac, 0xbe9bc951, 0x3e1ff480, 0x3d1472b4, 0x3e9a2e82, +]; + +fn from_bits(bits: &[u32; 288]) -> Mat { + Mat::from_fn(12, 12, |i, j| { + let k = 2 * (i + 12 * j); + c32::new(f32::from_bits(bits[k]), f32::from_bits(bits[k + 1])) + }) +} + +fn to_c32(m: &Mat) -> Mat { + Mat::from_fn(m.nrows(), m.ncols(), |i, j| { + c32::new(m[(i, j)].re as f32, m[(i, j)].im as f32) + }) +} + +fn to_c64(m: &Mat) -> Mat { + Mat::from_fn(m.nrows(), m.ncols(), |i, j| { + c64::new(m[(i, j)].re as f64, m[(i, j)].im as f64) + }) +} + +/// returns the largest left and right relative residuals over all columns, and +/// the smallest `|beta| / (|alpha| + |beta|)` +fn residuals( + a: &Mat, + b: &Mat, + up: impl Fn(&T) -> c64, +) -> (f64, f64, f64) { + let n = a.nrows(); + let (mut a1, mut b1) = (a.clone(), b.clone()); + let mut s = Diag::::zeros(n); + let mut beta = Diag::::zeros(n); + let mut ul = Mat::::zeros(n, n); + let mut ur = Mat::::zeros(n, n); + let params = Default::default(); + let mut mem = MemBuffer::new(gevd_scratch::( + n, + ComputeEigenvectors::Yes, + ComputeEigenvectors::Yes, + Par::Seq.shape(), + params, + )); + gevd_cplx( + a1.as_mut(), + b1.as_mut(), + s.as_mut(), + beta.as_mut(), + Some(ul.as_mut()), + Some(ur.as_mut()), + Par::Seq, + MemStack::new(&mut mem), + params, + ) + .unwrap(); + let a = Mat::from_fn(n, n, |i, j| up(&a[(i, j)])); + let b = Mat::from_fn(n, n, |i, j| up(&b[(i, j)])); + let (na, nb) = (a.norm_l2(), b.norm_l2()); + let (mut left, mut right, mut beta_min) = (0.0f64, 0.0f64, f64::INFINITY); + for j in 0..n { + let (al, be) = (up(&s[j]), up(&beta[j])); + beta_min = beta_min.min(be.norm() / (al.norm() + be.norm())); + let u = Mat::from_fn(n, 1, |i, _| up(&ul[(i, j)])); + let v = Mat::from_fn(n, 1, |i, _| up(&ur[(i, j)])); + let lu = u.adjoint() * &a * Scale(be) - u.adjoint() * &b * Scale(al); + let rv = &a * &v * Scale(be) - &b * &v * Scale(al); + let den = be.norm() * na + al.norm() * nb; + left = left.max(lu.norm_l2() / (u.norm_l2() * den)); + right = right.max(rv.norm_l2() / (v.norm_l2() * den)); + } + (left, right, beta_min) +} + +fn check(a: Mat, b: Mat) { + let n = a.nrows() as f64; + let up32 = |x: &c32| c64::new(x.re as f64, x.im as f64); + let (l, r, beta_min) = residuals(&to_c32(&a), &to_c32(&b), up32); + let tol = 10.0 * n * f32::EPSILON as f64; + assert!(l < tol, "c32 left residual {l:e} >= {tol:e}"); + assert!(r < tol, "c32 right residual {r:e} >= {tol:e}"); + // `B` is singular up to the rounding of its entries to f32 + assert!(beta_min < tol, "c32 smallest beta {beta_min:e} >= {tol:e}"); + let (l, r, beta_min) = residuals(&a, &b, |x: &c64| *x); + let tol = 10.0 * n * f64::EPSILON; + assert!(l < tol, "c64 left residual {l:e} >= {tol:e}"); + assert!(r < tol, "c64 right residual {r:e} >= {tol:e}"); + assert!(beta_min < 10.0 * n * f32::EPSILON as f64); +} + +#[test] +fn left_eigenvectors_of_issue_324_pencil() { + let (a, b) = (from_bits(&A_BITS), from_bits(&B_BITS)); + check(to_c64(&a), to_c64(&b)); +} From f4948f499b8cc9331cb4a9868cd7d3b56e90e89a Mon Sep 17 00:00:00 2001 From: Hiroshi Shinaoka Date: Mon, 9 Nov 2026 14:21:43 +0900 Subject: [PATCH 2/5] test(gevd): make the #324 regression test fail on a NaN residual, and fix the size claim `residuals()` in `faer/tests/gevd_infinite_eigenvalue.rs` folded each column into its maxima with `f64::max`, which returns the other operand when one is `NaN`. An all-zero left eigenvector column gives `0 / 0`, so the `NaN` was dropped and the column read as perfectly accurate: with `ul` forced to zeros the test passed. Each column now asserts, before it is folded in, that both residual numerators and `|beta| / (|alpha| + |beta|)` are finite and that `||u||`, `||v||` and the residual denominator are non-zero, reporting the column index and values. The same injection now fails at `column 0: left eigenvector has norm 0e0`. The `n < 15` size claim was wrong. `hessenberg_to_qz_blocked` calls `hessenberg_to_qz_unblocked` and returns when `n < nmin`, and calls it again at the end on the whole range when it does not, so the branch is reachable at every `n`. `nmin` is `Ord::max(15, blocking_threshold)` and the default `blocking_threshold` is 75, so with the default parameters the whole iteration is unblocked for `n < 75`. `CHANGELOG.md`, `PATCHSET.md` and the test's module comment now say that. No change to the rotation logic. --- CHANGELOG.md | 3 +- PATCHSET.md | 2 +- faer/tests/gevd_infinite_eigenvalue.rs | 55 ++++++++++++++++++++------ 3 files changed, 46 insertions(+), 14 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 9a1f539a..393b51d7 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,7 +13,8 @@ Changes since `t4a-faer` 0.32.0 (tag `faer-v0.32.0`). Numbers are tensor4all/fae ## Correctness - `gevd_cplx` returned wrong left eigenvectors for some pencils with an infinite eigenvalue (a zero on the diagonal of `B` during the QZ - iteration, `n < 15`), while the eigenvalues and the right eigenvectors were right. The rotations that chase that zero to the bottom were + iteration's unblocked path, which the default parameters run for `n < max(15, blocking_threshold)`, i.e. `n < 75`, and which the blocked path + also ends in), while the eigenvalues and the right eigenvectors were right. The rotations that chase that zero to the bottom were applied to the rows of `Q` instead of its columns. Whether a pencil reaches that branch depends on rounding, so the same input could fail on x86-64 with AVX2 and pass on aarch64 (#324). diff --git a/PATCHSET.md b/PATCHSET.md index 5a1529e3..da6fca6a 100644 --- a/PATCHSET.md +++ b/PATCHSET.md @@ -193,7 +193,7 @@ gate has changed and two kinds of adoption are now in the ledger. | `norm_max` propagates NaN (#316) | `linalg::reductions::norm_max` (and `MatRef`/`ColRef`/`RowRef::norm_max`) skipped NaN entries: the SIMD kernel accumulates with `simd.max`, a `select(acc > x, acc, x)` that keeps the non-NaN operand, and the pairwise halves and the scalar loop combine with `fmax` (IEEE `maxNum`). It returned the largest finite magnitude, and 0 for an all-NaN input, while `norm_l1`/`norm_l2` returned NaN. It now returns NaN when any entry, or either part of a complex entry, is NaN, as LAPACK's `xLANGE('M')` and Julia's `norm(x, Inf)` do. The SIMD kernel keeps its max accumulators and adds a second set that sums the magnitudes: every term is `>= 0`, `+inf` or NaN, so the sum is NaN exactly when an entry is NaN, and overflow gives `inf`, not NaN. The pairwise combine propagates NaN, and the scalar loop returns NaN at the first NaN entry. An empty matrix still returns 0. This changes behaviour: internal callers that take `norm_max` of a matrix holding NaN (QZ, tridiagonal EVD, bidiagonal SVD, Bunch-Kaufman scaling) now see NaN instead of a finite norm. Both files are byte-identical to the upstream sync base `0ecd846`, so upstream faer carries the same behaviour there; whether later upstream changed it is not checked, and feedback is deferred. | EXPERIMENT | tensor4all/faer#316 | Correctness only, no timing claim. The SIMD inner loop does one extra vector add per batch. New tests in `linalg::reductions::norm_max::tests`: f32/f64/c32/c64 with NaN in a real entry and in the real or imaginary part of a complex entry, at the first, middle and last position, and all-NaN, across shapes on both sides of the SIMD tail and the pairwise split. Each case goes through the contiguous (SIMD) path and a row-stride-2 view (scalar loop), and through the `MatRef`/`ColRef`/`RowRef` methods. Finite and `inf` inputs give the same result on both paths, and `norm_l1`/`norm_l2` agree on NaN. Removing the SIMD sum line or the scalar NaN return fails `test_norm_max_propagates_nan`. | | `stats` column/row reductions with a reversed output axis (#317) | `col_mean`/`col_varm`, and through them `row_mean`/`row_varm`, normalize a negative row stride with `mat.reverse_rows()` before the kernel runs. Rows are the output axis of a column reduction, and `out` (and, for `varm`, `col_mean`) were not reversed with them. Entry `i` was written to `nrows - 1 - i`, and `varm` subtracted the mean of the mirrored row. All four internal functions (`col_{mean,varm}_{ignore,propagate}`) now reverse `out`, and `col_mean` for `varm`, together with `mat`. Reversing the reduced axis (`reverse_cols` in `col_*`) was already harmless. Public API unchanged. `meanvar.rs` is byte-identical to the upstream sync base `0ecd846`, so upstream faer carries the same defect there; later upstream is not checked, and feedback is deferred. | EXPERIMENT | tensor4all/faer#317 | Correctness only, no timing claim. New tests in `stats::meanvar::tests`: `col_*` on `reverse_rows()` and `row_*` on `reverse_cols()`, mean and varm, `Propagate` and `Ignore`, f32/f64/c32/c64, on column-major (fallback kernel) and row-major (row-major SIMD kernel) storage with one NaN. Each is compared with the same matrix copied to an owned one. The issue's reproducer is pinned as well. Removing the output reversal from any one of the four functions fails both tests. | | Complex `abs_max` drops a NaN in the real part (#320) | `ComplexField::simd_abs_max` takes the larger of the two part magnitudes with an unordered `re > im` select, and that select returns `im` when `re` is NaN: a NaN real part read as the finite imaginary magnitude, while a NaN imaginary part did return NaN. The imaginary side already lands on `im`, so the real side takes an explicit `!(re == re)` mask in the generic `Complex` SIMD path, and the scalar fallbacks of `ComplexImpl` / `ComplexImpl` test both parts before the comparison. No in-tree type reaches either path (native `c32` / `c64` go through the real coercion in `norm_max`), so no dispatch, allocation or numerical result on an existing input changes. | KEEP | tensor4all/faer#320 | Correctness; the `faer-traits` unit test `complex_abs_max_propagates_nan` drives `SimdCtx::abs_max` on `Complex` and fails before / passes after (real-part NaN, imaginary-part NaN, both, and three finite controls). No performance claim: two mask ops per SIMD chunk on a path no in-tree type reaches. | -| Complex QZ: a zero chased down `B`'s diagonal updates the columns of `Q` (#324) | When `hessenberg_to_qz_unblocked` (the complex QZ path for `n < 15`) finds `T[(j, j)]` negligible with `j == ilo` or a negligible `H[(j, j - 1)]`, it chases the zero to the bottom with rotations on rows `jch, jch + 1` of `H` and `T`. `Q` must take the adjoint rotation on its columns, as the QZ sweep, the other zero-chase branch and the blocked path do (LAPACK `zhgeqz`: `ZROT(N, Q(1,JCH), 1, Q(1,JCH+1), 1, C, DCONJG(S))`); it rotated rows `jch, jch + 1` of `Q` with `s` instead. `Q` stayed unitary but `A = Q S Z^H` failed (`O(1)`), so every left eigenvector of `gevd_cplx` was wrong while the eigenvalues and the right eigenvectors, which use only `S`, `T` and `Z`, were right. It now rotates columns `jch, jch + 1` with `conj(s)`. | KEEP | tensor4all/faer#324 | Correctness only, no timing and no performance claim: the same number of rotations on a column pair instead of a row pair. Whether the branch runs depends on rounding: the issue's 12 x 12 `c32` pencil takes it on x86-64 with AVX2/FMA (left residual `1.962e-1` before, `2.079e-7` after; right `1.889e-7` both) and not on aarch64 or on x86-64 without AVX (Rosetta), where it passes before and after; `c64` does not take it. Tests: `linalg::gevd::qz_cplx::tests::unblocked_zero_chase_updates_q_columns` gives the unblocked iteration an exact zero at `j == ilo` and below a zero subdiagonal at `c32`/`c64`, `n = 2, 3, 6, 12`, and checks `Q`, `Z` unitary and `A = Q S Z^H`, `B = Q T Z^H` to `10 n eps` (fails before on aarch64: `2.51` against `7.0e-15`); `faer/tests/gevd_infinite_eigenvalue.rs` embeds the issue's bit patterns and checks the left and right residuals at `c32` and `c64` (fails before on x86-64 AVX2 only). The same code is in upstream faer `main` `5833f4e` at `faer/src/linalg/gevd/qz_cplx/mod.rs:123-127`, from `57e25ac` (generalized evd initial impl); verified by reading that source, not by running it, and no `upstream-reports/` entry is added (the deferred-feedback decision of `AGENTS.md`). | +| Complex QZ: a zero chased down `B`'s diagonal updates the columns of `Q` (#324) | When `hessenberg_to_qz_unblocked` (the unblocked complex QZ path, which the default parameters run for `n < max(15, blocking_threshold)`, i.e. `n < 75`, and which the blocked path also calls at the end) finds `T[(j, j)]` negligible with `j == ilo` or a negligible `H[(j, j - 1)]`, it chases the zero to the bottom with rotations on rows `jch, jch + 1` of `H` and `T`. `Q` must take the adjoint rotation on its columns, as the QZ sweep, the other zero-chase branch and the blocked path do (LAPACK `zhgeqz`: `ZROT(N, Q(1,JCH), 1, Q(1,JCH+1), 1, C, DCONJG(S))`); it rotated rows `jch, jch + 1` of `Q` with `s` instead. `Q` stayed unitary but `A = Q S Z^H` failed (`O(1)`), so every left eigenvector of `gevd_cplx` was wrong while the eigenvalues and the right eigenvectors, which use only `S`, `T` and `Z`, were right. It now rotates columns `jch, jch + 1` with `conj(s)`. | KEEP | tensor4all/faer#324 | Correctness only, no timing and no performance claim: the same number of rotations on a column pair instead of a row pair. Whether the branch runs depends on rounding: the issue's 12 x 12 `c32` pencil takes it on x86-64 with AVX2/FMA (left residual `1.962e-1` before, `2.079e-7` after; right `1.889e-7` both) and not on aarch64 or on x86-64 without AVX (Rosetta), where it passes before and after; `c64` does not take it. Tests: `linalg::gevd::qz_cplx::tests::unblocked_zero_chase_updates_q_columns` gives the unblocked iteration an exact zero at `j == ilo` and below a zero subdiagonal at `c32`/`c64`, `n = 2, 3, 6, 12`, and checks `Q`, `Z` unitary and `A = Q S Z^H`, `B = Q T Z^H` to `10 n eps` (fails before on aarch64: `2.51` against `7.0e-15`); `faer/tests/gevd_infinite_eigenvalue.rs` embeds the issue's bit patterns and checks the left and right residuals at `c32` and `c64` (fails before on x86-64 AVX2 only); before folding a column into the maxima it asserts both residual numerators and `|beta|/(|alpha|+|beta|)` finite and `||u||`, `||v||` and the residual denominator non-zero, with the column number and values in the message, because `f64::max`/`f64::min` drop a `NaN` and the `0/0` of an all-zero left eigenvector column would otherwise read as a perfect residual (an all-zero `ul` passed the test before these asserts and fails after: `column 0: left eigenvector has norm 0e0, alpha = Complex { re: -2.6592135429382324, im: 2.4517483711242676 }, beta = Complex { re: 0.0, im: 0.0 }`). The same code is in upstream faer `main` `5833f4e` at `faer/src/linalg/gevd/qz_cplx/mod.rs:123-127`, from `57e25ac` (generalized evd initial impl); verified by reading that source, not by running it, and no `upstream-reports/` entry is added (the deferred-feedback decision of `AGENTS.md`). | ## Ground rules diff --git a/faer/tests/gevd_infinite_eigenvalue.rs b/faer/tests/gevd_infinite_eigenvalue.rs index 9086e95a..1b4aa970 100644 --- a/faer/tests/gevd_infinite_eigenvalue.rs +++ b/faer/tests/gevd_infinite_eigenvalue.rs @@ -1,14 +1,15 @@ //! faer#324: `gevd_cplx` returned wrong left eigenvectors for a pencil with an //! infinite eigenvalue, while its eigenvalues and right eigenvectors were right. //! -//! When the complex QZ iteration (`hessenberg_to_qz_unblocked`, the path for -//! `n < 15`) finds a zero on the diagonal of `B` at `j == ilo` or below a -//! negligible subdiagonal of `A`, it chases that zero to the bottom with -//! rotations applied to rows `jch, jch + 1` of `A` and `B`. `Q` must then -//! take the adjoint rotation on columns `jch, jch + 1`; it rotated rows -//! instead. `Q` stayed unitary but `A = Q S Z^H` no longer held, so every -//! left eigenvector `Q * y` was wrong. `Z` was untouched, so the right -//! eigenvectors were fine. +//! When the complex QZ iteration's unblocked path (`hessenberg_to_qz_unblocked`, +//! which the default parameters run for `n < max(15, blocking_threshold)`, i.e. +//! `n < 75`, and which the blocked path also calls at the end) finds a zero on +//! the diagonal of `B` at `j == ilo` or below a negligible subdiagonal of `A`, +//! it chases that zero to the bottom with rotations applied to rows +//! `jch, jch + 1` of `A` and `B`. `Q` must then take the adjoint rotation on +//! columns `jch, jch + 1`; it rotated rows instead. `Q` stayed unitary but +//! `A = Q S Z^H` no longer held, so every left eigenvector `Q * y` was wrong. +//! `Z` was untouched, so the right eigenvectors were fine. //! //! Whether the iteration takes that branch depends on rounding: the issue's //! `c32` pencil takes it on x86-64 with AVX2/FMA and not on aarch64 or on @@ -20,7 +21,8 @@ //! The test checks, for every column, the relative residuals //! `||beta u^H A - alpha u^H B|| / (||u|| (|beta| ||A|| + |alpha| ||B||))` and //! `||beta A v - alpha B v|| / (||v|| (|beta| ||A|| + |alpha| ||B||))`, -//! computed in `c64`, against `10 * n * eps` of the scalar type. +//! computed in `c64`, against `10 * n * eps` of the scalar type. A `NaN` +//! residual or `beta` fails the test instead of being dropped by `f64::max`. use faer::diag::Diag; use faer::dyn_stack::{MemBuffer, MemStack}; use faer::linalg::evd::ComputeEigenvectors; @@ -191,14 +193,43 @@ fn residuals( let (mut left, mut right, mut beta_min) = (0.0f64, 0.0f64, f64::INFINITY); for j in 0..n { let (al, be) = (up(&s[j]), up(&beta[j])); - beta_min = beta_min.min(be.norm() / (al.norm() + be.norm())); + // `f64::max` and `f64::min` drop a NaN operand, so a `0/0` residual + // would leave the maximum untouched and read as a perfect one. + let ratio = be.norm() / (al.norm() + be.norm()); + assert!( + ratio.is_finite(), + "column {j}: |beta|/(|alpha|+|beta|) = {ratio:e}, alpha = {al:?}, beta = {be:?}" + ); + beta_min = beta_min.min(ratio); let u = Mat::from_fn(n, 1, |i, _| up(&ul[(i, j)])); let v = Mat::from_fn(n, 1, |i, _| up(&ur[(i, j)])); let lu = u.adjoint() * &a * Scale(be) - u.adjoint() * &b * Scale(al); let rv = &a * &v * Scale(be) - &b * &v * Scale(al); let den = be.norm() * na + al.norm() * nb; - left = left.max(lu.norm_l2() / (u.norm_l2() * den)); - right = right.max(rv.norm_l2() / (v.norm_l2() * den)); + let (nlu, nrv) = (lu.norm_l2(), rv.norm_l2()); + let (nu, nv) = (u.norm_l2(), v.norm_l2()); + assert!( + nlu.is_finite(), + "column {j}: left residual numerator {nlu:e}, alpha = {al:?}, beta = {be:?}" + ); + assert!( + nrv.is_finite(), + "column {j}: right residual numerator {nrv:e}, alpha = {al:?}, beta = {be:?}" + ); + assert!( + nu > 0.0, + "column {j}: left eigenvector has norm {nu:e}, alpha = {al:?}, beta = {be:?}" + ); + assert!( + nv > 0.0, + "column {j}: right eigenvector has norm {nv:e}, alpha = {al:?}, beta = {be:?}" + ); + assert!( + den > 0.0, + "column {j}: residual denominator {den:e}, ||A|| = {na:e}, ||B|| = {nb:e}, alpha = {al:?}, beta = {be:?}" + ); + left = left.max(nlu / (nu * den)); + right = right.max(nrv / (nv * den)); } (left, right, beta_min) } From cdee0a4e3d61bcf09ad2232c85d349bae77a282f Mon Sep 17 00:00:00 2001 From: Hiroshi Shinaoka Date: Mon, 9 Nov 2026 14:34:05 +0900 Subject: [PATCH 3/5] test(gevd): assert the folded residuals finite, not only the numerators `left` and `right` are accumulated with `f64::max`, which returns the other operand when one is `NaN`; a `NaN` therefore leaves the maximum where it was and reads as a perfect residual. Each column now computes its two residuals, asserts both finite with the numerator, norm and denominator in the message, and only then folds them in. `ul = Mat::::zeros(n, n)` after the `gevd_cplx` call turns the left assert red and the mirror mutation on `ur` turns the right one red; both were green before. --- PATCHSET.md | 2 +- faer/tests/gevd_infinite_eigenvalue.rs | 16 ++++++++++++++-- 2 files changed, 15 insertions(+), 3 deletions(-) diff --git a/PATCHSET.md b/PATCHSET.md index 3969aef9..19298919 100644 --- a/PATCHSET.md +++ b/PATCHSET.md @@ -194,7 +194,7 @@ gate has changed and two kinds of adoption are now in the ledger. | `stats` column/row reductions with a reversed output axis (#317) | `col_mean`/`col_varm`, and through them `row_mean`/`row_varm`, normalize a negative row stride with `mat.reverse_rows()` before the kernel runs. Rows are the output axis of a column reduction, and `out` (and, for `varm`, `col_mean`) were not reversed with them. Entry `i` was written to `nrows - 1 - i`, and `varm` subtracted the mean of the mirrored row. All four internal functions (`col_{mean,varm}_{ignore,propagate}`) now reverse `out`, and `col_mean` for `varm`, together with `mat`. Reversing the reduced axis (`reverse_cols` in `col_*`) was already harmless. Public API unchanged. `meanvar.rs` is byte-identical to the upstream sync base `0ecd846`, so upstream faer carries the same defect there; later upstream is not checked, and feedback is deferred. | EXPERIMENT | tensor4all/faer#317 | Correctness only, no timing claim. New tests in `stats::meanvar::tests`: `col_*` on `reverse_rows()` and `row_*` on `reverse_cols()`, mean and varm, `Propagate` and `Ignore`, f32/f64/c32/c64, on column-major (fallback kernel) and row-major (row-major SIMD kernel) storage with one NaN. Each is compared with the same matrix copied to an owned one. The issue's reproducer is pinned as well. Removing the output reversal from any one of the four functions fails both tests. | | Explicit sparse column ordering stops paying for COLAMD (#322) | `factorize_symbolic_lu_with_column_ordering` and `factorize_symbolic_qr_with_column_ordering` build their internal `MemBuffer` from `StackReq::or(colamd::order_scratch::(m, n, A_nnz), ..)`, but COLAMD runs only in the `None` branch, so a caller that brings the ordering still allocates a COLAMD-sized workspace on every symbolic analysis (dyn-stack's `MemBuffer` takes it uninitialized, so the cost is the allocation, not a memset). Both functions include `order_scratch` only when `column_ordering.is_none()` now; the `None` request and the analysis request are unchanged. Both entry points are this fork's own (#239, `3889bd1`), so upstream faer does not carry the defect and no `upstream-reports/` entry is added. | KEEP | tensor4all/faer#322 | Allocation only, no numerical change: the two functions' permutations, symbolic structures and every downstream numeric path are untouched, and no benchmark was run. The issue's reproducer (n = 2000, nnz = 7990, identity ordering, `--release`, counting global allocator) goes from 332712 B largest allocation to 239928 B for simplicial LU and QR: the 332712 B was exactly `colamd::order_scratch`, and 239928 B is the analysis request (4 + 5 n-element arrays, the transpose arrays, and the largest simplicial arm), while the supernodal arm stays at 468080 B, which is the supernodal symbolic workspace rather than the COLAMD one. `faer/tests/sparse_symbolic_colamd_workspace.rs` pins both directions on a 100 x 100 sample: with an explicit ordering the largest allocation stays below `order_scratch`, without one it reaches it; the first fails before the change and both pass after. `cargo test -p t4a-faer -j6` passes except `lq_no_alloc::lq_allocates_nothing_60x300`, which fails identically on `origin/main` (dense LQ, unrelated), plus the process-wide allocation counters in `no_step_alloc` and `triangular_solve::aliased_rhs::staged_solve_is_bitwise_the_in_place_one`, which pass standalone and flake only inside the whole-suite parallel run. | | Complex `abs_max` drops a NaN in the real part (#320) | `ComplexField::simd_abs_max` takes the larger of the two part magnitudes with an unordered `re > im` select, and that select returns `im` when `re` is NaN: a NaN real part read as the finite imaginary magnitude, while a NaN imaginary part did return NaN. The imaginary side already lands on `im`, so the real side takes an explicit `!(re == re)` mask in the generic `Complex` SIMD path, and the scalar fallbacks of `ComplexImpl` / `ComplexImpl` test both parts before the comparison. No in-tree type reaches either path (native `c32` / `c64` go through the real coercion in `norm_max`), so no dispatch, allocation or numerical result on an existing input changes. | KEEP | tensor4all/faer#320 | Correctness; the `faer-traits` unit test `complex_abs_max_propagates_nan` drives `SimdCtx::abs_max` on `Complex` and fails before / passes after (real-part NaN, imaginary-part NaN, both, and three finite controls). No performance claim: two mask ops per SIMD chunk on a path no in-tree type reaches. | -| Complex QZ: a zero chased down `B`'s diagonal updates the columns of `Q` (#324) | When `hessenberg_to_qz_unblocked` (the unblocked complex QZ path, which the default parameters run for `n < max(15, blocking_threshold)`, i.e. `n < 75`, and which the blocked path also calls at the end) finds `T[(j, j)]` negligible with `j == ilo` or a negligible `H[(j, j - 1)]`, it chases the zero to the bottom with rotations on rows `jch, jch + 1` of `H` and `T`. `Q` must take the adjoint rotation on its columns, as the QZ sweep, the other zero-chase branch and the blocked path do (LAPACK `zhgeqz`: `ZROT(N, Q(1,JCH), 1, Q(1,JCH+1), 1, C, DCONJG(S))`); it rotated rows `jch, jch + 1` of `Q` with `s` instead. `Q` stayed unitary but `A = Q S Z^H` failed (`O(1)`), so every left eigenvector of `gevd_cplx` was wrong while the eigenvalues and the right eigenvectors, which use only `S`, `T` and `Z`, were right. It now rotates columns `jch, jch + 1` with `conj(s)`. | KEEP | tensor4all/faer#324 | Correctness only, no timing and no performance claim: the same number of rotations on a column pair instead of a row pair. Whether the branch runs depends on rounding: the issue's 12 x 12 `c32` pencil takes it on x86-64 with AVX2/FMA (left residual `1.962e-1` before, `2.079e-7` after; right `1.889e-7` both) and not on aarch64 or on x86-64 without AVX (Rosetta), where it passes before and after; `c64` does not take it. Tests: `linalg::gevd::qz_cplx::tests::unblocked_zero_chase_updates_q_columns` gives the unblocked iteration an exact zero at `j == ilo` and below a zero subdiagonal at `c32`/`c64`, `n = 2, 3, 6, 12`, and checks `Q`, `Z` unitary and `A = Q S Z^H`, `B = Q T Z^H` to `10 n eps` (fails before on aarch64: `2.51` against `7.0e-15`); `faer/tests/gevd_infinite_eigenvalue.rs` embeds the issue's bit patterns and checks the left and right residuals at `c32` and `c64` (fails before on x86-64 AVX2 only); before folding a column into the maxima it asserts both residual numerators and `|beta|/(|alpha|+|beta|)` finite and `||u||`, `||v||` and the residual denominator non-zero, with the column number and values in the message, because `f64::max`/`f64::min` drop a `NaN` and the `0/0` of an all-zero left eigenvector column would otherwise read as a perfect residual (an all-zero `ul` passed the test before these asserts and fails after: `column 0: left eigenvector has norm 0e0, alpha = Complex { re: -2.6592135429382324, im: 2.4517483711242676 }, beta = Complex { re: 0.0, im: 0.0 }`). The same code is in upstream faer `main` `5833f4e` at `faer/src/linalg/gevd/qz_cplx/mod.rs:123-127`, from `57e25ac` (generalized evd initial impl); verified by reading that source, not by running it, and no `upstream-reports/` entry is added (the deferred-feedback decision of `AGENTS.md`). | +| Complex QZ: a zero chased down `B`'s diagonal updates the columns of `Q` (#324) | When `hessenberg_to_qz_unblocked` (the unblocked complex QZ path, which the default parameters run for `n < max(15, blocking_threshold)`, i.e. `n < 75`, and which the blocked path also calls at the end) finds `T[(j, j)]` negligible with `j == ilo` or a negligible `H[(j, j - 1)]`, it chases the zero to the bottom with rotations on rows `jch, jch + 1` of `H` and `T`. `Q` must take the adjoint rotation on its columns, as the QZ sweep, the other zero-chase branch and the blocked path do (LAPACK `zhgeqz`: `ZROT(N, Q(1,JCH), 1, Q(1,JCH+1), 1, C, DCONJG(S))`); it rotated rows `jch, jch + 1` of `Q` with `s` instead. `Q` stayed unitary but `A = Q S Z^H` failed (`O(1)`), so every left eigenvector of `gevd_cplx` was wrong while the eigenvalues and the right eigenvectors, which use only `S`, `T` and `Z`, were right. It now rotates columns `jch, jch + 1` with `conj(s)`. | KEEP | tensor4all/faer#324 | Correctness only, no timing and no performance claim: the same number of rotations on a column pair instead of a row pair. Whether the branch runs depends on rounding: the issue's 12 x 12 `c32` pencil takes it on x86-64 with AVX2/FMA (left residual `1.962e-1` before, `2.079e-7` after; right `1.889e-7` both) and not on aarch64 or on x86-64 without AVX (Rosetta), where it passes before and after; `c64` does not take it. Tests: `linalg::gevd::qz_cplx::tests::unblocked_zero_chase_updates_q_columns` gives the unblocked iteration an exact zero at `j == ilo` and below a zero subdiagonal at `c32`/`c64`, `n = 2, 3, 6, 12`, and checks `Q`, `Z` unitary and `A = Q S Z^H`, `B = Q T Z^H` to `10 n eps` (fails before on aarch64: `2.51` against `7.0e-15`); `faer/tests/gevd_infinite_eigenvalue.rs` embeds the issue's bit patterns and checks the left and right residuals at `c32` and `c64` (fails before on x86-64 AVX2 only); before folding a column into the maxima it asserts the left and right residuals, both numerators and `|beta|/(|alpha|+|beta|)` finite and `||u||`, `||v||` and the residual denominator non-zero, with the column number and values in the message, because `f64::max`/`f64::min` drop a `NaN` and the `0/0` of an all-zero eigenvector column would otherwise read as a perfect residual. Measured: with `ul = Mat::::zeros(n, n)` after the `gevd_cplx` call the test passed before these asserts and fails after with `column 0: left eigenvector has norm 0e0, alpha = Complex { re: -2.6592135429382324, im: 2.4517483711242676 }, beta = Complex { re: 0.0, im: 0.0 }`, and the same mutation on `ur` fails with `column 0: right eigenvector has norm 0e0`. The dispatch claim was measured too, with a temporary `eprintln!` at both `hessenberg_to_qz_unblocked` call sites and default `params`: a 20 x 20 pencil prints `n=20 nmin=75 blocking_threshold=75 -> unblocked`, and a 100 x 100 one prints `n=100 nmin=75 -> final unblocked call after the blocked loop`. The same code is in upstream faer `main` `5833f4e` at `faer/src/linalg/gevd/qz_cplx/mod.rs:123-127`, from `57e25ac` (generalized evd initial impl); verified by reading that source, not by running it, and no `upstream-reports/` entry is added (the deferred-feedback decision of `AGENTS.md`). | ## Ground rules diff --git a/faer/tests/gevd_infinite_eigenvalue.rs b/faer/tests/gevd_infinite_eigenvalue.rs index 1b4aa970..971b0fd2 100644 --- a/faer/tests/gevd_infinite_eigenvalue.rs +++ b/faer/tests/gevd_infinite_eigenvalue.rs @@ -228,8 +228,20 @@ fn residuals( den > 0.0, "column {j}: residual denominator {den:e}, ||A|| = {na:e}, ||B|| = {nb:e}, alpha = {al:?}, beta = {be:?}" ); - left = left.max(nlu / (nu * den)); - right = right.max(nrv / (nv * den)); + let (lj, rj) = (nlu / (nu * den), nrv / (nv * den)); + // `f64::max` returns the other operand when one is `NaN`, so a `NaN` + // folded into `left`/`right` would leave it unchanged and read as a + // perfect residual instead of failing. + assert!( + lj.is_finite(), + "column {j}: left residual {lj:e} (numerator {nlu:e}, ||u|| {nu:e}, denominator {den:e})" + ); + assert!( + rj.is_finite(), + "column {j}: right residual {rj:e} (numerator {nrv:e}, ||v|| {nv:e}, denominator {den:e})" + ); + left = left.max(lj); + right = right.max(rj); } (left, right, beta_min) } From de39e61a3eab4c57c1a08d5e77a3ad2681a127a1 Mon Sep 17 00:00:00 2001 From: Hiroshi Shinaoka Date: Mon, 9 Nov 2026 14:45:15 +0900 Subject: [PATCH 4/5] test(gevd): divide the residual in two steps so it cannot overflow The residual was `numerator / (||v|| * denominator)`. For a column whose eigenvector norm is large enough the product overflows to infinity and the quotient is 0, which passes `is_finite` and reads as a perfect residual. Replacing a column's eigenvectors with `k * e0`, `k = 1.125 * (f64::MAX / den)`, left the test green with `nu * den = inf` and the residual recorded as `0e0` against a staged `2.978e-1`. Two divisions `(numerator / ||v||) / denominator` make the same mutation fail at `c32 left residual 2.977833357222273e-1 >= 1.430511474609375e-5`. --- PATCHSET.md | 2 +- faer/tests/gevd_infinite_eigenvalue.rs | 5 ++++- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/PATCHSET.md b/PATCHSET.md index 19298919..8ab43ac2 100644 --- a/PATCHSET.md +++ b/PATCHSET.md @@ -194,7 +194,7 @@ gate has changed and two kinds of adoption are now in the ledger. | `stats` column/row reductions with a reversed output axis (#317) | `col_mean`/`col_varm`, and through them `row_mean`/`row_varm`, normalize a negative row stride with `mat.reverse_rows()` before the kernel runs. Rows are the output axis of a column reduction, and `out` (and, for `varm`, `col_mean`) were not reversed with them. Entry `i` was written to `nrows - 1 - i`, and `varm` subtracted the mean of the mirrored row. All four internal functions (`col_{mean,varm}_{ignore,propagate}`) now reverse `out`, and `col_mean` for `varm`, together with `mat`. Reversing the reduced axis (`reverse_cols` in `col_*`) was already harmless. Public API unchanged. `meanvar.rs` is byte-identical to the upstream sync base `0ecd846`, so upstream faer carries the same defect there; later upstream is not checked, and feedback is deferred. | EXPERIMENT | tensor4all/faer#317 | Correctness only, no timing claim. New tests in `stats::meanvar::tests`: `col_*` on `reverse_rows()` and `row_*` on `reverse_cols()`, mean and varm, `Propagate` and `Ignore`, f32/f64/c32/c64, on column-major (fallback kernel) and row-major (row-major SIMD kernel) storage with one NaN. Each is compared with the same matrix copied to an owned one. The issue's reproducer is pinned as well. Removing the output reversal from any one of the four functions fails both tests. | | Explicit sparse column ordering stops paying for COLAMD (#322) | `factorize_symbolic_lu_with_column_ordering` and `factorize_symbolic_qr_with_column_ordering` build their internal `MemBuffer` from `StackReq::or(colamd::order_scratch::(m, n, A_nnz), ..)`, but COLAMD runs only in the `None` branch, so a caller that brings the ordering still allocates a COLAMD-sized workspace on every symbolic analysis (dyn-stack's `MemBuffer` takes it uninitialized, so the cost is the allocation, not a memset). Both functions include `order_scratch` only when `column_ordering.is_none()` now; the `None` request and the analysis request are unchanged. Both entry points are this fork's own (#239, `3889bd1`), so upstream faer does not carry the defect and no `upstream-reports/` entry is added. | KEEP | tensor4all/faer#322 | Allocation only, no numerical change: the two functions' permutations, symbolic structures and every downstream numeric path are untouched, and no benchmark was run. The issue's reproducer (n = 2000, nnz = 7990, identity ordering, `--release`, counting global allocator) goes from 332712 B largest allocation to 239928 B for simplicial LU and QR: the 332712 B was exactly `colamd::order_scratch`, and 239928 B is the analysis request (4 + 5 n-element arrays, the transpose arrays, and the largest simplicial arm), while the supernodal arm stays at 468080 B, which is the supernodal symbolic workspace rather than the COLAMD one. `faer/tests/sparse_symbolic_colamd_workspace.rs` pins both directions on a 100 x 100 sample: with an explicit ordering the largest allocation stays below `order_scratch`, without one it reaches it; the first fails before the change and both pass after. `cargo test -p t4a-faer -j6` passes except `lq_no_alloc::lq_allocates_nothing_60x300`, which fails identically on `origin/main` (dense LQ, unrelated), plus the process-wide allocation counters in `no_step_alloc` and `triangular_solve::aliased_rhs::staged_solve_is_bitwise_the_in_place_one`, which pass standalone and flake only inside the whole-suite parallel run. | | Complex `abs_max` drops a NaN in the real part (#320) | `ComplexField::simd_abs_max` takes the larger of the two part magnitudes with an unordered `re > im` select, and that select returns `im` when `re` is NaN: a NaN real part read as the finite imaginary magnitude, while a NaN imaginary part did return NaN. The imaginary side already lands on `im`, so the real side takes an explicit `!(re == re)` mask in the generic `Complex` SIMD path, and the scalar fallbacks of `ComplexImpl` / `ComplexImpl` test both parts before the comparison. No in-tree type reaches either path (native `c32` / `c64` go through the real coercion in `norm_max`), so no dispatch, allocation or numerical result on an existing input changes. | KEEP | tensor4all/faer#320 | Correctness; the `faer-traits` unit test `complex_abs_max_propagates_nan` drives `SimdCtx::abs_max` on `Complex` and fails before / passes after (real-part NaN, imaginary-part NaN, both, and three finite controls). No performance claim: two mask ops per SIMD chunk on a path no in-tree type reaches. | -| Complex QZ: a zero chased down `B`'s diagonal updates the columns of `Q` (#324) | When `hessenberg_to_qz_unblocked` (the unblocked complex QZ path, which the default parameters run for `n < max(15, blocking_threshold)`, i.e. `n < 75`, and which the blocked path also calls at the end) finds `T[(j, j)]` negligible with `j == ilo` or a negligible `H[(j, j - 1)]`, it chases the zero to the bottom with rotations on rows `jch, jch + 1` of `H` and `T`. `Q` must take the adjoint rotation on its columns, as the QZ sweep, the other zero-chase branch and the blocked path do (LAPACK `zhgeqz`: `ZROT(N, Q(1,JCH), 1, Q(1,JCH+1), 1, C, DCONJG(S))`); it rotated rows `jch, jch + 1` of `Q` with `s` instead. `Q` stayed unitary but `A = Q S Z^H` failed (`O(1)`), so every left eigenvector of `gevd_cplx` was wrong while the eigenvalues and the right eigenvectors, which use only `S`, `T` and `Z`, were right. It now rotates columns `jch, jch + 1` with `conj(s)`. | KEEP | tensor4all/faer#324 | Correctness only, no timing and no performance claim: the same number of rotations on a column pair instead of a row pair. Whether the branch runs depends on rounding: the issue's 12 x 12 `c32` pencil takes it on x86-64 with AVX2/FMA (left residual `1.962e-1` before, `2.079e-7` after; right `1.889e-7` both) and not on aarch64 or on x86-64 without AVX (Rosetta), where it passes before and after; `c64` does not take it. Tests: `linalg::gevd::qz_cplx::tests::unblocked_zero_chase_updates_q_columns` gives the unblocked iteration an exact zero at `j == ilo` and below a zero subdiagonal at `c32`/`c64`, `n = 2, 3, 6, 12`, and checks `Q`, `Z` unitary and `A = Q S Z^H`, `B = Q T Z^H` to `10 n eps` (fails before on aarch64: `2.51` against `7.0e-15`); `faer/tests/gevd_infinite_eigenvalue.rs` embeds the issue's bit patterns and checks the left and right residuals at `c32` and `c64` (fails before on x86-64 AVX2 only); before folding a column into the maxima it asserts the left and right residuals, both numerators and `|beta|/(|alpha|+|beta|)` finite and `||u||`, `||v||` and the residual denominator non-zero, with the column number and values in the message, because `f64::max`/`f64::min` drop a `NaN` and the `0/0` of an all-zero eigenvector column would otherwise read as a perfect residual. Measured: with `ul = Mat::::zeros(n, n)` after the `gevd_cplx` call the test passed before these asserts and fails after with `column 0: left eigenvector has norm 0e0, alpha = Complex { re: -2.6592135429382324, im: 2.4517483711242676 }, beta = Complex { re: 0.0, im: 0.0 }`, and the same mutation on `ur` fails with `column 0: right eigenvector has norm 0e0`. The dispatch claim was measured too, with a temporary `eprintln!` at both `hessenberg_to_qz_unblocked` call sites and default `params`: a 20 x 20 pencil prints `n=20 nmin=75 blocking_threshold=75 -> unblocked`, and a 100 x 100 one prints `n=100 nmin=75 -> final unblocked call after the blocked loop`. The same code is in upstream faer `main` `5833f4e` at `faer/src/linalg/gevd/qz_cplx/mod.rs:123-127`, from `57e25ac` (generalized evd initial impl); verified by reading that source, not by running it, and no `upstream-reports/` entry is added (the deferred-feedback decision of `AGENTS.md`). | +| Complex QZ: a zero chased down `B`'s diagonal updates the columns of `Q` (#324) | When `hessenberg_to_qz_unblocked` (the unblocked complex QZ path, which the default parameters run for `n < max(15, blocking_threshold)`, i.e. `n < 75`, and which the blocked path also calls at the end) finds `T[(j, j)]` negligible with `j == ilo` or a negligible `H[(j, j - 1)]`, it chases the zero to the bottom with rotations on rows `jch, jch + 1` of `H` and `T`. `Q` must take the adjoint rotation on its columns, as the QZ sweep, the other zero-chase branch and the blocked path do (LAPACK `zhgeqz`: `ZROT(N, Q(1,JCH), 1, Q(1,JCH+1), 1, C, DCONJG(S))`); it rotated rows `jch, jch + 1` of `Q` with `s` instead. `Q` stayed unitary but `A = Q S Z^H` failed (`O(1)`), so every left eigenvector of `gevd_cplx` was wrong while the eigenvalues and the right eigenvectors, which use only `S`, `T` and `Z`, were right. It now rotates columns `jch, jch + 1` with `conj(s)`. | KEEP | tensor4all/faer#324 | Correctness only, no timing and no performance claim: the same number of rotations on a column pair instead of a row pair. Whether the branch runs depends on rounding: the issue's 12 x 12 `c32` pencil takes it on x86-64 with AVX2/FMA (left residual `1.962e-1` before, `2.079e-7` after; right `1.889e-7` both) and not on aarch64 or on x86-64 without AVX (Rosetta), where it passes before and after; `c64` does not take it. Tests: `linalg::gevd::qz_cplx::tests::unblocked_zero_chase_updates_q_columns` gives the unblocked iteration an exact zero at `j == ilo` and below a zero subdiagonal at `c32`/`c64`, `n = 2, 3, 6, 12`, and checks `Q`, `Z` unitary and `A = Q S Z^H`, `B = Q T Z^H` to `10 n eps` (fails before on aarch64: `2.51` against `7.0e-15`); `faer/tests/gevd_infinite_eigenvalue.rs` embeds the issue's bit patterns and checks the left and right residuals at `c32` and `c64` (fails before on x86-64 AVX2 only); before folding a column into the maxima it asserts the left and right residuals, both numerators and `|beta|/(|alpha|+|beta|)` finite and `||u||`, `||v||` and the residual denominator non-zero (each residual as `(numerator / ||u||) / denominator`, because the single-division form overflows the product to infinity for a column with a huge norm and the quotient is then 0), with the column number and values in the message, because `f64::max`/`f64::min` drop a `NaN` and the `0/0` of an all-zero eigenvector column would otherwise read as a perfect residual. Measured: with `ul = Mat::::zeros(n, n)` after the `gevd_cplx` call the test passed before these asserts and fails after with `column 0: left eigenvector has norm 0e0, alpha = Complex { re: -2.6592135429382324, im: 2.4517483711242676 }, beta = Complex { re: 0.0, im: 0.0 }`, and the same mutation on `ur` fails with `column 0: right eigenvector has norm 0e0`; replacing each column's eigenvectors with `k * e0`, `k = 1.125 * (f64::MAX / den)`, passed with the single-division form (`nu * den = inf`, residual recorded as `0e0` against a staged `2.978e-1`) and fails with `c32 left residual 2.977833357222273e-1 >= 1.430511474609375e-5` after it. The dispatch claim was measured too, with a temporary `eprintln!` at both `hessenberg_to_qz_unblocked` call sites and default `params`: a 20 x 20 pencil prints `n=20 nmin=75 blocking_threshold=75 -> unblocked`, and a 100 x 100 one prints `n=100 nmin=75 -> final unblocked call after the blocked loop`. The same code is in upstream faer `main` `5833f4e` at `faer/src/linalg/gevd/qz_cplx/mod.rs:123-127`, from `57e25ac` (generalized evd initial impl); verified by reading that source, not by running it, and no `upstream-reports/` entry is added (the deferred-feedback decision of `AGENTS.md`). | ## Ground rules diff --git a/faer/tests/gevd_infinite_eigenvalue.rs b/faer/tests/gevd_infinite_eigenvalue.rs index 971b0fd2..00cbf68d 100644 --- a/faer/tests/gevd_infinite_eigenvalue.rs +++ b/faer/tests/gevd_infinite_eigenvalue.rs @@ -228,7 +228,10 @@ fn residuals( den > 0.0, "column {j}: residual denominator {den:e}, ||A|| = {na:e}, ||B|| = {nb:e}, alpha = {al:?}, beta = {be:?}" ); - let (lj, rj) = (nlu / (nu * den), nrv / (nv * den)); + // Two divisions, not `nlu / (nu * den)`: the product overflows to + // infinity for a column with a huge norm, and the quotient is then 0, + // which reads as a perfect residual. + let (lj, rj) = ((nlu / nu) / den, (nrv / nv) / den); // `f64::max` returns the other operand when one is `NaN`, so a `NaN` // folded into `left`/`right` would leave it unchanged and read as a // perfect residual instead of failing. From 8e96c95baf3ccaee4e5c3688e994d11bcc4c74d4 Mon Sep 17 00:00:00 2001 From: Hiroshi Shinaoka Date: Mon, 9 Nov 2026 15:58:55 +0900 Subject: [PATCH 5/5] test(gevd): normalize the eigenpair before forming the residual `residuals()` built each column's numerator out of the eigenvector and the homogeneous pair as they come. The residual is invariant under both scalings, so the value was right whenever it survived, but two columns made it read as a perfect residual instead: - `alpha = f64::MAX / 2, beta = 0, u = v = 1e-300 * e0`: the denominator `|beta| ||A|| + |alpha| ||B||` overflows to `inf`, and `x / inf` is `0`; - `alpha = 1e-150, beta = 0, u = v = 1e-200 * e0`: the numerator underflows to `0` before the division. Both passed every guard and both are wrong by `O(1)` once normalized (`2.977833e-1` left, `2.941382e-1` right). Dividing afterwards cannot recover information the numerator never had, so the eigenvector is now divided by its norm and `(alpha, beta)` by `max(|alpha|, |beta|)` before the residual is formed, the two norms, that scale and the denominator are required to be finite and positive, and the residual is one division again because `||u|| = 1` and `|alpha|, |beta| <= 1` make the numerator at most the denominator. The two columns are now the test `residual_oracle_rejects_unnormalized_eigenpairs`, which the unnormalized oracle fails with `left residual 0e0 < 1.43e-5` and `right residual 0e0 < 1.43e-5` for both. The columns of the issue's pencil read the same as before (`c32` 1.772195055005881e-7 / 1.3466182542302985e-7, `c64` 3.838923294839156e-16 / 2.9362584762522423e-16 against `3.864850603564751e-16` / `2.9009055001556393e-16` unnormalized; tolerance 2.7e-14). --- PATCHSET.md | 2 +- faer/tests/gevd_infinite_eigenvalue.rs | 142 ++++++++++++++++++------- 2 files changed, 104 insertions(+), 40 deletions(-) diff --git a/PATCHSET.md b/PATCHSET.md index 8ab43ac2..021c6479 100644 --- a/PATCHSET.md +++ b/PATCHSET.md @@ -194,7 +194,7 @@ gate has changed and two kinds of adoption are now in the ledger. | `stats` column/row reductions with a reversed output axis (#317) | `col_mean`/`col_varm`, and through them `row_mean`/`row_varm`, normalize a negative row stride with `mat.reverse_rows()` before the kernel runs. Rows are the output axis of a column reduction, and `out` (and, for `varm`, `col_mean`) were not reversed with them. Entry `i` was written to `nrows - 1 - i`, and `varm` subtracted the mean of the mirrored row. All four internal functions (`col_{mean,varm}_{ignore,propagate}`) now reverse `out`, and `col_mean` for `varm`, together with `mat`. Reversing the reduced axis (`reverse_cols` in `col_*`) was already harmless. Public API unchanged. `meanvar.rs` is byte-identical to the upstream sync base `0ecd846`, so upstream faer carries the same defect there; later upstream is not checked, and feedback is deferred. | EXPERIMENT | tensor4all/faer#317 | Correctness only, no timing claim. New tests in `stats::meanvar::tests`: `col_*` on `reverse_rows()` and `row_*` on `reverse_cols()`, mean and varm, `Propagate` and `Ignore`, f32/f64/c32/c64, on column-major (fallback kernel) and row-major (row-major SIMD kernel) storage with one NaN. Each is compared with the same matrix copied to an owned one. The issue's reproducer is pinned as well. Removing the output reversal from any one of the four functions fails both tests. | | Explicit sparse column ordering stops paying for COLAMD (#322) | `factorize_symbolic_lu_with_column_ordering` and `factorize_symbolic_qr_with_column_ordering` build their internal `MemBuffer` from `StackReq::or(colamd::order_scratch::(m, n, A_nnz), ..)`, but COLAMD runs only in the `None` branch, so a caller that brings the ordering still allocates a COLAMD-sized workspace on every symbolic analysis (dyn-stack's `MemBuffer` takes it uninitialized, so the cost is the allocation, not a memset). Both functions include `order_scratch` only when `column_ordering.is_none()` now; the `None` request and the analysis request are unchanged. Both entry points are this fork's own (#239, `3889bd1`), so upstream faer does not carry the defect and no `upstream-reports/` entry is added. | KEEP | tensor4all/faer#322 | Allocation only, no numerical change: the two functions' permutations, symbolic structures and every downstream numeric path are untouched, and no benchmark was run. The issue's reproducer (n = 2000, nnz = 7990, identity ordering, `--release`, counting global allocator) goes from 332712 B largest allocation to 239928 B for simplicial LU and QR: the 332712 B was exactly `colamd::order_scratch`, and 239928 B is the analysis request (4 + 5 n-element arrays, the transpose arrays, and the largest simplicial arm), while the supernodal arm stays at 468080 B, which is the supernodal symbolic workspace rather than the COLAMD one. `faer/tests/sparse_symbolic_colamd_workspace.rs` pins both directions on a 100 x 100 sample: with an explicit ordering the largest allocation stays below `order_scratch`, without one it reaches it; the first fails before the change and both pass after. `cargo test -p t4a-faer -j6` passes except `lq_no_alloc::lq_allocates_nothing_60x300`, which fails identically on `origin/main` (dense LQ, unrelated), plus the process-wide allocation counters in `no_step_alloc` and `triangular_solve::aliased_rhs::staged_solve_is_bitwise_the_in_place_one`, which pass standalone and flake only inside the whole-suite parallel run. | | Complex `abs_max` drops a NaN in the real part (#320) | `ComplexField::simd_abs_max` takes the larger of the two part magnitudes with an unordered `re > im` select, and that select returns `im` when `re` is NaN: a NaN real part read as the finite imaginary magnitude, while a NaN imaginary part did return NaN. The imaginary side already lands on `im`, so the real side takes an explicit `!(re == re)` mask in the generic `Complex` SIMD path, and the scalar fallbacks of `ComplexImpl` / `ComplexImpl` test both parts before the comparison. No in-tree type reaches either path (native `c32` / `c64` go through the real coercion in `norm_max`), so no dispatch, allocation or numerical result on an existing input changes. | KEEP | tensor4all/faer#320 | Correctness; the `faer-traits` unit test `complex_abs_max_propagates_nan` drives `SimdCtx::abs_max` on `Complex` and fails before / passes after (real-part NaN, imaginary-part NaN, both, and three finite controls). No performance claim: two mask ops per SIMD chunk on a path no in-tree type reaches. | -| Complex QZ: a zero chased down `B`'s diagonal updates the columns of `Q` (#324) | When `hessenberg_to_qz_unblocked` (the unblocked complex QZ path, which the default parameters run for `n < max(15, blocking_threshold)`, i.e. `n < 75`, and which the blocked path also calls at the end) finds `T[(j, j)]` negligible with `j == ilo` or a negligible `H[(j, j - 1)]`, it chases the zero to the bottom with rotations on rows `jch, jch + 1` of `H` and `T`. `Q` must take the adjoint rotation on its columns, as the QZ sweep, the other zero-chase branch and the blocked path do (LAPACK `zhgeqz`: `ZROT(N, Q(1,JCH), 1, Q(1,JCH+1), 1, C, DCONJG(S))`); it rotated rows `jch, jch + 1` of `Q` with `s` instead. `Q` stayed unitary but `A = Q S Z^H` failed (`O(1)`), so every left eigenvector of `gevd_cplx` was wrong while the eigenvalues and the right eigenvectors, which use only `S`, `T` and `Z`, were right. It now rotates columns `jch, jch + 1` with `conj(s)`. | KEEP | tensor4all/faer#324 | Correctness only, no timing and no performance claim: the same number of rotations on a column pair instead of a row pair. Whether the branch runs depends on rounding: the issue's 12 x 12 `c32` pencil takes it on x86-64 with AVX2/FMA (left residual `1.962e-1` before, `2.079e-7` after; right `1.889e-7` both) and not on aarch64 or on x86-64 without AVX (Rosetta), where it passes before and after; `c64` does not take it. Tests: `linalg::gevd::qz_cplx::tests::unblocked_zero_chase_updates_q_columns` gives the unblocked iteration an exact zero at `j == ilo` and below a zero subdiagonal at `c32`/`c64`, `n = 2, 3, 6, 12`, and checks `Q`, `Z` unitary and `A = Q S Z^H`, `B = Q T Z^H` to `10 n eps` (fails before on aarch64: `2.51` against `7.0e-15`); `faer/tests/gevd_infinite_eigenvalue.rs` embeds the issue's bit patterns and checks the left and right residuals at `c32` and `c64` (fails before on x86-64 AVX2 only); before folding a column into the maxima it asserts the left and right residuals, both numerators and `|beta|/(|alpha|+|beta|)` finite and `||u||`, `||v||` and the residual denominator non-zero (each residual as `(numerator / ||u||) / denominator`, because the single-division form overflows the product to infinity for a column with a huge norm and the quotient is then 0), with the column number and values in the message, because `f64::max`/`f64::min` drop a `NaN` and the `0/0` of an all-zero eigenvector column would otherwise read as a perfect residual. Measured: with `ul = Mat::::zeros(n, n)` after the `gevd_cplx` call the test passed before these asserts and fails after with `column 0: left eigenvector has norm 0e0, alpha = Complex { re: -2.6592135429382324, im: 2.4517483711242676 }, beta = Complex { re: 0.0, im: 0.0 }`, and the same mutation on `ur` fails with `column 0: right eigenvector has norm 0e0`; replacing each column's eigenvectors with `k * e0`, `k = 1.125 * (f64::MAX / den)`, passed with the single-division form (`nu * den = inf`, residual recorded as `0e0` against a staged `2.978e-1`) and fails with `c32 left residual 2.977833357222273e-1 >= 1.430511474609375e-5` after it. The dispatch claim was measured too, with a temporary `eprintln!` at both `hessenberg_to_qz_unblocked` call sites and default `params`: a 20 x 20 pencil prints `n=20 nmin=75 blocking_threshold=75 -> unblocked`, and a 100 x 100 one prints `n=100 nmin=75 -> final unblocked call after the blocked loop`. The same code is in upstream faer `main` `5833f4e` at `faer/src/linalg/gevd/qz_cplx/mod.rs:123-127`, from `57e25ac` (generalized evd initial impl); verified by reading that source, not by running it, and no `upstream-reports/` entry is added (the deferred-feedback decision of `AGENTS.md`). | +| Complex QZ: a zero chased down `B`'s diagonal updates the columns of `Q` (#324) | When `hessenberg_to_qz_unblocked` (the unblocked complex QZ path, which the default parameters run for `n < max(15, blocking_threshold)`, i.e. `n < 75`, and which the blocked path also calls at the end) finds `T[(j, j)]` negligible with `j == ilo` or a negligible `H[(j, j - 1)]`, it chases the zero to the bottom with rotations on rows `jch, jch + 1` of `H` and `T`. `Q` must take the adjoint rotation on its columns, as the QZ sweep, the other zero-chase branch and the blocked path do (LAPACK `zhgeqz`: `ZROT(N, Q(1,JCH), 1, Q(1,JCH+1), 1, C, DCONJG(S))`); it rotated rows `jch, jch + 1` of `Q` with `s` instead. `Q` stayed unitary but `A = Q S Z^H` failed (`O(1)`), so every left eigenvector of `gevd_cplx` was wrong while the eigenvalues and the right eigenvectors, which use only `S`, `T` and `Z`, were right. It now rotates columns `jch, jch + 1` with `conj(s)`. | KEEP | tensor4all/faer#324 | Correctness only, no timing and no performance claim: the same number of rotations on a column pair instead of a row pair. Whether the branch runs depends on rounding: the issue's 12 x 12 `c32` pencil takes it on x86-64 with AVX2/FMA (left residual `1.962e-1` before, `2.079e-7` after; right `1.889e-7` both) and not on aarch64 or on x86-64 without AVX (Rosetta), where it passes before and after; `c64` does not take it. Tests: `linalg::gevd::qz_cplx::tests::unblocked_zero_chase_updates_q_columns` gives the unblocked iteration an exact zero at `j == ilo` and below a zero subdiagonal at `c32`/`c64`, `n = 2, 3, 6, 12`, and checks `Q`, `Z` unitary and `A = Q S Z^H`, `B = Q T Z^H` to `10 n eps` (fails before on aarch64: `2.51` against `7.0e-15`); `faer/tests/gevd_infinite_eigenvalue.rs` embeds the issue's bit patterns and checks the left and right residuals at `c32` and `c64` (fails before on x86-64 AVX2 only); before folding a column into the maxima it normalizes the eigenvector by its norm and `(alpha, beta)` by `max(|alpha|, |beta|)`, which leaves the residual unchanged, and requires the two eigenvector norms, that scale and the residual denominator to be finite and positive; a numerator formed from the values as they come loses the residual instead, and no later division brings it back: `u = v = 1e-200 * e0` with `alpha = 1e-150` underflows it to `0`, and `alpha = f64::MAX / 2` with `beta = 0` sends the denominator to `inf`, whose quotient is `0` as well. It then asserts both residuals finite, with the column number and values in the message, because `f64::max`/`f64::min` drop a `NaN` and the `0/0` of an all-zero eigenvector column would otherwise read as a perfect residual. Measured: with `ul = Mat::::zeros(n, n)` after the `gevd_cplx` call the test passed before these asserts and fails after with `column 0: left eigenvector has norm 0e0, alpha = Complex { re: -2.6592135429382324, im: 2.4517483711242676 }, beta = Complex { re: 0.0, im: 0.0 }`, and the same mutation on `ur` fails with `column 0: right eigenvector has norm 0e0`; replacing each column's eigenvectors with `k * e0`, `k = 1.125 * (f64::MAX / den)`, passed with the single-division form (`nu * den = inf`, residual recorded as `0e0` against a staged `2.978e-1`) and fails with `c32 left residual 2.977833357222273e-1 >= 1.430511474609375e-5` after it; and `residual_oracle_rejects_unnormalized_eigenpairs` hands the oracle the review's two columns directly (`alpha = f64::MAX / 2, beta = 0, u = v = 1e-300 * e0` and `alpha = 1e-150, beta = 0, u = v = 1e-200 * e0`): on the oracle without normalization both read `left residual 0e0 < 1.430511474609375e-5` and `right residual 0e0 < 1.430511474609375e-5`, and with it `2.977833357222273e-1` and `2.941381743992797e-1`. The dispatch claim was measured too, with a temporary `eprintln!` at both `hessenberg_to_qz_unblocked` call sites and default `params`: a 20 x 20 pencil prints `n=20 nmin=75 blocking_threshold=75 -> unblocked`, and a 100 x 100 one prints `n=100 nmin=75 -> final unblocked call after the blocked loop`. The same code is in upstream faer `main` `5833f4e` at `faer/src/linalg/gevd/qz_cplx/mod.rs:123-127`, from `57e25ac` (generalized evd initial impl); verified by reading that source, not by running it, and no `upstream-reports/` entry is added (the deferred-feedback decision of `AGENTS.md`). | ## Ground rules diff --git a/faer/tests/gevd_infinite_eigenvalue.rs b/faer/tests/gevd_infinite_eigenvalue.rs index 00cbf68d..5b4a5fb2 100644 --- a/faer/tests/gevd_infinite_eigenvalue.rs +++ b/faer/tests/gevd_infinite_eigenvalue.rs @@ -21,7 +21,8 @@ //! The test checks, for every column, the relative residuals //! `||beta u^H A - alpha u^H B|| / (||u|| (|beta| ||A|| + |alpha| ||B||))` and //! `||beta A v - alpha B v|| / (||v|| (|beta| ||A|| + |alpha| ||B||))`, -//! computed in `c64`, against `10 * n * eps` of the scalar type. A `NaN` +//! computed in `c64` against `10 * n * eps` of the scalar type, with `u`, `v` +//! and `(alpha, beta)` normalized before the numerator is formed. A `NaN` //! residual or `beta` fails the test instead of being dropped by `f64::max`. use faer::diag::Diag; use faer::dyn_stack::{MemBuffer, MemStack}; @@ -154,6 +155,76 @@ fn to_c64(m: &Mat) -> Mat { }) } +/// relative left and right residuals of one eigenpair, `j` is only for messages. +/// +/// The eigenvector and the homogeneous pair are normalized here, before the +/// residual is formed. The residual is invariant under both scalings, but a +/// numerator built from the values as they come can lose it entirely, and no +/// later division brings it back: `u = 1e-200 * e0` with `alpha = 1e-150` +/// underflows the numerator to `0`, and `alpha = f64::MAX / 2` with `beta = 0` +/// sends the denominator to `inf`, whose quotient is `0` as well. Both pass for +/// a column that is wrong by `O(1)`. +fn column_residual( + a: &Mat, + b: &Mat, + al: c64, + be: c64, + u: &Mat, + v: &Mat, + j: usize, +) -> (f64, f64) { + let (na, nb) = (a.norm_l2(), b.norm_l2()); + let (nu, nv) = (u.norm_l2(), v.norm_l2()); + // the representative scale of the homogeneous pair: at least one of the + // normalized components is `1`, and neither can overflow + let sc = al.norm().max(be.norm()); + assert!( + nu.is_finite() && nu > 0.0, + "column {j}: left eigenvector has norm {nu:e}, alpha = {al:?}, beta = {be:?}" + ); + assert!( + nv.is_finite() && nv > 0.0, + "column {j}: right eigenvector has norm {nv:e}, alpha = {al:?}, beta = {be:?}" + ); + assert!( + sc.is_finite() && sc > 0.0, + "column {j}: eigenvalue scale {sc:e}, alpha = {al:?}, beta = {be:?}" + ); + let (al, be) = (al / sc, be / sc); + let (u, v) = (u / nu, v / nv); + let lu = u.adjoint() * a * Scale(be) - u.adjoint() * b * Scale(al); + let rv = a * &v * Scale(be) - b * &v * Scale(al); + let den = be.norm() * na + al.norm() * nb; + assert!( + den.is_finite() && den > 0.0, + "column {j}: residual denominator {den:e}, ||A|| = {na:e}, ||B|| = {nb:e}, alpha = {al:?}, beta = {be:?}" + ); + let (nlu, nrv) = (lu.norm_l2(), rv.norm_l2()); + assert!( + nlu.is_finite(), + "column {j}: left residual numerator {nlu:e}, alpha = {al:?}, beta = {be:?}" + ); + assert!( + nrv.is_finite(), + "column {j}: right residual numerator {nrv:e}, alpha = {al:?}, beta = {be:?}" + ); + // one division is enough: `u` and `v` are unit and `|alpha|, |beta| <= 1`, so + // the numerator is at most `den` and neither side can overflow + let (lj, rj) = (nlu / den, nrv / den); + // `f64::max` returns the other operand when one is `NaN`, so a `NaN` + // folded into `left`/`right` would leave it unchanged and read as a + // perfect residual instead of failing. + assert!( + lj.is_finite(), + "column {j}: left residual {lj:e} (numerator {nlu:e}, denominator {den:e})" + ); + assert!( + rj.is_finite(), + "column {j}: right residual {rj:e} (numerator {nrv:e}, denominator {den:e})" + ); + (lj, rj) +} + /// returns the largest left and right relative residuals over all columns, and /// the smallest `|beta| / (|alpha| + |beta|)` fn residuals( @@ -189,7 +260,6 @@ fn residuals( .unwrap(); let a = Mat::from_fn(n, n, |i, j| up(&a[(i, j)])); let b = Mat::from_fn(n, n, |i, j| up(&b[(i, j)])); - let (na, nb) = (a.norm_l2(), b.norm_l2()); let (mut left, mut right, mut beta_min) = (0.0f64, 0.0f64, f64::INFINITY); for j in 0..n { let (al, be) = (up(&s[j]), up(&beta[j])); @@ -203,46 +273,10 @@ fn residuals( beta_min = beta_min.min(ratio); let u = Mat::from_fn(n, 1, |i, _| up(&ul[(i, j)])); let v = Mat::from_fn(n, 1, |i, _| up(&ur[(i, j)])); - let lu = u.adjoint() * &a * Scale(be) - u.adjoint() * &b * Scale(al); - let rv = &a * &v * Scale(be) - &b * &v * Scale(al); - let den = be.norm() * na + al.norm() * nb; - let (nlu, nrv) = (lu.norm_l2(), rv.norm_l2()); - let (nu, nv) = (u.norm_l2(), v.norm_l2()); - assert!( - nlu.is_finite(), - "column {j}: left residual numerator {nlu:e}, alpha = {al:?}, beta = {be:?}" - ); - assert!( - nrv.is_finite(), - "column {j}: right residual numerator {nrv:e}, alpha = {al:?}, beta = {be:?}" - ); - assert!( - nu > 0.0, - "column {j}: left eigenvector has norm {nu:e}, alpha = {al:?}, beta = {be:?}" - ); - assert!( - nv > 0.0, - "column {j}: right eigenvector has norm {nv:e}, alpha = {al:?}, beta = {be:?}" - ); - assert!( - den > 0.0, - "column {j}: residual denominator {den:e}, ||A|| = {na:e}, ||B|| = {nb:e}, alpha = {al:?}, beta = {be:?}" - ); - // Two divisions, not `nlu / (nu * den)`: the product overflows to - // infinity for a column with a huge norm, and the quotient is then 0, - // which reads as a perfect residual. - let (lj, rj) = ((nlu / nu) / den, (nrv / nv) / den); + let (lj, rj) = column_residual(&a, &b, al, be, &u, &v, j); // `f64::max` returns the other operand when one is `NaN`, so a `NaN` // folded into `left`/`right` would leave it unchanged and read as a // perfect residual instead of failing. - assert!( - lj.is_finite(), - "column {j}: left residual {lj:e} (numerator {nlu:e}, ||u|| {nu:e}, denominator {den:e})" - ); - assert!( - rj.is_finite(), - "column {j}: right residual {rj:e} (numerator {nrv:e}, ||v|| {nv:e}, denominator {den:e})" - ); left = left.max(lj); right = right.max(rj); } @@ -270,3 +304,33 @@ fn left_eigenvectors_of_issue_324_pencil() { let (a, b) = (from_bits(&A_BITS), from_bits(&B_BITS)); check(to_c64(&a), to_c64(&b)); } + +/// the oracle itself, on eigenpairs that are wrong by `O(1)` once normalized. +/// +/// A column with `alpha = f64::MAX / 2, beta = 0` sends the denominator to +/// `inf`, and one with `alpha = 1e-150, u = v = 1e-200 * e0` underflows the +/// numerator to `0`; both leave `0` as the residual, which passes. The residual +/// is invariant under both scalings, so normalizing the eigenvector and the +/// homogeneous pair first loses nothing, and keeps the numerator from being +/// formed out of numbers that cannot carry it. +#[test] +fn residual_oracle_rejects_unnormalized_eigenpairs() { + let (a, b) = (to_c64(&from_bits(&A_BITS)), to_c64(&from_bits(&B_BITS))); + let n = a.nrows(); + let e0 = |k: f64| Mat::from_fn(n, 1, |i, _| c64::new(if i == 0 { k } else { 0.0 }, 0.0)); + let tol = 10.0 * n as f64 * f32::EPSILON as f64; + let mut bad = Vec::new(); + for (al, k) in [(f64::MAX / 2.0, 1e-300), (1e-150, 1e-200)] { + let (al, be) = (c64::new(al, 0.0), c64::new(0.0, 0.0)); + let (u, v) = (e0(k), e0(k)); + let (l, r) = column_residual(&a, &b, al, be, &u, &v, 0); + for (side, res) in [("left", l), ("right", r)] { + if res < tol { + bad.push(format!( + "{side} residual {res:e} < {tol:e} for alpha = {al:?}, beta = {be:?}, u = v = {k:e} * e0" + )); + } + } + } + assert!(bad.is_empty(), "{}", bad.join("\n")); +}