evd: the real Schur deflation-window sort always ends (#144)

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
 

faer/src/linalg/evd/schur/real_schur.rs

@@ -1210,56 +1210,52 @@ fn aggressive_early_deflation<T: RealField>(
 		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::<f32>::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::<f32>::zeros(1, n - 1);
+			hessenberg_in_place(
+				a.as_mut(),
+				householder.as_mut(),
+				Par::Seq,
+				MemStack::new(&mut MemBuffer::new(
+					hessenberg_in_place_scratch::<f32>(
+						n,
+						1,
+						ParShape::Seq,
+						Default::default(),
+					),
+				)),
+				Default::default(),
+			);
+			let h = Mat::<f32>::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::<f32>::identity(n, n);
+			let mut wr = Col::<f32>::zeros(n);
+			let mut wi = Col::<f32>::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::<f32>(
+						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::<f32>::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}"
+				);
+			}
+		}
+	}
 }

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).