From 5372d9bd985d2bc7daccdb29749de04cd31343a4 Mon Sep 17 00:00:00 2001 From: Hiroshi Shinaoka Date: Thu, 29 Oct 2026 20:18:34 +0900 Subject: [PATCH] evd: the real Schur deflation-window sort always ends (#144) aggressive_early_deflation sorted the window with a bubble sort ported from LAPACK dlaqr3, but over the deflated part (dlaqr3's initial I = NS + 1 only sets the first pass's end) and with each next window set from the last swap's position. schur_swap reports success even when schur_swap_22 rejects an ill-conditioned swap, so the same pair was compared forever: Mat::eigen of the f32 matrix J_n + 1e-6 G scaled to 2^-100 did not return at n = 128 and 257. The sort now follows dlaqr3: it sorts the undeflated part infqr..ns, ends each pass at kend, one before where the previous pass stopped, and swaps with schur_swap_checked (the dlaexc ports for 1 x 2 and 2 x 1), leaving a rejected pair in place. Upstream faer 0.24.4 has the same code; the upstream report is faer-029. Closes #144 Co-Authored-By: Claude Opus 5.5 --- PATCHSET.md | 1 + faer/src/linalg/evd/schur/real_schur.rs | 181 ++++- upstream-reports/README.md | 1 + .../Cargo.lock | 761 ++++++++++++++++++ .../Cargo.toml | 13 + .../README.md | 52 ++ .../src/main.rs | 45 ++ 7 files changed, 1013 insertions(+), 41 deletions(-) create mode 100644 upstream-reports/faer-029-real-schur-aed-sort-hang/Cargo.lock create mode 100644 upstream-reports/faer-029-real-schur-aed-sort-hang/Cargo.toml create mode 100644 upstream-reports/faer-029-real-schur-aed-sort-hang/README.md create mode 100644 upstream-reports/faer-029-real-schur-aed-sort-hang/src/main.rs diff --git a/PATCHSET.md b/PATCHSET.md index b2cb793c..79ff781d 100644 --- a/PATCHSET.md +++ b/PATCHSET.md @@ -133,6 +133,7 @@ gate has changed and two kinds of adoption are now in the ledger. | Team dispatch cost read from threadteam | `utils::cost::team_section_cost` and `MEASURED_MAX_WIDTH` (faer's own linear extrapolation above width 8, #97) are deleted: the extrapolation lives in threadteam (`WidthCosts` 0.1.1, linear above a row's largest anchor; 0.1.2, a measured Zen2 team row for `x86_64-avx2` at widths 2 to 64). `ParShape::dispatch_cost` reads `DispatchCosts::cost(Team, w)`; `par_pays_time`, `team::kernel::kernel_width` and the tiled team product's floor ask `TeamExec::dispatch_cost(w)`. `kernel_width` searches from the full width down instead of stopping at the first narrow width that does not pay (measured rows fall with the width). Requirement `threadteam = "0.1.2"`. Decisions change at widths 7 and 8 for rows whose largest anchor is 6, at width 8 on `x86_64-avx512` (its 0.98 us anchor) and at every width on `x86_64-avx2`. | EXPERIMENT (pending the maintainer's acceptance decision) | tensor4all/threadteam-rs#6, tensor4all/faer#91 | No timing (gate inputs only; acceptance run left to the maintainer). Tests: `utils::cost::tests::the_team_cost_is_threadteams_table_at_every_width`, `utils::team::tests::kernel_width_finds_the_widest_paying_width_on_a_falling_row`. Record: [worklog](docs/worklogs/2026-10-29-team-cost-from-threadteam.md). | | Every region at the widest width that pays | The full-width gates (`par_pays`, `ParThreshold::region_pays`, `ParShape::pays_full` at the call sites, the block apply's `n_tasks`) chose between the token's whole degree and one thread, so on a wide team a region that pays at 8 or 16 ran serially. Every region now asks `utils::cost::par_width` / `ParThreshold::region_width` for the widest width that pays (every width searched, falling rows included; the GEMM kernel's `kernel_width` is the same rule) and runs on a sub-team of that width (`utils::thread::with_width`, a split into one child: no worker dispatch). Splits whose task count is capped (triangular solve jobs of 64 columns, block apply columns, `upgrade_householder_factor`'s rows per participant) search inside their cap. The factorizations that latched `par = Par::Seq` for the rest of the call (Hessenberg, tridiagonal, unblocked bidiagonal, full pivoting LU, Bunch-Kaufman, column pivoting QR) decide per step; `Never` still makes the whole call sequential; the post-loop compact-factor upgrades of the Hessenberg and tridiagonal reductions now ask their own gate (they ran under the last step's latch, i.e. sequentially). The blocked bidiagonal panel gate is a floor (sequential or the whole `par`), so the panel's products choose their own widths. The Bunch-Kaufman blocked swap loop, split over the whole degree without a gate, now has one (row swaps at the LU swap rate). `Fixed` keeps its rules (the whole degree above the floor). `par_pays`, `region_pays` and the Seq-only GEMM predicate are removed; `par-trace` records the chosen width. | EXPERIMENT (pending the maintainer's acceptance run on broccoli / evo) | tensor4all/faer#91, tensor4all/faer#95 | No timing here (performance evidence is the monitor's: broccoli 1/8/16/32/64T, evo). Tests: `utils::team::tests::{par_width_takes_the_widest_width_that_pays, with_width_runs_on_a_sub_team_of_that_width, the_families_narrow_their_regions_and_match_sequential}`, the rewritten `region_rules_follow_the_serial_time`, `fixed_and_never_keep_their_meaning`, `par_width_reads_the_degree_of_the_token`. Record: [worklog](docs/worklogs/2026-10-29-width-per-region.md). | | No heap allocation inside a parallel region (#111) | Every region on the team route allocated on the heap. `for_each_child`, behind every `for_each` and `for_each_with_par`, collected all the items into a `Vec`, one entry per column for `par_col_iter_mut`, and then a slot `Vec`. The full-pivoting LU step and both Bunch-Kaufman full-pivoting updates built an `nthreads`-long vector of pivot candidates per step, and `join_sized` boxed its two closures. The allocations were suspected of causing #109's floor-level loss (`full_piv_lu f64 128` at 8T on evo, +6.6 to +10.4%). The evo re-measurement of this change shows they were not: the cell is at 1.0025 of `main` with them gone. The change stands on its own as a cleanup of the mechanism. Now faer's dispatch and these workspaces allocate nothing. An iterator adapter that allocates in its own `split_at`, such as rayon's `map_with`, still does. `for_each_child` splits the iterator's producer on the leader into the blocks (`split_at`, remainder-relative) and each block's leader iterates its piece. Up to 64 blocks, with a producer of at most 256 bytes, are dealt from an inline array in one fork; wider or heavier splits halve recursively, one section per level. The pivot candidates come from the call's scratch: `lu_in_place_scratch` and `cholesky_in_place_scratch` grow by `par.degree()` entries. `join_sized` takes two typed slots. | KEEP (evo: no regression; `full_piv_lu f64 128` 8T at 1.0025 of `main`) | tensor4all/faer#111, tensor4all/faer#95 | Test `faer/tests/no_step_alloc.rs`: its own binary with a process-wide counting allocator. On teams of 4, 8 and 70 (the last takes the recursive path), with the regions forced parallel and a warm-up first, it counts zero allocations for full-pivoting LU 96, Bunch-Kaufman with full pivoting 96 and `join_sized`, and checks that each one dispatched. On `origin/main` the same test counts 190 allocations in one full-pivoting LU 96 call at width 4. Record: [worklog](docs/worklogs/2026-10-29-no-region-allocation.md). | +| Real Schur: the deflation-window sort always ends (#144) | `aggressive_early_deflation` sorts the window's eigenvalues with a bubble sort ported from LAPACK `dlaqr3`, but it sorted the deflated part (`ns..jw`, from `dlaqr3`'s initial `I = NS + 1`, which only sets the first pass's end) and set each next window from the last swap's position. `schur_swap` reports success even when `schur_swap_22` rejects an ill-conditioned 2 x 2 swap, so the pass kept recording the same position and the sort never ended: `Mat::eigen` of the `f32` matrix `J_n + 1e-6 G` scaled to `2^-100` hangs at n = 128 and 257. The sort now follows `dlaqr3`: it sorts the undeflated part `infqr..ns`, ends each pass at `kend`, one before where the previous pass stopped, and swaps with `schur_swap_checked` (the `dlaexc` ports for 1 x 2 and 2 x 1), leaving a rejected pair in place. Upstream has the same code: faer 0.24.4 `src/linalg/evd/schur/real_schur.rs:921-969`. Not reported upstream. | KEEP | tensor4all/faer#144 | Correctness: the sort ends after at most `ns` passes. Test `tests::aed_sort_terminates_when_a_swap_fails` (the reproducer at n = 128 and 257 through `hessenberg_schur_real` under a 120 s timeout; it times out on `origin/main`; residual within LAPACK's `smlnum` at `2^-100` and within `10 n eps |H|` at `2^-80`, orthogonality `10 n eps`). Upstream report: [`upstream-reports/faer-029-real-schur-aed-sort-hang`](upstream-reports/faer-029-real-schur-aed-sort-hang/). | ## Ground rules diff --git a/faer/src/linalg/evd/schur/real_schur.rs b/faer/src/linalg/evd/schur/real_schur.rs index 55adc493..a2bf7936 100644 --- a/faer/src/linalg/evd/schur/real_schur.rs +++ b/faer/src/linalg/evd/schur/real_schur.rs @@ -1210,56 +1210,52 @@ fn aggressive_early_deflation( ns -= infqr; return (ns, nd); } + // Sort the undeflated eigenvalues `tw[infqr..ns]` by decreasing magnitude, + // as LAPACK `dlaqr3` does. Each pass ends at `kend`, one before where the + // previous pass stopped, so the window shrinks by at least one per pass and + // the sort ends even when a swap leaves a pair out of order. + let ev = |tw: &MatMut<'_, T>, i: usize, n: usize| -> T { + if n == 1 { + tw[(i, i)].abs() + } else { + tw[(i, i)].abs() + + tw[(i + 1, i)].abs().sqrt() * tw[(i, i + 1)].abs().sqrt() + } + }; + // Size of the block starting at `i` when the window ends at `kend`. + let block = |tw: &MatMut<'_, T>, i: usize, kend: usize| -> usize { + if i >= kend || tw[(i + 1, i)] == zero() { 1 } else { 2 } + }; let mut sorted = false; - let mut sorting_window_size = jw as isize; - while !sorted { + let mut i = ns; + while !sorted && i > infqr { sorted = true; - let mut ilst = 0isize; - let mut i1 = ns; - while i1 as isize + 1 < sorting_window_size { - let mut n1 = 1; - if tw[(i1 + 1, i1)] != zero() { - n1 = 2; - } - if i1 + n1 == jw { - ilst -= n1 as isize; - break; - } - let i2 = i1 + n1; - let mut n2 = 1; - if i2 + 1 < jw && (tw[(i2 + 1, i2)] != zero()) { - n2 = 2; - } - let (ev1, ev2); - if n1 == 1 { - ev1 = tw[(i1, i1)].abs(); - } else { - ev1 = tw[(i1, i1)].abs() - + tw[(i1 + 1, i1)].abs().sqrt() - * tw[(i1, i1 + 1)].abs().sqrt(); - } - if n2 == 1 { - ev2 = tw[(i2, i2)].abs(); - } else { - ev2 = tw[(i2, i2)].abs() - + tw[(i2 + 1, i2)].abs().sqrt() - * tw[(i2, i2 + 1)].abs().sqrt(); - } - if ev1 >= ev2 { - i1 = i2; + let kend = i - 1; + i = infqr; + let mut k = i + block(&tw, i, kend); + while k <= kend { + let n1 = k - i; + let n2 = block(&tw, k, kend); + if ev(&tw, i, n1) >= ev(&tw, k, n2) { + i = k; } else { sorted = false; - let ierr = - schur_swap(tw.rb_mut(), Some(v.rb_mut()), i1, n1, n2); + let ierr = schur_swap_checked( + tw.rb_mut(), + Some(v.rb_mut()), + i, + n1, + n2, + true, + ); if ierr == 0 { - i1 += n2; + i += n2; } else { - i1 = i2; + i = k; } - ilst = i1 as isize; } + k = i + block(&tw, i, kend); } - sorting_window_size = ilst; } let mut i = 0; while i < jw { @@ -3273,4 +3269,107 @@ mod tests { } } } + /// tensor4all/faer#144: the deflation-window sort must finish when a swap + /// is rejected. `J_n + 1e-6 G` (Jordan block with eigenvalue 1), divided by + /// its largest entry and scaled by `2^-100`, made the sort compare the same + /// pair forever in `f32` at n = 128 and 257. Each run is bounded by a + /// timeout and checked for a similarity: at `2^-100` the input is below + /// LAPACK's negligibility level `smlnum = safmin * n / eps`, so the residual + /// is held to that absolute level; the same matrix at `2^-80` is held to + /// the usual relative bound. + #[test] + fn aed_sort_terminates_when_a_swap_fails() { + use crate::linalg::evd::hessenberg::{ + hessenberg_in_place, hessenberg_in_place_scratch, + }; + use crate::linalg::evd::{ + hessenberg_schur_real, hessenberg_schur_scratch, + }; + use rand::prelude::*; + use rand_distr::StandardNormal; + use std::sync::mpsc; + use std::time::Duration; + + /// Returns `max |h - q t q^T|`, `max |h|` and `max |q^T q - I|`. + fn run(n: usize, exponent: i32) -> (f32, f32, f32) { + let rng = &mut StdRng::seed_from_u64(112); + let mut a = Mat::::from_fn(n, n, |i, j| { + let g: f32 = rng.sample(StandardNormal); + let jordan = if i == j || i + 1 == j { 1.0 } else { 0.0 }; + jordan + 1e-6 * g + }); + let scale = 2.0f32.powi(exponent) / a.norm_max(); + a *= Scale(scale); + let mut householder = Mat::::zeros(1, n - 1); + hessenberg_in_place( + a.as_mut(), + householder.as_mut(), + Par::Seq, + MemStack::new(&mut MemBuffer::new( + hessenberg_in_place_scratch::( + n, + 1, + ParShape::Seq, + Default::default(), + ), + )), + Default::default(), + ); + let h = Mat::::from_fn(n, n, |i, j| { + if i > j + 1 { 0.0 } else { a[(i, j)] } + }); + let mut t = h.clone(); + let mut q = Mat::::identity(n, n); + let mut wr = Col::::zeros(n); + let mut wi = Col::::zeros(n); + hessenberg_schur_real( + t.as_mut(), + Some(q.as_mut()), + wr.as_mut(), + wi.as_mut(), + Par::Seq, + MemStack::new(&mut MemBuffer::new( + hessenberg_schur_scratch::( + n, + true, + ParShape::Seq, + Default::default(), + ), + )), + Default::default(), + ) + .unwrap(); + assert!(t.is_all_finite() && q.is_all_finite()); + let residual = (&h - &q * &t * q.transpose()).norm_max(); + let orthogonality = + (q.transpose() * &q - Mat::::identity(n, n)).norm_max(); + (residual, h.norm_max(), orthogonality) + } + + for n in [128, 257] { + let n_eps = n as f32 * f32::EPSILON; + let smlnum = f32::MIN_POSITIVE * (n as f32 / f32::EPSILON); + for exponent in [-100, -80] { + let (send, recv) = mpsc::channel(); + std::thread::spawn(move || { + send.send(run(n, exponent)).unwrap() + }); + let (residual, norm, orthogonality) = recv + .recv_timeout(Duration::from_secs(120)) + .unwrap_or_else(|_| { + panic!("n = {n}, 2^{exponent}: real Schur did not finish") + }); + let bound = if exponent == -100 { + smlnum + } else { + 10.0 * n_eps * norm + }; + assert!(residual <= bound, "n = {n}, 2^{exponent}: residual {residual}, bound {bound}"); + assert!( + orthogonality <= 10.0 * n_eps, + "n = {n}, 2^{exponent}: orthogonality {orthogonality}" + ); + } + } + } } diff --git a/upstream-reports/README.md b/upstream-reports/README.md index 75eed75a..41c9f798 100644 --- a/upstream-reports/README.md +++ b/upstream-reports/README.md @@ -42,6 +42,7 @@ The program prints the expected and the actual value, and exits with status 1 wh | [faer-026](faer-026-svd-qr-two-by-two-no-convergence/) | The SVD returns `NoConvergence` on a 12 x 10 complex matrix | `faer` 0.24.4 | no convergence | tensor4all/faer#84 (`f83ffb0`) | not yet | 2026-10-29 | | [faer-027](faer-027-qr-near-dependent-column-skip/) | Unpivoted QR skips a nearly dependent column of a tall matrix | `faer` 0.24.4 | inaccurate value | tensor4all/faer#1 (`2e6ea52`, `d16929b`) | not yet | 2026-10-29 | | [faer-028](faer-028-sparse-lu-numeric-singular/) | Sparse LU factors a numerically singular matrix without an error | `faer` 0.24.4 | NaN | tensor4all/faer#132 | not yet | 2026-10-29 | +| [faer-029](faer-029-real-schur-aed-sort-hang/) | Real Schur QR does not return when the deflation-window sort cannot swap a pair | `faer` 0.24.4 | hang | tensor4all/faer#144 | not yet | 2026-10-29 | | [gemm-x86-001](gemm-x86-001-diag-negative-stride/) | A product with a diagonal gives a wrong value when the left matrix is stored backwards | `private-gemm-x86` 0.1.22 | wrong value | tensor4all/t4a-private-gemm-x86#5 (`84ca6db`) | not yet | 2026-10-29 | "Fixed upstream" was checked against upstream `main`: `faer` at `0ecd846` (2026-10-22) and `private-gemm-x86` at `71872b9` (its 0.1.22 release). diff --git a/upstream-reports/faer-029-real-schur-aed-sort-hang/Cargo.lock b/upstream-reports/faer-029-real-schur-aed-sort-hang/Cargo.lock new file mode 100644 index 00000000..73a5be5a --- /dev/null +++ b/upstream-reports/faer-029-real-schur-aed-sort-hang/Cargo.lock @@ -0,0 +1,761 @@ +# This file is automatically @generated by Cargo. +# It is not intended for manual editing. +version = 4 + +[[package]] +name = "autocfg" +version = "1.5.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f2032f911046de80f0a198e0901378627c33f59ea0ac00e363d481118bd70a53" + +[[package]] +name = "bitflags" +version = "2.13.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3ded4057c258ba199e2d26386d3af3780957ecaee6c4ef4041c6b4b8b97c0b06" + +[[package]] +name = "bytemuck" +version = "1.25.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "95832e849adfb21180ccb6826a99da14e5d266ae5c2e668e1602cf234f153797" +dependencies = [ + "bytemuck_derive", +] + +[[package]] +name = "bytemuck_derive" +version = "1.12.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6a1f896587b6f2c069c73d2f0913e2d590c3990285cd2f0b6aa02b786b4c679c" +dependencies = [ + "proc-macro2", + "quote", + "syn 3.0.9", +] + +[[package]] +name = "byteorder" +version = "1.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1fd0f2584146f6f2ef48085050886acf353beff7305ebd1ae69500e27c67f64b" + +[[package]] +name = "cfg-if" +version = "1.0.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4e7648175b45a9a48536d676f68d918270699102aa8dab5496df06904c914600" + +[[package]] +name = "crunchy" +version = "0.2.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "460fbee9c2c2f33933d720630a6a0bac33ba7053db5344fac858d4b8952d77d5" + +[[package]] +name = "defer" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "930c7171c8df9fb1782bdf9b918ed9ed2d33d1d22300abb754f9085bc48bf8e8" + +[[package]] +name = "dyn-stack" +version = "0.13.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1c4713e43e2886ba72b8271aa66c93d722116acf7a75555cce11dcde84388fe8" +dependencies = [ + "bytemuck", + "dyn-stack-macros", +] + +[[package]] +name = "dyn-stack-macros" +version = "0.1.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e1d926b4d407d372f141f93bb444696142c29d32962ccbd3531117cf3aa0bfa9" + +[[package]] +name = "enum-as-inner" +version = "0.6.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a1e6a265c649f3f5979b601d26f1d05ada116434c87741c9493cb56218f76cbc" +dependencies = [ + "heck", + "proc-macro2", + "quote", + "syn 2.0.119", +] + +[[package]] +name = "equator" +version = "0.2.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c35da53b5a021d2484a7cc49b2ac7f2d840f8236a286f84202369bd338d761ea" +dependencies = [ + "equator-macro 0.2.1", +] + +[[package]] +name = "equator" +version = "0.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "02da895aab06bbebefb6b2595f6d637b18c9ff629b4cd840965bb3164e4194b0" +dependencies = [ + "equator-macro 0.6.0", +] + +[[package]] +name = "equator-macro" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3bf679796c0322556351f287a51b49e48f7c4986e727b5dd78c972d30e2e16cc" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.119", +] + +[[package]] +name = "equator-macro" +version = "0.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2b14b339eb76d07f052cdbad76ca7c1310e56173a138095d3bf42a23c06ef5d8" + +[[package]] +name = "faer" +version = "0.24.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5ab6df3dd147fe8d702a288b95bcd8fcc499ab572fc80da6828f60cd4d524d67" +dependencies = [ + "bytemuck", + "dyn-stack", + "equator 0.6.0", + "faer-traits", + "gemm", + "generativity", + "libm", + "nano-gemm", + "num-complex", + "num-traits", + "private-gemm-x86", + "pulp", + "reborrow", +] + +[[package]] +name = "faer-029-real-schur-aed-sort-hang" +version = "0.0.0" +dependencies = [ + "faer", + "rand", + "rand_distr", +] + +[[package]] +name = "faer-traits" +version = "0.24.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b87d23ed7ab1f26c0cba0e5b9e061a796fbb7dc170fa8bee6970055a1308bb0f" +dependencies = [ + "bytemuck", + "dyn-stack", + "generativity", + "libm", + "num-complex", + "num-traits", + "pulp", + "qd", + "reborrow", +] + +[[package]] +name = "gemm" +version = "0.19.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "aa0673db364b12263d103b68337a68fbecc541d6f6b61ba72fe438654709eacb" +dependencies = [ + "dyn-stack", + "gemm-c32", + "gemm-c64", + "gemm-common", + "gemm-f16", + "gemm-f32", + "gemm-f64", + "num-complex", + "num-traits", + "paste", + "raw-cpuid", + "seq-macro", +] + +[[package]] +name = "gemm-c32" +version = "0.19.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "086936dbdcb99e37aad81d320f98f670e53c1e55a98bee70573e83f95beb128c" +dependencies = [ + "dyn-stack", + "gemm-common", + "num-complex", + "num-traits", + "paste", + "raw-cpuid", + "seq-macro", +] + +[[package]] +name = "gemm-c64" +version = "0.19.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "20c8aeeeec425959bda4d9827664029ba1501a90a0d1e6228e48bef741db3a3f" +dependencies = [ + "dyn-stack", + "gemm-common", + "num-complex", + "num-traits", + "paste", + "raw-cpuid", + "seq-macro", +] + +[[package]] +name = "gemm-common" +version = "0.19.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "88027625910cc9b1085aaaa1c4bc46bb3a36aad323452b33c25b5e4e7c8e2a3e" +dependencies = [ + "bytemuck", + "dyn-stack", + "half", + "libm", + "num-complex", + "num-traits", + "once_cell", + "paste", + "pulp", + "raw-cpuid", + "seq-macro", + "sysctl", +] + +[[package]] +name = "gemm-f16" +version = "0.19.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e3df7a55202e6cd6739d82ae3399c8e0c7e1402859b30e4cb780e61525d9486e" +dependencies = [ + "dyn-stack", + "gemm-common", + "gemm-f32", + "half", + "num-complex", + "num-traits", + "paste", + "raw-cpuid", + "seq-macro", +] + +[[package]] +name = "gemm-f32" +version = "0.19.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "02e0b8c9da1fbec6e3e3ab2ce6bc259ef18eb5f6f0d3e4edf54b75f9fd41a81c" +dependencies = [ + "dyn-stack", + "gemm-common", + "num-complex", + "num-traits", + "paste", + "raw-cpuid", + "seq-macro", +] + +[[package]] +name = "gemm-f64" +version = "0.19.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "056131e8f2a521bfab322f804ccd652520c79700d81209e9d9275bbdecaadc6a" +dependencies = [ + "dyn-stack", + "gemm-common", + "num-complex", + "num-traits", + "paste", + "raw-cpuid", + "seq-macro", +] + +[[package]] +name = "generativity" +version = "1.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d2c81fb5260e37854d09d5c87183309fd8c555b75289427884b25660bc87a85e" + +[[package]] +name = "getrandom" +version = "0.3.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "899def5c37c4fd7b2664648c28120ecec138e4d395b459e5ca34f9cce2dd77fd" +dependencies = [ + "cfg-if", + "libc", + "r-efi", + "wasip2", +] + +[[package]] +name = "half" +version = "2.7.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6ea2d84b969582b4b1864a92dc5d27cd2b77b622a8d79306834f1be5ba20d84b" +dependencies = [ + "bytemuck", + "cfg-if", + "crunchy", + "num-traits", + "zerocopy", +] + +[[package]] +name = "heck" +version = "0.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2304e00983f87ffb38b55b444b5e3b60a884b5d30c0fca7d82fe33449bbe55ea" + +[[package]] +name = "hermit-abi" +version = "0.5.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e17592d60ebacc7d5e169f4663c5f84f9161cc90328abcfe8456f41e4dfcb284" + +[[package]] +name = "interpol" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "eb58032ba748f4010d15912a1855a8a0b1ba9eaad3395b0c171c09b3b356ae50" +dependencies = [ + "proc-macro2", + "quote", + "syn 1.0.109", +] + +[[package]] +name = "libc" +version = "0.2.191" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "65cd87c9a6d60536580b316ad985726104730a5151c99845f2490428ef91240a" + +[[package]] +name = "libm" +version = "0.2.16" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b6d2cec3eae94f9f509c767b45932f1ada8350c4bdb85af2fcab4a3c14807981" + +[[package]] +name = "nano-gemm" +version = "0.2.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9e04345dc84b498ff89fe0d38543d1f170da9e43a2c2bcee73a0f9069f72d081" +dependencies = [ + "equator 0.2.2", + "nano-gemm-c32", + "nano-gemm-c64", + "nano-gemm-codegen", + "nano-gemm-core", + "nano-gemm-f32", + "nano-gemm-f64", + "num-complex", +] + +[[package]] +name = "nano-gemm-c32" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0775b1e2520e64deee8fc78b7732e3091fb7585017c0b0f9f4b451757bbbc562" +dependencies = [ + "nano-gemm-codegen", + "nano-gemm-core", + "num-complex", +] + +[[package]] +name = "nano-gemm-c64" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9af49a20d58816e6b5ee65f64142e50edb5eba152678d4bb7377fcbf63f8437a" +dependencies = [ + "nano-gemm-codegen", + "nano-gemm-core", + "num-complex", +] + +[[package]] +name = "nano-gemm-codegen" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6cc8d495c791627779477a2cf5df60049f5b165342610eb0d76bee5ff5c5d74c" + +[[package]] +name = "nano-gemm-core" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d998dfa644de87a0f8660e5ea511d7cb5c33b5a2d9847b7af57a2565105089f0" + +[[package]] +name = "nano-gemm-f32" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "879d962e79bc8952e4ad21ca4845a21132540ed3f5e01184b2ff7f720e666523" +dependencies = [ + "nano-gemm-codegen", + "nano-gemm-core", +] + +[[package]] +name = "nano-gemm-f64" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b9a513473dce7dc00c7e7c318481ca4494034e76997218d8dad51bd9f007a815" +dependencies = [ + "nano-gemm-codegen", + "nano-gemm-core", +] + +[[package]] +name = "num-complex" +version = "0.4.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "73f88a1307638156682bada9d7604135552957b7818057dcef22705b4d509495" +dependencies = [ + "bytemuck", + "num-traits", +] + +[[package]] +name = "num-traits" +version = "0.2.19" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "071dfc062690e90b734c0b2273ce72ad0ffa95f0c74596bc250dcfd960262841" +dependencies = [ + "autocfg", + "libm", +] + +[[package]] +name = "num_cpus" +version = "1.17.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "91df4bbde75afed763b708b7eee1e8e7651e02d97f6d5dd763e89367e957b23b" +dependencies = [ + "hermit-abi", + "libc", +] + +[[package]] +name = "once_cell" +version = "1.21.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9f7c3e4beb33f85d45ae3e3a1792185706c8e16d043238c593331cc7cd313b50" + +[[package]] +name = "paste" +version = "1.0.15" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "57c0d7b74b563b49d38dae00a0c37d4d6de9b432382b2892f0574ddcae73fd0a" + +[[package]] +name = "ppv-lite86" +version = "0.2.21" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "85eae3c4ed2f50dcfe72643da4befc30deadb458a9b590d720cde2f2b1e97da9" +dependencies = [ + "zerocopy", +] + +[[package]] +name = "private-gemm-x86" +version = "0.1.22" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bba0a7831d1297712e7e7b33759b0c837941e79196dd05df506d5f36f7f78619" +dependencies = [ + "defer", + "interpol", + "num_cpus", + "raw-cpuid", +] + +[[package]] +name = "proc-macro2" +version = "1.0.107" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "985e7ec9bb745e6ce6535b544d84d6cd6f7ad8bd711c398938ae983b91a766d9" +dependencies = [ + "unicode-ident", +] + +[[package]] +name = "pulp" +version = "0.22.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "046aa45b989642ec2e4717c8e72d677b13edd831a4d3b6cf37d9a3e54912496a" +dependencies = [ + "bytemuck", + "cfg-if", + "libm", + "num-complex", + "paste", + "pulp-wasm-simd-flag", + "raw-cpuid", + "reborrow", + "version_check", +] + +[[package]] +name = "pulp-wasm-simd-flag" +version = "0.1.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1d8f70e07b9c3962945a74e59ca1c511bba65b6419468acc217c457d93f3c740" + +[[package]] +name = "qd" +version = "0.8.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "15f1304a5aecdcfe9ee72fbba90aa37b3aa067a69d14cb7f3d9deada0be7c07c" +dependencies = [ + "bytemuck", + "libm", + "num-traits", + "pulp", +] + +[[package]] +name = "quote" +version = "1.0.47" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1fbf4db142a473a8d80c26bbf18454ed458bf8d26c8219c331daecfdbd079001" +dependencies = [ + "proc-macro2", +] + +[[package]] +name = "r-efi" +version = "5.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "69cdb34c158ceb288df11e18b4bd39de994f6657d83847bdffdbd7f346754b0f" + +[[package]] +name = "rand" +version = "0.9.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b9ef1d0d795eb7d84685bca4f72f3649f064e6641543d3a8c415898726a57b41" +dependencies = [ + "rand_chacha", + "rand_core", +] + +[[package]] +name = "rand_chacha" +version = "0.9.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d3022b5f1df60f26e1ffddd6c66e8aa15de382ae63b3a0c1bfc0e4d3e3f325cb" +dependencies = [ + "ppv-lite86", + "rand_core", +] + +[[package]] +name = "rand_core" +version = "0.9.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "76afc826de14238e6e8c374ddcc1fa19e374fd8dd986b0d2af0d02377261d83c" +dependencies = [ + "getrandom", +] + +[[package]] +name = "rand_distr" +version = "0.5.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6a8615d50dcf34fa31f7ab52692afec947c4dd0ab803cc87cb3b0b4570ff7463" +dependencies = [ + "num-traits", + "rand", +] + +[[package]] +name = "raw-cpuid" +version = "11.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "498cd0dc59d73224351ee52a95fee0f1a617a2eae0e7d9d720cc622c73a54186" +dependencies = [ + "bitflags", +] + +[[package]] +name = "reborrow" +version = "0.5.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "03251193000f4bd3b042892be858ee50e8b3719f2b08e5833ac4353724632430" + +[[package]] +name = "same-file" +version = "1.0.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "93fc1dc3aaa9bfed95e02e6eadabb4baf7e3078b0bd1b4d7b6b0b68378900502" +dependencies = [ + "winapi-util", +] + +[[package]] +name = "seq-macro" +version = "0.3.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1bc711410fbe7399f390ca1c3b60ad0f53f80e95c5eb935e52268a0e2cd49acc" + +[[package]] +name = "syn" +version = "1.0.109" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "72b64191b275b66ffe2469e8af2c1cfe3bafa67b529ead792a6d0160888b4237" +dependencies = [ + "proc-macro2", + "quote", + "unicode-ident", +] + +[[package]] +name = "syn" +version = "2.0.119" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "872831b642d1a07999a962a351ed35b955ea2cfc8f3862091e2a240a84f17297" +dependencies = [ + "proc-macro2", + "quote", + "unicode-ident", +] + +[[package]] +name = "syn" +version = "3.0.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d78c8dee4c7bf0e14673097256fed6142ce9d3b85a408189d07482442145823b" +dependencies = [ + "proc-macro2", + "quote", + "unicode-ident", +] + +[[package]] +name = "sysctl" +version = "0.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "01198a2debb237c62b6826ec7081082d951f46dbb64b0e8c7649a452230d1dfc" +dependencies = [ + "bitflags", + "byteorder", + "enum-as-inner", + "libc", + "thiserror", + "walkdir", +] + +[[package]] +name = "thiserror" +version = "1.0.69" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b6aaf5339b578ea85b50e080feb250a3e8ae8cfcdff9a461c9ec2904bc923f52" +dependencies = [ + "thiserror-impl", +] + +[[package]] +name = "thiserror-impl" +version = "1.0.69" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4fee6c4efc90059e10f81e6d42c60a18f76588c3d74cb83a0b242a2b6c7504c1" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.119", +] + +[[package]] +name = "unicode-ident" +version = "1.0.27" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a2c754d6c33795a1c324727428e5a7dedb5b06195f9890bdbcba760d3e246563" + +[[package]] +name = "version_check" +version = "0.9.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a" + +[[package]] +name = "walkdir" +version = "2.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "29790946404f91d9c5d06f9874efddea1dc06c5efe94541a7d6863108e3a5e4b" +dependencies = [ + "same-file", + "winapi-util", +] + +[[package]] +name = "wasip2" +version = "1.0.4+wasi-0.2.12" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b67efb37e106e55ce722a510d6b5f9c17f083e5fc79afc2badeb12cc313d9487" +dependencies = [ + "wit-bindgen", +] + +[[package]] +name = "winapi-util" +version = "0.1.11" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c2a7b1c03c876122aa43f3020e6c3c3ee5c05081c9a00739faf7503aeba10d22" +dependencies = [ + "windows-sys", +] + +[[package]] +name = "windows-link" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f0805222e57f7521d6a62e36fa9163bc891acd422f971defe97d64e70d0a4fe5" + +[[package]] +name = "windows-sys" +version = "0.61.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ae137229bcbd6cdf0f7b80a31df61766145077ddf49416a728b02cb3921ff3fc" +dependencies = [ + "windows-link", +] + +[[package]] +name = "wit-bindgen" +version = "0.57.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1ebf944e87a7c253233ad6766e082e3cd714b5d03812acc24c318f549614536e" + +[[package]] +name = "zerocopy" +version = "0.8.63" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e5fe1f8f1b06191a00962174c61aa5005e0bb391a6d80d07e24d115c01a92ed8" +dependencies = [ + "zerocopy-derive", +] + +[[package]] +name = "zerocopy-derive" +version = "0.8.63" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "863ad3ac83293fb4d740aedbfdc9240dd8d1a50c1099acd76ce80ce7c7230c7f" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.119", +] diff --git a/upstream-reports/faer-029-real-schur-aed-sort-hang/Cargo.toml b/upstream-reports/faer-029-real-schur-aed-sort-hang/Cargo.toml new file mode 100644 index 00000000..e87c2d04 --- /dev/null +++ b/upstream-reports/faer-029-real-schur-aed-sort-hang/Cargo.toml @@ -0,0 +1,13 @@ +[package] +name = "faer-029-real-schur-aed-sort-hang" +version = "0.0.0" +edition = "2021" +publish = false + +# Not part of the t4a-faer workspace: this reproduces the published crate on its own. +[workspace] + +[dependencies] +faer = { version = "=0.24.4", default-features = false, features = ["std", "linalg"] } +rand = "0.9" +rand_distr = "0.5" diff --git a/upstream-reports/faer-029-real-schur-aed-sort-hang/README.md b/upstream-reports/faer-029-real-schur-aed-sort-hang/README.md new file mode 100644 index 00000000..e1e83b64 --- /dev/null +++ b/upstream-reports/faer-029-real-schur-aed-sort-hang/README.md @@ -0,0 +1,52 @@ +# faer-029: real Schur QR does not return when the deflation-window sort cannot swap a pair + +**What happens:** `Mat::eigen` never returns for some finite `f32` matrices; the thread spins in the real Schur QR. +**Which input:** the near-Jordan matrix `J_n + 1e-6 G` (Jordan block with eigenvalue 1, `G` standard normal from `StdRng` seed 112), divided by its largest entry and scaled by `2^-100`, at n = 128 and n = 257. +**Which version:** `faer` 0.24.4. +**Added to this list:** 2026-10-29. + +## Reproduce + +```sh +cargo run --release +``` + +The program runs `Mat::eigen` on that matrix for n = 128 and n = 257 in a separate thread and waits 60 s for each. + +| n | expected | `faer` 0.24.4 | +|---|---|---| +| 128 | returns | **no return after 60 s** | +| 257 | returns | **no return after 60 s** | + +The same matrices at scales `2^0` to `2^-90` return. +A backtrace taken while it spins ends in `schur_swap_22` called from `aggressive_early_deflation`. + +## Cause + +`aggressive_early_deflation` sorts the eigenvalues of the deflation window with a bubble sort (`src/linalg/evd/schur/real_schur.rs:921-969`). +It is a port of LAPACK `dlaqr3`, with two differences: + +* `dlaqr3` sorts the undeflated part `T(INFQR+1:NS)`. The port sorts from `ns` to `jw`, the deflated part: it took `dlaqr3`'s initial `I = NS + 1`, which only sets the first pass's end `KEND = I - 1`, as the start of the sort. +* `dlaqr3` ends each pass at `KEND`, one before where the previous pass stopped, so the window shrinks by at least one per pass and the sort ends whatever the swaps do. The port sets the next window from the position of the last swap. + +`schur_swap_22` rejects an ill-conditioned swap of two 2 x 2 blocks and leaves the matrix unchanged, but `schur_swap` reports success either way. +The pass then records the same swap position as the next window, the pair is still out of order, and the next pass makes the same attempt, forever. + +## How the fork fixed it + +The sort follows `dlaqr3`: it sorts the undeflated part and bounds each pass by `KEND`. +It swaps with the checked swap, which uses the `dlaexc` ports for the 1 x 2 and 2 x 1 cases, and a rejected swap leaves the pair in place, as `dlaqr3` does after `DTREXC` fails. +With that change the fork returns for both sizes within a few seconds. +At this scale the input lies below LAPACK's negligibility level `smlnum = safmin * n / eps`, so the Schur form reproduces it only to within `smlnum` in absolute terms, as LAPACK's `?HSEQR` does for unscaled input; the same matrix scaled to `2^-80` gives the usual small relative residual. + +Fork regression test, in `faer/src/linalg/evd/schur/real_schur.rs`: `tests::aed_sort_terminates_when_a_swap_fails`. + +## Affected versions + +Confirmed by running this program: `faer` 0.24.4. +Not fixed upstream: the same code is in upstream `main` at `0ecd846` (2026-10-22). + +## References + +* Fork fix: tensor4all/faer#144 (PR to follow). +* Where it was found: tensor4all/faer#112, the LAPACK comparison sweep. diff --git a/upstream-reports/faer-029-real-schur-aed-sort-hang/src/main.rs b/upstream-reports/faer-029-real-schur-aed-sort-hang/src/main.rs new file mode 100644 index 00000000..04714d68 --- /dev/null +++ b/upstream-reports/faer-029-real-schur-aed-sort-hang/src/main.rs @@ -0,0 +1,45 @@ +// `Mat::eigen` of an f32 near-Jordan matrix, J_n + 1e-6 G (Jordan block with eigenvalue 1, +// G standard normal from StdRng seed 112), divided by its largest entry and scaled by 2^-100. +// All entries are finite, normal f32 values. +// The call should return in well under a second; the program gives each size 60 s. +use faer::Mat; +use rand::prelude::*; +use rand_distr::StandardNormal; +use std::sync::mpsc; +use std::time::{Duration, Instant}; + +fn input(n: usize) -> Mat { + let rng = &mut StdRng::seed_from_u64(112); + let a = Mat::::from_fn(n, n, |i, j| { + let g: f32 = rng.sample(StandardNormal); + let jordan = if i == j || i + 1 == j { 1.0 } else { 0.0 }; + jordan + 1e-6 * g + }); + let max = a.norm_max(); + Mat::from_fn(n, n, |i, j| a[(i, j)] / max * 2f32.powi(-100)) +} + +fn main() { + let mut bad = false; + for n in [128, 257] { + let a = input(n); + let (send, recv) = mpsc::channel(); + let start = Instant::now(); + std::thread::spawn(move || { + let ok = a.eigen().is_ok(); + let _ = send.send(ok); + }); + match recv.recv_timeout(Duration::from_secs(60)) { + Ok(ok) => println!("n = {n}: returned after {:.2} s (Ok: {ok})", start.elapsed().as_secs_f64()), + Err(_) => { + println!("n = {n}: no return after 60 s"); + bad = true; + } + } + } + println!("expected: both calls return"); + if bad { + println!("WRONG: Mat::eigen does not return"); + std::process::exit(1); + } +}