par: no heap allocation inside a parallel region

PATCHSET.md

@@ -130,6 +130,7 @@ gate has changed and two kinds of adoption are now in the ledger.
 | Threadteam is the only parallel route | `Par::Rayon`, `Par::rayon`, `ParShape::Rayon`, `ParShape::rayon` and `Par::__Lifetime` are removed, with the Rayon arms of `thread::join_raw`/`for_each`/`for_each_with_par`/`split`, of `ParShape::dispatch_cost`/`min_dispatch_cost`, of `par_pays_time`/`par_can_pay_time` and of the GEMM route; `parallelism_degree` folds into `Par::degree` and `is_team` goes. The Rayon cost model in `utils::cost` (`outside_cost_of`, `in_pool_cost_of`, `width_model`, `rayon_cost`, `min_rayon_cost`, `caller_in_pool`, the outside/in-pool anchors) and its tests are deleted; `ParShape::dispatch_cost` reads the costs of the team token the shape came from. faer's own serial `Rates` table and its `HostClass` key stay, and the `PAYS`/`pays_at` rule has one copy in `utils::cost` with a test against `threadteam`'s. `spindle` leaves `Cargo.toml` (no production code opened a group; the last users were tests and the two example harnesses). New: `faer::team::with_team(width, |par| ..)` owns a `threadteam::Pool` for a closure and hands out `Par::Team` (a width of one hands out `Par::Seq`), which is the caller-side answer #59 was waiting for and what the benchmark harness and the threshold probe use instead of `Par::rayon(n)`. `threadteam` becomes the feature that compiles the route (on by default, pulling `rayon` as its substrate); with it off there is no `Par::Team` and no team module, so a no_std build has `Par::Seq` only and nothing pulls std in through Rayon. 81 lib-test and 15 example sites converted to `with_team`; the tests that only asserted the removed route are deleted with the route. | KEEP | tensor4all/faer#85, tensor4all/faer#59 | Mechanism, not a speed-up: the acceptance measurement pairs the branch against `origin/main` with both rows on the team route, and the claim is only that deleting the second route did not regress the team route. Verification: `cargo test -p t4a-faer --lib`, `cargo test -p t4a-faer --doc`, `cargo check` for `--no-default-features`, `--no-default-features --features std`, `--no-default-features --features std,rayon`, default and `std,rayon,threadteam`, and `cargo +nightly-2025-01-08 run --profile nostd` in `faer-no-std-test` (the inherited nostd job). Record: [worklog](docs/worklogs/2026-10-29-threadteam-only.md). |
 | 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). |
 
 ## Ground rules
 

docs/worklogs/2026-10-29-no-region-allocation.md

