diff --git a/PATCHSET.md b/PATCHSET.md index b60193b4..c7128220 100644 --- a/PATCHSET.md +++ b/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 diff --git a/docs/worklogs/2026-10-29-no-region-allocation.md b/docs/worklogs/2026-10-29-no-region-allocation.md new file mode 100644 index 00000000..86a2ec5c --- /dev/null +++ b/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

` 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`. +- **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. diff --git a/faer/src/linalg/cholesky/bunch_kaufman/factor.rs b/faer/src/linalg/cholesky/bunch_kaufman/factor.rs index dc2b18b5..daf12d2a 100644 --- a/faer/src/linalg/cholesky/bunch_kaufman/factor.rs +++ b/faer/src/linalg/cholesky/bunch_kaufman/factor.rs @@ -140,8 +140,11 @@ fn rank_1_update_and_argmax( 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( 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( ); }, ); - 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( 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( } 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( d00.copy(), d11.copy(), d10.copy(), + best, par, ) }); @@ -1158,12 +1170,16 @@ pub fn cholesky_in_place_scratch( params: Spec, ) -> StackReq { let params = params.config; - let _ = par; let mut bs = params.block_size; if bs < 2 || dim <= bs { bs = 0; } - StackReq::new::(dim).and(temp_mat_scratch::(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::(dim).and( + temp_mat_scratch::(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::())); + lblt_full_piv(A, subdiag, pivots, best.as_mut(), par, params); } else { lblt_blocked(A, subdiag, pivots, bs, rook, diagonal, par, stack); } diff --git a/faer/src/linalg/lu/full_pivoting/factor.rs b/faer/src/linalg/lu/full_pivoting/factor.rs index 5e5d1d02..b4aba20a 100644 --- a/faer/src/linalg/lu/full_pivoting/factor.rs +++ b/faer/src/linalg/lu/full_pivoting/factor.rs @@ -342,10 +342,13 @@ fn lu_in_place_unblocked( A: MatMut<'_, T>, row_trans: &mut [usize], col_trans: &mut [usize], + best: &mut [(usize, usize, T::Real)], par: Par, transpose: bool, params: Spec, ) -> 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( 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::>(); + // 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( 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( par: ParShape, params: Spec, ) -> StackReq { - _ = par; _ = params; let size = Ord::min(nrows, ncols); - StackReq::new::(size).array(2) + // the transpositions, and one pivot candidate per block of a parallel step + StackReq::new::(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::())); + 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, diff --git a/faer/src/utils/mod.rs b/faer/src/utils/mod.rs index 9ab7c8b6..e523785a 100644 --- a/faer/src/utils/mod.rs +++ b/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>; 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) + } + } }); } diff --git a/faer/src/utils/team.rs b/faer/src/utils/team.rs index 27751dc0..94b8adea 100644 --- a/faer/src/utils/team.rs +++ b/faer/src/utils/team.rs @@ -12,7 +12,6 @@ //! Nothing here starts a Rayon job: the team members are threadteam members that spin while the //! reservation is held, and the section dispatch is a write to their slots. -use alloc::vec::Vec; use core::cell::UnsafeCell; use rayon::iter::plumbing::{Producer, ProducerCallback}; use rayon::iter::{IndexedParallelIterator, IntoParallelIterator}; @@ -207,25 +206,6 @@ pub(crate) fn run(team: TeamExec<'_>, n: usize, f: &(dyn Fn(usize, usize) + Sync } } -struct Slots(Vec>>); -// SAFETY: every slot is taken by exactly one participant. -unsafe impl Sync for Slots {} -impl Slots { - /// # Safety - /// Slot `j` is taken by one participant only. - unsafe fn take(&self, j: usize) -> T { - unsafe { (*self.0[j].get()).take() }.unwrap() - } -} - -struct Collect<'a, T>(&'a mut Vec); -impl ProducerCallback for Collect<'_, T> { - type Output = (); - fn callback>(self, producer: P) { - self.0.extend(producer.into_iter()); - } -} - /// A flat item loop on the team: the items of `iter` are dealt out in contiguous static blocks /// over `min(n_jobs, items, team size)` parts, each part processed by the leader of a sub-team. pub(crate) fn for_each( @@ -241,6 +221,11 @@ pub(crate) fn for_each( } /// [`for_each`] whose closure also receives the sub-team of its block, lent for the call. +/// +/// The items are not collected: the producer of `iter` is split on the leader into the blocks +/// (`split_at`), and each block is iterated by the leader of its part, so the dispatch allocates +/// nothing whatever the team's width (#111). An iterator adapter can still allocate in its own +/// `split_at` (rayon's `map_with` clones its state). pub(crate) fn for_each_child( team: TeamExec<'_>, n_jobs: usize, @@ -250,24 +235,102 @@ pub(crate) fn for_each_child( I: IntoParallelIterator, I::Iter: IndexedParallelIterator, { - let mut items = Vec::new(); - // `with_producer` walks the splitting machinery without starting any Rayon job. - iter.into_par_iter().with_producer(Collect(&mut items)); - let len = items.len(); + let iter = iter.into_par_iter(); + let len = iter.len(); if len == 0 { return; } - let slots = Slots(items.into_iter().map(|x| UnsafeCell::new(Some(x))).collect()); let jobs = n_jobs.max(1).min(len).min(team.size().get()); stats::hist(&stats::SPLIT, jobs); - team.fork(jobs, |c, mut child| { - let (start, n) = super::thread::par_split_indices(len, c, jobs); - for j in start..start + n { - // SAFETY: block `c` owns slots `start..start + n` alone. - let item = unsafe { slots.take(j) }; - f(item, child.reborrow()); + struct Deal<'t, 'f, F> { + team: TeamExec<'t>, + len: usize, + jobs: usize, + f: &'f F, + } + impl Fn(T, TeamExec<'c>)> ProducerCallback for Deal<'_, '_, F> { + type Output = (); + fn callback>(self, producer: P) { + // the producer may borrow from this callback's frame: every block is run and joined + // before it returns + deal(self.team, producer, self.len, self.jobs, self.f) } - }); + } + iter.with_producer(Deal { team, len, jobs, f: &f }); +} + +/// Blocks a single fork deals out of an inline array: up to this many, and only while the array +/// stays small (`FLAT_BYTES` per slot, 16 KiB in all), since a producer can carry state of any +/// size. Wider or heavier splits halve recursively instead, which costs a section per level +/// and allocates nothing either. +const FLAT_JOBS: usize = 64; +const FLAT_BYTES: usize = 256; + +/// One producer piece, taken by one participant. +struct Piece

(UnsafeCell>); +// SAFETY: each piece is taken by exactly one participant (the one-taker invariant below). +unsafe impl Sync for Piece

{} +impl

Piece

{ + /// # Safety + /// Called once per piece, by the participant that owns it. + unsafe fn take(&self) -> P { + unsafe { (*self.0.get()).take() }.expect("each piece is taken once") + } +} + +/// Runs the `len` items of `producer` as `jobs` contiguous blocks (the boundaries of +/// `par_split_indices`), block `c` on the leader of the `c`-th part of `team`. +fn deal, F: Sync + for<'c> Fn(T, TeamExec<'c>)>( + mut team: TeamExec<'_>, + producer: P, + len: usize, + jobs: usize, + f: &F, +) { + if jobs <= 1 { + for item in producer.into_iter() { + f(item, team.reborrow()); + } + return; + } + if jobs <= FLAT_JOBS && core::mem::size_of::>() <= FLAT_BYTES { + let mut pieces: [Piece

; FLAT_JOBS] = core::array::from_fn(|_| Piece(UnsafeCell::new(None))); + let mut rest = producer; + for (c, piece) in pieces.iter_mut().enumerate().take(jobs - 1) { + // split relative to what is left: the earlier blocks are already taken off + let (_, n) = super::thread::par_split_indices(len, c, jobs); + let (head, tail) = rest.split_at(n); + *piece.0.get_mut() = Some(head); + rest = tail; + } + *pieces[jobs - 1].0.get_mut() = Some(rest); + team.fork(jobs, |c, mut child| { + // SAFETY: the fork hands index `c` to one child only + let piece = unsafe { pieces[c].take() }; + for item in piece.into_iter() { + f(item, child.reborrow()); + } + }); + return; + } + // two halves, each with the participants of its blocks + let left = jobs / 2; + let (left_len, _) = super::thread::par_split_indices(len, left, jobs); + let size = team.size().get(); + let (left_size, _) = super::thread::par_split_indices(size, left, jobs); + let (a, b) = producer.split_at(left_len); + let halves = [Piece(UnsafeCell::new(Some(a))), Piece(UnsafeCell::new(Some(b)))]; + note_split(2); + team.split(&[left_size, size - left_size], |c, child| { + // SAFETY: the split hands index `c` to one child only, and `c < 2` + let half = unsafe { halves[c].take() }; + if c == 0 { + deal(child, half, left_len, left, f) + } else { + deal(child, half, len - left_len, jobs - left, f) + } + }) + .expect("two non-empty parts within the team"); } struct SendPtr(*mut T); @@ -1064,6 +1127,90 @@ mod tests { assert_eq!(parallel, 0, "{s:?}"); } + /// The producer dealing of `for_each_child` (#111): every item exactly once, each on a child + /// of an even share of the team, on the flat path, on the recursive one (more blocks than + /// `FLAT_JOBS`, or a producer over `FLAT_BYTES`: `map_with` keeps its 512-byte state inline), + /// with item, block and team counts that do not divide each other. + #[test] + fn for_each_child_deals_every_item_once() { + use rayon::prelude::*; + for (width, len, n_jobs, big) in [ + (8usize, 1000usize, 7usize, false), + (8, 1000, 7, true), + (70, 1000, 70, false), + (70, 333, 69, true), + (5, 3, 5, false), + ] { + with_par(width, |par| { + let Par::Team(team) = par else { unreachable!() }; + let seen: Vec = (0..len).map(|_| AtomicUsize::new(0)).collect(); + let widths: Vec = (0..len).map(|_| AtomicUsize::new(0)).collect(); + let f = |i: usize, child: TeamExec<'_>| { + seen[i].fetch_add(1, Ordering::Relaxed); + widths[i].store(child.size().get(), Ordering::Relaxed); + }; + if big { + core::assert!(core::mem::size_of::<[u8; 512]>() > FLAT_BYTES); + for_each_child(team, n_jobs, (0..len).into_par_iter().map_with([0u8; 512], |_, i| i), f); + } else { + for_each_child(team, n_jobs, (0..len).into_par_iter(), f); + } + let jobs = n_jobs.min(len).min(width); + let (lo, hi) = (width / jobs, width.div_ceil(jobs)); + for i in 0..len { + core::assert_eq!(seen[i].load(Ordering::Relaxed), 1, "{width} {len} {n_jobs} {big}: item {i}"); + let w = widths[i].load(Ordering::Relaxed); + core::assert!(lo <= w && w <= hi, "{width} {len} {n_jobs} {big}: item {i} on a child of {w}"); + } + }); + } + } + + /// A panic in one block reaches the caller after every block has finished, the items are + /// dropped once each, and the team stays usable. + #[test] + fn for_each_child_unwinds_and_keeps_the_team() { + use rayon::prelude::*; + use std::panic::{catch_unwind, AssertUnwindSafe}; + struct Counted<'a>(&'a AtomicUsize); + impl Drop for Counted<'_> { + fn drop(&mut self) { + self.0.fetch_add(1, Ordering::Relaxed); + } + } + for width in [8usize, 70] { + with_par(width, |mut par| { + let dropped = AtomicUsize::new(0); + let made = AtomicUsize::new(0); + let r = catch_unwind(AssertUnwindSafe(|| { + let Par::Team(team) = par.reborrow() else { unreachable!() }; + for_each_child( + team, + width, + (0..1000usize).into_par_iter().map(|i| { + made.fetch_add(1, Ordering::Relaxed); + (i, Counted(&dropped)) + }), + |(i, _c), _| { + if i == 500 { + panic!("block failure"); + } + }, + ); + })); + core::assert!(r.is_err()); + core::assert_eq!(dropped.load(Ordering::Relaxed), made.load(Ordering::Relaxed)); + // the team still works + let Par::Team(team) = par.reborrow() else { unreachable!() }; + let hits = AtomicUsize::new(0); + for_each_child(team, width, (0..100usize).into_par_iter(), |_, _| { + hits.fetch_add(1, Ordering::Relaxed); + }); + core::assert_eq!(hits.load(Ordering::Relaxed), 100); + }); + } + } + /// `par_pays_time` is `ParShape::pays` without the shape. #[test] fn par_pays_time_is_the_shape_gate() { diff --git a/faer/tests/no_step_alloc.rs b/faer/tests/no_step_alloc.rs new file mode 100644 index 00000000..2e3c0332 --- /dev/null +++ b/faer/tests/no_step_alloc.rs @@ -0,0 +1,124 @@ +//! #111: a parallel region allocates nothing on the heap. Each counted call below runs on a held +//! team with its regions forced parallel (`ParThreshold::Fixed(0)`), after a warm-up, and the +//! allocator counts every allocation of the process (any thread) while the call runs. The team, +//! the matrices and the scratch are built outside the counted interval, and the check that the +//! regions did dispatch reads the team counters after it. +#![cfg(all(feature = "std", feature = "threadteam"))] +use faer::dyn_stack::{MemBuffer, MemStack}; +use faer::utils::cost::ParThreshold; +use faer::utils::team::stats; +use faer::{Auto, Mat, Par}; +use std::alloc::{GlobalAlloc, Layout, System}; +use std::sync::atomic::{AtomicBool, AtomicUsize, Ordering}; + +struct Counting; +static COUNTING: AtomicBool = AtomicBool::new(false); +static ALLOCATIONS: AtomicUsize = AtomicUsize::new(0); +unsafe impl GlobalAlloc for Counting { + unsafe fn alloc(&self, layout: Layout) -> *mut u8 { + if COUNTING.load(Ordering::Relaxed) { + ALLOCATIONS.fetch_add(1, Ordering::Relaxed); + } + unsafe { System.alloc(layout) } + } + unsafe fn alloc_zeroed(&self, layout: Layout) -> *mut u8 { + if COUNTING.load(Ordering::Relaxed) { + ALLOCATIONS.fetch_add(1, Ordering::Relaxed); + } + unsafe { System.alloc_zeroed(layout) } + } + unsafe fn realloc(&self, ptr: *mut u8, layout: Layout, new_size: usize) -> *mut u8 { + if COUNTING.load(Ordering::Relaxed) { + ALLOCATIONS.fetch_add(1, Ordering::Relaxed); + } + unsafe { System.realloc(ptr, layout, new_size) } + } + unsafe fn dealloc(&self, ptr: *mut u8, layout: Layout) { + unsafe { System.dealloc(ptr, layout) } + } +} +#[global_allocator] +static GLOBAL: Counting = Counting; + +/// Allocations during `f`, and the number of parallel sections and splits it dispatched. +fn counted(mut f: impl FnMut()) -> (usize, usize) { + f(); // warm-up: first-use allocations (thread locals, lazy statics) are not the region's + stats::enable(true); + let _ = stats::take(); + ALLOCATIONS.store(0, Ordering::SeqCst); + COUNTING.store(true, Ordering::SeqCst); + f(); + COUNTING.store(false, Ordering::SeqCst); + let n = ALLOCATIONS.load(Ordering::SeqCst); + stats::enable(false); + let s = stats::take(); + let parallel = s.run.iter().skip(2).sum::() + s.split.iter().skip(2).sum::(); + (n, parallel) +} + +fn sample(m: usize, n: usize, seed: u64) -> Mat { + let mut x = seed.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407); + Mat::from_fn(m, n, |_, _| { + x = x.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407); + ((x >> 11) as f64 / (1u64 << 53) as f64) - 0.5 + }) +} + +// one test function: the counter is process wide, so nothing else may run while it counts +#[test] +fn parallel_regions_allocate_nothing() { + use faer::linalg::cholesky::lblt::factor as lblt; + use faer::linalg::lu::full_pivoting::factor as flu; + for width in [4usize, 8, 70] { + faer::team::with_team(width, |mut par: Par<'_>| { + // full pivoting LU: the per-step pivot search (`for_each` over column blocks) + let n = 96; + let a0 = sample(n, n, 1); + let mut a = a0.clone(); + let (mut p, mut pi, mut q, mut qi) = (vec![0usize; n], vec![0usize; n], vec![0usize; n], vec![0usize; n]); + let params = flu::FullPivLuParams { par_threshold: ParThreshold::Fixed(0), ..>::auto() }; + let mut mem = MemBuffer::new(flu::lu_in_place_scratch::(n, n, par.shape(), params.into())); + let (allocs, parallel) = counted(|| { + a.copy_from(&a0); + flu::lu_in_place(a.as_mut(), &mut p, &mut pi, &mut q, &mut qi, par.reborrow(), MemStack::new(&mut mem), params.into()); + }); + assert!(parallel > 0, "full pivoting LU, width {width}: no parallel region"); + assert_eq!(allocs, 0, "full pivoting LU, width {width}"); + + // Bunch-Kaufman with full pivoting: rank-one and rank-two steps + let x = sample(n, n, 2); + let s0 = &x + x.transpose(); + let mut s = s0.clone(); + let mut sub = faer::diag::Diag::::zeros(n); + let params = lblt::LbltParams { + pivoting: lblt::PivotingStrategy::Full, + par_threshold: ParThreshold::Fixed(0), + ..>::auto() + }; + let mut mem = MemBuffer::new(lblt::cholesky_in_place_scratch::(n, par.shape(), params.into())); + let (allocs, parallel) = counted(|| { + s.copy_from(&s0); + lblt::cholesky_in_place(s.as_mut(), sub.as_mut(), &mut p, &mut pi, par.reborrow(), MemStack::new(&mut mem), params.into()); + }); + assert!(parallel > 0, "Bunch-Kaufman, width {width}: no parallel region"); + assert_eq!(allocs, 0, "Bunch-Kaufman, width {width}"); + // both update ranks ran: a 2x2 pivot leaves a nonzero subdiagonal entry, a 1x1 a zero + let subs: Vec = (0..n).map(|i| sub[i]).collect(); + assert!(subs.iter().any(|&x| x != 0.0) && subs.iter().any(|&x| x == 0.0), "Bunch-Kaufman, width {width}: a pivot size never occurred"); + + // `join_sized`, the LU swap stage's split, with closures that capture + let mut left = vec![0usize; 4]; + let mut right = vec![0usize; 4]; + let (allocs, parallel) = counted(|| { + faer::utils::thread::join_sized( + par.reborrow(), + width / 2, + |_| left.iter_mut().for_each(|x| *x += 1), + |_| right.iter_mut().for_each(|x| *x += 2), + ); + }); + assert!(parallel > 0, "join_sized, width {width}: no split"); + assert_eq!(allocs, 0, "join_sized, width {width}"); + }); + } +}