diff --git a/CHANGELOG.md b/CHANGELOG.md index 26b525a2..010c6ee3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -22,6 +22,12 @@ 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'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). + - `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 6a28e1c0..021c6479 100644 --- a/PATCHSET.md +++ b/PATCHSET.md @@ -194,6 +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 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/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..5b4a5fb2 --- /dev/null +++ b/faer/tests/gevd_infinite_eigenvalue.rs @@ -0,0 +1,336 @@ +//! 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'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 +//! 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, 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}; +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) + }) +} + +/// 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( + 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 (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])); + // `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 (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. + left = left.max(lj); + right = right.max(rj); + } + (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)); +} + +/// 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")); +}