@@ -0,0 +1,52 @@
+# No heap allocation inside a parallel region (tensor4all/faer#111)
+
+## Decisions
+- **What allocated, per region, on `origin/main` `ad23987`:**
+  - `utils::team::for_each_child`, behind every `for_each` and `for_each_with_par`, collected every item of the
+    iterator into a `Vec` (one entry per column for `par_col_iter_mut`), then built a second `Vec` of slots.
+  - The full-pivoting LU step built a `best` vector of `nthreads` pivot candidates, and so did both
+    Bunch-Kaufman full-pivoting updates (rank one and rank two).
+  - `utils::thread::join_sized` boxed its two closures, which runs for each LU swap stage.
+- **Why.** The monitor ruled that these allocations are a mechanism defect, not a host constant.
+  - They were also suspected of causing #109's floor-level loss: `full_piv_lu f64 128` at 8T on evo was +6.6 to
+    +10.4% after #109 ran its per-step regions down to the rule's floor.
+  - **That suspicion was wrong.** The evo re-measurement of this change (2026-10-29 07:54Z) puts the cell at
+    1.0025 of `main` with the allocations gone, and finds no regression elsewhere.
+  - The change stands on its own as a cleanup of the mechanism.
+  - The remaining candidates for the floor-level loss are the fixed per-section costs of the 19 steps that
+    narrow: the reduction merge, and the trailing block's columns moving between cores as the width changes.
+    Whether to chase them is decided later.
+- **`for_each_child`.**
+  - It takes the producer through `with_producer`, splits it on the leader into the blocks with `split_at`,
+    and stores the pieces. The splits are relative to what remains; the block boundaries are those of
+    `par_split_indices`.
+  - Each block's leader takes its piece and iterates it there, so `P::IntoIter` never has to be `Send`.
+  - Everything runs inside the producer callback, so producers that borrow from its frame stay valid:
+    threadteam joins every share before it returns.
+  - **Storage.** An inline array of 64 `Option<P>` slots is used when there are at most 64 blocks and a slot
+    is at most 256 bytes, so the array is at most 16 KiB. A producer can carry state of any size; rayon's
+    `map_with` stores its state by value. Otherwise the producer is halved recursively with `team.split`
+    and two slots per level, which costs one section per level and allocates nothing either.
+- **Pivot candidates.** They come from the call's scratch:
+  - `full_pivoting::factor::lu_in_place_scratch` and `lblt::factor::cholesky_in_place_scratch` grow by
+    `par.degree()` entries of `(usize, usize, Real)`;
+  - the buffer is passed down to the step loop and reused by every step.
+
+  In `cholesky_in_place_scratch` it sits in an `or` with the blocked variant's workspace, since a call uses
+  one or the other.
+- **`join_sized`.** Two typed slots and a match on the child index, instead of `Box<dyn FnOnce>`.
+- **Not changed.**
+  - The householder upgrade's `partials` vector. It is allocated per call, not per step.
+  - A fixed-cost term in the cost model. Whether one is still needed is decided after the re-measurement.
+
+## Verification conclusions and constraints
+- **Allocation test:** `faer/tests/no_step_alloc.rs` is its own test binary, so it can own a `#[global_allocator]`.
+  It counts every allocation of the process, on any thread, while one call runs. Setup: teams of 4, 8 and 70
+  (70 takes the recursive path), regions forced parallel with `ParThreshold::Fixed(0)`, a warm-up call first.
+  - It finds 0 allocations for full-pivoting LU 96, Bunch-Kaufman with full pivoting 96, and `join_sized`.
+  - It checks afterwards, outside the counted interval, that each call dispatched a parallel section.
+  - On `origin/main` the same test counts 190 allocations in one full-pivoting LU 96 call at width 4.
+- **Scope of the promise.** It covers faer's own dispatch and workspaces. An iterator adapter that allocates in
+  its own `split_at` (rayon's `map_with` clones its state) still does.
+- **Performance**, measured by the monitor on evo: no regression. `full_piv_lu f64 128` at 8T is at 1.0025 of
+  `main`, so #109's +7 to 10% on that cell remains.

faer/src/linalg/cholesky/bunch_kaufman/factor.rs

@@ -140,8 +140,11 @@ fn rank_1_update_and_argmax<T: ComplexField>(
 	A: MatMut<'_, T>,
 	L: ColRef<'_, T>,
 	d: T::Real,
+	best: &mut [(usize, usize, T::Real)],
 	par: Par,
 ) -> (usize, usize, T::Real) {
+	// the parallel arm is the only reader (absent without `rayon`)
+	let _ = &best;
 	with_dim!(N, A.nrows());
 	match par {
 		Par::Seq => rank_1_update_and_argmax_seq(
@@ -164,8 +167,9 @@ fn rank_1_update_and_argmax<T: ComplexField>(
 					1.0f64 - libm::sqrt(1.0f64 - idx_as_percent);
 				(col_start_percent * n as f64) as usize
 			};
-			let mut r =
-				alloc::vec![(0usize, 0usize, zero::< T::Real > ()); nthreads];
+			// one entry per block, from the call's scratch (no allocation per step, #111); every
+			// entry is written by its block
+			let r = &mut best[..nthreads];
 			crate::utils::thread::for_each(par, 
 				nthreads,
 				r.par_iter_mut().enumerate(),
@@ -182,7 +186,8 @@ fn rank_1_update_and_argmax<T: ComplexField>(
 					);
 				},
 			);
-			r.into_iter()
+			r.iter()
+				.map(|(i, j, v)| (*i, *j, v.copy()))
 				.max_by(|(_, _, a), (_, _, b)| {
 					if a == b {
 						core::cmp::Ordering::Equal
@@ -204,8 +209,11 @@ fn rank_2_update_and_argmax<'N, T: ComplexField>(
 	d00: T::Real,
 	d11: T::Real,
 	d10: T,
+	best: &mut [(usize, usize, T::Real)],
 	par: Par,
 ) -> (usize, usize, T::Real) {
+	// the parallel arm is the only reader (absent without `rayon`)
+	let _ = &best;
 	with_dim!(N, A.nrows());
 	match par {
 		Par::Seq => rank_2_update_and_argmax_seq(
@@ -232,8 +240,9 @@ fn rank_2_update_and_argmax<'N, T: ComplexField>(
 					1.0f64 - libm::sqrt(1.0f64 - idx_as_percent);
 				(col_start_percent * n as f64) as usize
 			};
-			let mut r =
-				alloc::vec![(0usize, 0usize, zero::< T::Real > ()); nthreads];
+			// one entry per block, from the call's scratch (no allocation per step, #111); every
+			// entry is written by its block
+			let r = &mut best[..nthreads];
 			crate::utils::thread::for_each(par, 
 				nthreads,
 				r.par_iter_mut().enumerate(),
@@ -254,7 +263,8 @@ fn rank_2_update_and_argmax<'N, T: ComplexField>(
 					);
 				},
 			);
-			r.into_iter()
+			r.iter()
+				.map(|(i, j, v)| (*i, *j, v.copy()))
 				.max_by(|(_, _, a), (_, _, b)| {
 					if a == b {
 						core::cmp::Ordering::Equal
@@ -272,6 +282,7 @@ fn lblt_full_piv<T: ComplexField>(
 	A: MatMut<'_, T>,
 	subdiag: DiagMut<'_, T>,
 	pivots: &mut [usize],
+	best: &mut [(usize, usize, T::Real)],
 	par: Par,
 	params: LbltParams,
 ) {
@@ -359,7 +370,7 @@ fn lblt_full_piv<T: ComplexField>(
 			}
 			if n != 0 {
 				(max_i, max_j, max_offdiag) = with_width(par.reborrow(), width, |par| {
-					rank_1_update_and_argmax(A.rb_mut(), L.rb(), diag, par)
+					rank_1_update_and_argmax(A.rb_mut(), L.rb(), diag, best, par)
 				});
 			}
 		} else {
@@ -389,6 +400,7 @@ fn lblt_full_piv<T: ComplexField>(
 						d00.copy(),
 						d11.copy(),
 						d10.copy(),
+						best,
 						par,
 					)
 				});
@@ -1158,12 +1170,16 @@ pub fn cholesky_in_place_scratch<I: Index, T: ComplexField>(
 	params: Spec<LbltParams, T>,
 ) -> StackReq {
 	let params = params.config;
-	let _ = par;
 	let mut bs = params.block_size;
 	if bs < 2 || dim <= bs {
 		bs = 0;
 	}
-	StackReq::new::<usize>(dim).and(temp_mat_scratch::<T>(dim, bs))
+	// the pivots, then the blocked variant's workspace or the full-pivoting variant's pivot
+	// candidates, one per block of a parallel step
+	StackReq::new::<usize>(dim).and(
+		temp_mat_scratch::<T>(dim, bs)
+			.or(StackReq::new::<(usize, usize, T::Real)>(par.degree())),
+	)
 }
 /// info about the result of the $LBL^\top$ factorization
 #[derive(Copy, Clone, Debug)]
@@ -1234,7 +1250,8 @@ pub fn cholesky_in_place<'out, I: Index, T: ComplexField>(
 		_ => (false, false),
 	};
 	if params.pivoting == PivotingStrategy::Full {
-		lblt_full_piv(A, subdiag, pivots, par, params);
+		let (mut best, _) = stack.make_with(par.degree(), |_| (0usize, 0usize, zero::<T::Real>()));
+		lblt_full_piv(A, subdiag, pivots, best.as_mut(), par, params);
 	} else {
 		lblt_blocked(A, subdiag, pivots, bs, rook, diagonal, par, stack);
 	}

faer/src/linalg/lu/full_pivoting/factor.rs

@@ -342,10 +342,13 @@ fn lu_in_place_unblocked<T: ComplexField>(
 	A: MatMut<'_, T>,
 	row_trans: &mut [usize],
 	col_trans: &mut [usize],
+	best: &mut [(usize, usize, T::Real)],
 	par: Par,
 	transpose: bool,
 	params: Spec<FullPivLuParams, T>,
 ) -> usize {
+	// the parallel arm is the only reader (absent without `rayon`)
+	let _ = &best;
 	let params = params.config;
 	let mut n_trans = 0;
 	let (m, n) = A.shape();
@@ -420,9 +423,9 @@ fn lu_in_place_unblocked<T: ComplexField>(
 				let nthreads = core::num::NonZeroUsize::new(par.degree()).unwrap();
 					use rayon::prelude::*;
 					let nthreads = nthreads.get();
-					let mut best = core::iter::repeat_with(|| (0, 0, zero()))
-						.take(nthreads)
-						.collect::<alloc::vec::Vec<_>>();
+					// one entry per block, from the call's scratch (no allocation per step, #111); every
+					// entry is written by its block
+					let best = &mut best[..nthreads];
 					let full_cols = A11.ncols();
 					crate::utils::thread::for_each(par, 
 						nthreads,
@@ -447,11 +450,11 @@ fn lu_in_place_unblocked<T: ComplexField>(
 					max_row = 0;
 					max_col = 0;
 					max_score = zero();
-					for (row, col, val) in best {
-						if val > max_score {
-							max_row = row;
-							max_col = col;
-							max_score = val;
+					for (row, col, val) in best.iter() {
+						if *val > max_score {
+							max_row = *row;
+							max_col = *col;
+							max_score = val.copy();
 						}
 					}
 				},
@@ -488,10 +491,12 @@ pub fn lu_in_place_scratch<I: Index, T: ComplexField>(
 	par: ParShape,
 	params: Spec<FullPivLuParams, T>,
 ) -> StackReq {
-	_ = par;
 	_ = params;
 	let size = Ord::min(nrows, ncols);
-	StackReq::new::<usize>(size).array(2)
+	// the transpositions, and one pivot candidate per block of a parallel step
+	StackReq::new::<usize>(size)
+		.array(2)
+		.and(StackReq::new::<(usize, usize, T::Real)>(par.degree()))
 }
 #[derive(Copy, Clone, Debug)]
 pub struct FullPivLuInfo {
@@ -542,13 +547,16 @@ pub fn lu_in_place<'out, I: Index, T: ComplexField>(
 	let size = Ord::min(M, N);
 	let (mut row_transpositions, stack) = stack.make_with(size, |_| 0);
 	let row_transpositions = row_transpositions.as_mut();
-	let (mut col_transpositions, _) = stack.make_with(size, |_| 0);
+	let (mut col_transpositions, stack) = stack.make_with(size, |_| 0);
 	let col_transpositions = col_transpositions.as_mut();
+	let (mut best, _) = stack.make_with(par.degree(), |_| (0usize, 0usize, zero::<T::Real>()));
+	let best = best.as_mut();
 	let n_transpositions = if mat.row_stride().abs() < mat.col_stride().abs() {
 		lu_in_place_unblocked(
 			mat,
 			row_transpositions,
 			col_transpositions,
+			best,
 			par,
 			false,
 			params,
@@ -558,6 +566,7 @@ pub fn lu_in_place<'out, I: Index, T: ComplexField>(
 			mat.transpose_mut(),
 			col_transpositions,
 			row_transpositions,
+			best,
 			par,
 			true,
 			params,

faer/src/utils/mod.rs

@@ -200,14 +200,17 @@ pub mod thread {
 			return;
 		}
 		let size_a = size_a.clamp(1, degree - 1);
-		// one slot type for two different closure types
-		let ops: [Slot<alloc::boxed::Box<dyn FnOnce(Par) + Send>>; 2] = [
-			Slot::new(alloc::boxed::Box::new(op_a)),
-			Slot::new(alloc::boxed::Box::new(op_b)),
-		];
+		// a typed slot per closure, so nothing is boxed (#111: no allocation per region)
+		let (a, b) = (Slot::new(op_a), Slot::new(op_b));
 		split(par, &[size_a, degree - size_a], &|c, child| {
 			// SAFETY: the split hands index `c` to one child only, and `c < 2`
-			unsafe { (ops[c].take())(child) }
+			unsafe {
+				if c == 0 {
+					(a.take())(child)
+				} else {
+					(b.take())(child)
+				}
+			}
 		});
 	}