diff --git a/PATCHSET.md b/PATCHSET.md index 4921c5f8..ef6efca4 100644 --- a/PATCHSET.md +++ b/PATCHSET.md @@ -124,6 +124,7 @@ gate has changed and two kinds of adoption are now in the ledger. | Serial path overhead of a team token | A QR 128 on a `Par::Team` that cannot pay was +11 % over `Par::Seq`: the 59 block applies asked every product and triangular step. The block apply asks one floor (`utils::cost::par_can_pay_time`, against `min_cost`, the cheapest width of the table) and runs the inner steps sequentially; `kernel_width` rejects against the same cheapest width instead of assuming a table non-increasing in width; gates read the token's table without a `ParShape` copy; the Rayon host table is cached. threadteam pinned at `7095893`. Fork-only (`Par::Team`). | KEEP (with the team executor) | | 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). | ## Ground rules diff --git a/docs/worklogs/2026-10-29-width-per-region.md b/docs/worklogs/2026-10-29-width-per-region.md new file mode 100644 index 00000000..f4245642 --- /dev/null +++ b/docs/worklogs/2026-10-29-width-per-region.md @@ -0,0 +1,88 @@ +# Every region at the widest width that pays (tensor4all/faer#91) + +## Decisions +- **One rule, one helper.** The maintainer and monitor decided the rule: a region runs at the widest width in + `2..=degree` that pays, `pays_at(estimate, sections, dispatch_cost(w))`, on a sub-team of that width. + - `utils::cost::par_width` (units at a rate) and `par_width_time` (a serial time) search the widths from + the full degree down and stop at the first width that pays. Every width is considered, so a falling row + (`aarch64-apple`) is safe. The cheap rejection against the cheapest width of the table comes first. + - `par_width_upto` and `par_width_time_upto` search inside a cap, for splits whose task count is bounded + by their shape. + - `ParThreshold::region_width` adds the parameter policy: `Never` gives 1, and `Fixed(floor)` gives the + whole degree when `units >= floor`. + - `kernel::kernel_width` (the x86 GEMM kernel's two sections) is the same function. +- **Narrowing is a single-child split.** `utils::thread::with_width(par, w, f)` gives: + - `f(par)` when `w >= degree`; + - `f(Par::Seq)` when `w < 2`; + - otherwise `team.split(&[w])`. A one-child split dispatches no worker share: the child is led by the + caller, and threadteam's `section` still runs its join guard and its panic handling. + + The helper lives in `utils::thread`, so callers compile without `threadteam`. +- **Only the region narrows.** Each site wraps only its own region in `with_width`, and the caller keeps its + whole `par` for later, larger work. +- **Per-step factorizations.** Hessenberg, tridiagonal, unblocked bidiagonal, full pivoting LU, Bunch-Kaufman + and column pivoting QR set `par = Par::Seq` the first time a step did not pay and kept it for the rest of + the call. They now decide per step. + - `Never` still makes the whole call sequential, through an explicit `par = Seq` at the top. + - Column pivoting QR's non-delayed product at the head of step k uses the width decided at step k-1, which + is the decision its comment says covers it. + - **Behaviour change:** the post-loop `upgrade_householder_factor` of the Hessenberg and tridiagonal + reductions ran under the last step's latch, which was always sequential by then. It now asks its own + gate. +- **Capped splits.** A split whose task count has an upper bound searches inside the cap; it does not clamp + the answer afterwards, which a non-monotone row would make wrong. + - Triangular solve: jobs of 64 columns. The region's width is searched separately and shared out among the + jobs, which keeps the nested parallelism the old code gave each job. + - The block apply's column split: `K` tasks, likewise. + - `upgrade_householder_factor`: four rows per participant, at both gates. A narrower outer width re-runs + the call on a sub-team of that width. +- **Blocked bidiagonal panel.** The gate stays a floor: sequential, or the whole `par`. Its estimate is the + remaining area, and narrowing the whole panel to it would cap the panel's two trailing products, which + choose their own widths. +- **The tiled product (the non-x86 route) takes the same rule.** `gemm_tiled` chose between the whole + team and one thread. It now asks `par_width_team`. When the tile grid uses fewer participants than the + chosen width, that smaller width must pay as well, because a row need not be monotone. +- **The `V^H V` reduction splits into exactly the chosen number of row ranges.** It uses balanced + `par_split_indices`, where ceiling-sized chunks could produce fewer tasks than the width the gate chose. +- **`Never` at the public boundaries.** Column pivoting QR and the bidiagonalization apply `Never` before + forwarding `par` to their compact-factor upgrades. Otherwise the upgrade's own gate could open a split + under `Never`, now that it caps by rows per participant instead of rejecting. +- **New gate on the Bunch-Kaufman blocked swap loop.** It split over the whole degree with no gate at all, + and the family test caught it. It now uses `par_width` at the LU swap rate. +- **Removed:** + - `par_pays`, `region_pays` and `region_pays_impl`; + - the GEMM route's full-width predicate at `matmul/mod.rs`. It was reached only after the team arm + returns, so it was always 1. + + `ParShape::pays_full` stays as public API. +- **Trace.** `par-trace` records the chosen width (a new `width` field), and `trace_pub` takes it too. +- **Counter.** `utils::team::stats` gains `width`: the regions narrowed to a sub-team, by width, up to 256. + The `run` histogram saturates at 16. + +## Verification conclusions and constraints +- **Tests, no timing:** + - `par_width_takes_the_widest_width_that_pays`: a rising row, a falling row, a cap inside a falling row. + - `with_width_runs_on_a_sub_team_of_that_width`: the participant count of a section run on the narrowed + token, and the parent's width afterwards. + - `the_families_narrow_their_regions_and_match_sequential`, on a team of 8 whose row makes width 4 the + widest that pays. Families: full and partial pivoting LU, column pivoting QR, Hessenberg, tridiagonal, + bidiagonal, Bunch-Kaufman, triangular solve and matvecs. For each: + - regions narrowed to 4, and some of the family's own splits ran at 4; + - no split wider than 4; + - the outputs equal the sequential ones within `1e-10` relative. + - The family test runs every converted region in its own counter window: the triangular solve and both + matvecs separately, Bunch-Kaufman with full pivoting (the default `PartialDiag` does not reach the + converted steps), the block apply, the compact-factor upgrade, the sparse product and the tiled product. + - `never_reaches_the_compact_factor_upgrades`: a 33 x 9 column pivoting QR at block 9 and a + bidiagonalization under `Never`, on a cheap row, open no section. It fails without the boundary guards. + - The rewritten gate tests: `region_rules_follow_the_serial_time` (partial pivoting LU), + `fixed_and_never_keep_their_meaning`, `par_width_reads_the_degree_of_the_token`. +- **Builds:** default, `par-trace`, and `--no-default-features --features std`. The last has the same + warnings as the base. +- **Not changed (follow-ups):** + - The Hessenberg step estimate is `(n-k-1)^2`. The work it gates after the first step is the `A02` matvecs, + `k (n-k-1)`, which grows while the estimate shrinks. Re-modelling it needs calibration. + - The regions with no gate or a degree-wide split that #94 owns: `rank_update`, the SIMD GEMM loops, the + schoolbook product, native `K == 1`. + - The section counts passed by some gates are inherited models, not dispatch counts. +- **Performance evidence is the monitor's:** broccoli at 1/8/16/32/64 threads, and evo. diff --git a/faer/src/linalg/cholesky/bunch_kaufman/factor.rs b/faer/src/linalg/cholesky/bunch_kaufman/factor.rs index d07b9da9..dc2b18b5 100644 --- a/faer/src/linalg/cholesky/bunch_kaufman/factor.rs +++ b/faer/src/linalg/cholesky/bunch_kaufman/factor.rs @@ -1,5 +1,6 @@ use crate::internal_prelude::*; use crate::utils::cost::{ParThreshold, rate}; +use crate::utils::thread::with_width; use crate::{assert, perm}; use linalg::matmul::triangular::BlockStructure; const TOP_BIT: usize = 1 << (usize::BITS - 1); @@ -280,6 +281,10 @@ fn lblt_full_piv( let mut A = A; let mut subdiag = subdiag.column_vector_mut(); let mut par = par; + // `Never` keeps the whole call sequential (the steps decide their own width otherwise) + if matches!(params.par_threshold, ParThreshold::Never) { + par = Par::Seq; + } let n = A.nrows(); let scale_fwd = &A.norm_max(); let scale_bwd = &scale_fwd.recip(); @@ -331,14 +336,8 @@ fn lblt_full_piv( let rem = n - npiv; // One section per step (the rank-1 or rank-2 update fused with the pivot search); // the unit is an element of the trailing block. - if !params.par_threshold.region_pays::( - par.shape(), - rem * rem, - rate::bunch_kaufman(), - 1, - ) { - par = Par::Seq; - } + // The step runs at the widest width that pays. + let width = params.par_threshold.region_width::(&par, rem * rem, rate::bunch_kaufman(), 1); if i0 != 0 { swap_self_adjoint(A.rb_mut(), 0, i0); perm::swap_rows_idx(Aprev.rb_mut(), 0, i0); @@ -359,8 +358,9 @@ fn lblt_full_piv( A[(i, i)] = (A[(i, i)].real() - &diag * L[i].abs2()).to_cplx(); } if n != 0 { - (max_i, max_j, max_offdiag) = - rank_1_update_and_argmax(A.rb_mut(), L.rb(), diag, par.reborrow()); + (max_i, max_j, max_offdiag) = with_width(par.reborrow(), width, |par| { + rank_1_update_and_argmax(A.rb_mut(), L.rb(), diag, par) + }); } } else { let a00 = A[(0, 0)].real(); @@ -380,16 +380,18 @@ fn lblt_full_piv( let (mut L0, mut L1) = L.two_cols_mut(0, 1); let n = A.nrows(); if n != 0 { - (max_i, max_j, max_offdiag) = rank_2_update_and_argmax( - A.rb_mut(), - L0.rb(), - L1.rb(), - d.copy(), - d00.copy(), - d11.copy(), - d10.copy(), - par.reborrow(), - ); + (max_i, max_j, max_offdiag) = with_width(par.reborrow(), width, |par| { + rank_2_update_and_argmax( + A.rb_mut(), + L0.rb(), + L1.rb(), + d.copy(), + d00.copy(), + d11.copy(), + d10.copy(), + par, + ) + }); } for j in 0..n { let x0 = &L0[j].copy(); @@ -761,35 +763,45 @@ fn lblt_blocked( let pivots = &pivots[k..next]; let A_left = A.rb_mut().get_mut(.., ..k); if A_left.ncols() > 0 { - match par { - Par::Seq => { - for mut col in A_left.col_iter_mut() { - for (i, &j) in core::iter::zip(k..next, pivots) { - let j = j & !TOP_BIT; - linalg::lu::partial_pivoting::factor::swap_elems( - col.rb_mut(), - i, - j, - ); - } - } - }, - #[cfg(feature = "rayon")] - _ => { - let nthreads = core::num::NonZeroUsize::new(par.degree()).unwrap(); - let nthreads = nthreads.get(); - crate::utils::thread::for_each(par.reborrow(), - nthreads, - A_left.par_col_iter_mut(), - |mut col| { + // one section over the columns left of the block, at the widest width that pays for + // its row swaps (it was split over the whole degree, whatever its size) + let width = crate::utils::cost::par_width::( + &par, + A_left.ncols().saturating_mul(next - k), + rate::lu_swap(), + 1, + ); + with_width(par.reborrow(), width, |par| { + match par { + Par::Seq => { + for mut col in A_left.col_iter_mut() { for (i, &j) in core::iter::zip(k..next, pivots) { let j = j & !TOP_BIT; - linalg::lu::partial_pivoting::factor::swap_elems(col.rb_mut(), i, j); + linalg::lu::partial_pivoting::factor::swap_elems( + col.rb_mut(), + i, + j, + ); } - }, - ); - }, - } + } + }, + #[cfg(feature = "rayon")] + _ => { + let nthreads = core::num::NonZeroUsize::new(par.degree()).unwrap(); + let nthreads = nthreads.get(); + crate::utils::thread::for_each(par, + nthreads, + A_left.par_col_iter_mut(), + |mut col| { + for (i, &j) in core::iter::zip(k..next, pivots) { + let j = j & !TOP_BIT; + linalg::lu::partial_pivoting::factor::swap_elems(col.rb_mut(), i, j); + } + }, + ); + }, + } + }); } k = next; } diff --git a/faer/src/linalg/evd/hessenberg.rs b/faer/src/linalg/evd/hessenberg.rs index f4ac8d1a..c19b0442 100644 --- a/faer/src/linalg/evd/hessenberg.rs +++ b/faer/src/linalg/evd/hessenberg.rs @@ -1,5 +1,6 @@ use crate::assert; use crate::internal_prelude::*; +use crate::utils::thread::with_width; use crate::utils::cost::{ParThreshold, rate}; use linalg::householder::{self, HouseholderInfo}; use linalg::matmul::triangular::BlockStructure; @@ -249,6 +250,10 @@ fn hessenberg_rearranged_unblocked( let mut A = A; let mut H = H; let mut par = par; + // `Never` keeps the whole call sequential (the steps decide their own width otherwise) + if matches!(params.par_threshold, ParThreshold::Never) { + par = Par::Seq; + } { let mut H = H.rb_mut().row_mut(0); alloca!('stack: { @@ -289,20 +294,14 @@ fn hessenberg_rearranged_unblocked( |uz!(a, u, z): Zip!(&mut T, &T, &T)| *a -= u * &y1 + z, ); } - { - // Two sections per step: the two matrix-vector products (`w0` and the rank-one - // update of `A02`, or `z2` and `y2` on the first step); the unit is an element - // of the trailing block. + // Two sections per step: the two matrix-vector products (`w0` and the rank-one + // update of `A02`, or `z2` and `y2` on the first step); the unit is an element of + // the trailing block. The step's products run at the widest width that pays; the + // next step decides again (`A02` grows while the trailing block shrinks). + let width = { let n = n - k - 1; - if !params.par_threshold.region_pays::( - par.shape(), - n * n, - rate::hessenberg(), - 2, - ) { - par = Par::Seq; - } - } + params.par_threshold.region_width::(&par, n * n, rate::hessenberg(), 2) + }; if k + 1 == n { break; } @@ -337,22 +336,26 @@ fn hessenberg_rearranged_unblocked( y2.copy_from(v2.rb()); z2.copy_from(w2.rb().col(0)); } else { - matmul( - z2.rb_mut().as_mat_mut(), - Accum::Replace, - A22.rb(), - x2.as_mat(), - one(), - par.reborrow(), - ); - matmul( - y2.rb_mut().as_mat_mut(), - Accum::Replace, - x2.adjoint().as_mat(), - A22.rb(), - one(), - par.reborrow(), - ); + with_width(par.reborrow(), width, |par| { + matmul( + z2.rb_mut().as_mat_mut(), + Accum::Replace, + A22.rb(), + x2.as_mat(), + one(), + par, + ); + }); + with_width(par.reborrow(), width, |par| { + matmul( + y2.rb_mut().as_mat_mut(), + Accum::Replace, + x2.adjoint().as_mat(), + A22.rb(), + one(), + par, + ); + }); } let u2 = x2; let ref b = dot::inner_prod( @@ -376,22 +379,26 @@ fn hessenberg_rearranged_unblocked( .mul_real(&tau_inv); z!(&mut A12, u2.transpose()) .for_each(|uz!(a, u): Zip!(&mut T, &T)| *a -= dot * u.conj()); - matmul( - w0.rb_mut().col_mut(0).as_mat_mut(), - Accum::Replace, - A02.rb(), - u2.as_mat(), - one(), - par.reborrow(), - ); - matmul( - A02.rb_mut(), - Accum::Add, - w0.rb().col(0).as_mat(), - u2.adjoint().as_mat(), - (-&tau_inv).to_cplx(), - par.reborrow(), - ); + with_width(par.reborrow(), width, |par| { + matmul( + w0.rb_mut().col_mut(0).as_mat_mut(), + Accum::Replace, + A02.rb(), + u2.as_mat(), + one(), + par, + ); + }); + with_width(par.reborrow(), width, |par| { + matmul( + A02.rb_mut(), + Accum::Add, + w0.rb().col(0).as_mat(), + u2.adjoint().as_mat(), + (-&tau_inv).to_cplx(), + par, + ); + }); A21[0] = beta; } } diff --git a/faer/src/linalg/evd/tridiag.rs b/faer/src/linalg/evd/tridiag.rs index 28d1c85f..db25efd4 100644 --- a/faer/src/linalg/evd/tridiag.rs +++ b/faer/src/linalg/evd/tridiag.rs @@ -1,6 +1,7 @@ use crate::assert; use crate::internal_prelude::*; use crate::utils::cost::{ParThreshold, rate}; +use crate::utils::thread::with_width; use linalg::householder; use linalg::matmul::triangular::BlockStructure; use linalg::matmul::{self, dot}; @@ -285,6 +286,10 @@ pub fn tridiag_in_place( let mut A = A; let mut H = householder; let mut par = par; + // `Never` keeps the whole call sequential (the steps decide their own width otherwise) + if matches!(params.par_threshold, ParThreshold::Never) { + par = Par::Seq; + } let n = A.nrows(); let b = H.nrows(); assert!(H.ncols() == n.saturating_sub(1)); @@ -323,15 +328,8 @@ pub fn tridiag_in_place( } let rem = n - k - 1; // One section per step: the symmetric fused update over column blocks; the unit is - // an element of the triangle. - if !params.par_threshold.region_pays::( - par.shape(), - rem * rem / 2, - rate::tridiag(), - 1, - ) { - par = Par::Seq; - } + // an element of the triangle. The step runs at the widest width that pays. + let width = params.par_threshold.region_width::(&par, rem * rem / 2, rate::tridiag(), 1); let k1 = k + 1; let tau_inv; { @@ -362,128 +360,134 @@ pub fn tridiag_in_place( }, ); w2.copy_from(y2.rb()); - match par { - Par::Seq => { - let mut z2 = z2.rb_mut().col_mut(0); - tridiag_fused_op( - A22.rb_mut(), - y2.rb_mut(), - z2.rb_mut(), - w2.rb(), - w2.rb(), - u2.rb(), - u2.rb(), - x2.rb(), - x2.rb(), - tau_inv.to_cplx(), - simd_align(k1 + 1), - ); - z!(y2.rb_mut(), z2.rb()).for_each( - |uz!(y, z): Zip!(&mut T, &T)| *y += z, - ); - }, - #[cfg(feature = "rayon")] - _ => { - let nthreads = core::num::NonZeroUsize::new(par.degree()).unwrap(); - use rayon::prelude::*; - let nthreads = nthreads.get(); - let mut z2 = z2.rb_mut().subcols_mut(0, nthreads); - let n2 = A22.ncols(); - assert!((n2 as u64) < (1u64 << 50)); - let idx_to_col_start = |idx: usize| { - let idx_as_percent = - idx as f64 / nthreads as f64; - let col_start_percent = 1.0f64 - - libm::sqrt(1.0f64 - idx_as_percent); - (col_start_percent * n2 as f64) as usize - }; - { - let A22 = A22.rb(); - let y2 = y2.rb(); - let ref f = tau_inv.to_cplx::(); - crate::utils::thread::for_each(par.reborrow(), - nthreads, - z2.rb_mut().par_col_iter_mut().enumerate(), - |(idx, mut z2)| { - let first = idx_to_col_start(idx); - let last_col = - idx_to_col_start(idx + 1); - let nrows = n2 - first; - let ncols = last_col - first; - let mut A = unsafe { - A22.rb() - .subcols(first, ncols) - .subrows(first, nrows) - .const_cast() - }; - { - let y2 = unsafe { - y2.subrows(first, ncols) + with_width(par.reborrow(), width, |par| { + match par { + Par::Seq => { + let mut z2 = z2.rb_mut().col_mut(0); + tridiag_fused_op( + A22.rb_mut(), + y2.rb_mut(), + z2.rb_mut(), + w2.rb(), + w2.rb(), + u2.rb(), + u2.rb(), + x2.rb(), + x2.rb(), + tau_inv.to_cplx(), + simd_align(k1 + 1), + ); + z!(y2.rb_mut(), z2.rb()).for_each( + |uz!(y, z): Zip!(&mut T, &T)| *y += z, + ); + }, + #[cfg(feature = "rayon")] + _ => { + let nthreads = core::num::NonZeroUsize::new(par.degree()).unwrap(); + use rayon::prelude::*; + let nthreads = nthreads.get(); + let mut z2 = z2.rb_mut().subcols_mut(0, nthreads); + let n2 = A22.ncols(); + assert!((n2 as u64) < (1u64 << 50)); + let idx_to_col_start = |idx: usize| { + let idx_as_percent = + idx as f64 / nthreads as f64; + let col_start_percent = 1.0f64 + - libm::sqrt(1.0f64 - idx_as_percent); + (col_start_percent * n2 as f64) as usize + }; + { + let A22 = A22.rb(); + let y2 = y2.rb(); + let ref f = tau_inv.to_cplx::(); + crate::utils::thread::for_each(par, + nthreads, + z2.rb_mut().par_col_iter_mut().enumerate(), + |(idx, mut z2)| { + let first = idx_to_col_start(idx); + let last_col = + idx_to_col_start(idx + 1); + let nrows = n2 - first; + let ncols = last_col - first; + let mut A = unsafe { + A22.rb() + .subcols(first, ncols) + .subrows(first, nrows) .const_cast() }; - let mut z2 = z2 - .rb_mut() - .subrows_mut(first, nrows); - let ry2 = - w2.rb().subrows(first, ncols); - let rz2 = - w2.rb().subrows(first, nrows); - let u0 = u2.subrows(first, nrows); - let u1 = u2.subrows(first, ncols); - let u2 = - x2.rb().subrows(first, ncols); - let v2 = - x2.rb().subrows(first, nrows); - tridiag_fused_op( - A.rb_mut(), - y2, - z2.rb_mut(), - ry2, - rz2, - u0, - u1, - u2, - v2, - f.copy(), - n.next_power_of_two() - - (k1 + 1) - first, - ); - } - z2.rb_mut() - .subrows_mut(0, first) - .fill(zero()); - }, - ); - } - for z2 in z2.rb().col_iter() { - z!(y2.rb_mut(), z2.rb()) - .for_each(|uz!(y, z)| *y += z); - } - }, - } + { + let y2 = unsafe { + y2.subrows(first, ncols) + .const_cast() + }; + let mut z2 = z2 + .rb_mut() + .subrows_mut(first, nrows); + let ry2 = + w2.rb().subrows(first, ncols); + let rz2 = + w2.rb().subrows(first, nrows); + let u0 = u2.subrows(first, nrows); + let u1 = u2.subrows(first, ncols); + let u2 = + x2.rb().subrows(first, ncols); + let v2 = + x2.rb().subrows(first, nrows); + tridiag_fused_op( + A.rb_mut(), + y2, + z2.rb_mut(), + ry2, + rz2, + u0, + u1, + u2, + v2, + f.copy(), + n.next_power_of_two() + - (k1 + 1) - first, + ); + } + z2.rb_mut() + .subrows_mut(0, first) + .fill(zero()); + }, + ); + } + for z2 in z2.rb().col_iter() { + z!(y2.rb_mut(), z2.rb()) + .for_each(|uz!(y, z)| *y += z); + } + }, + } + }); } else { - matmul::triangular::matmul( - y2.rb_mut(), - BlockStructure::Rectangular, - Accum::Replace, - A22.rb(), - BlockStructure::TriangularLower, - x2.rb(), - BlockStructure::Rectangular, - tau_inv.to_cplx(), - par.reborrow(), - ); - matmul::triangular::matmul( - y2.rb_mut(), - BlockStructure::Rectangular, - Accum::Add, - A22.rb().adjoint(), - BlockStructure::StrictTriangularUpper, - x2.rb(), - BlockStructure::Rectangular, - tau_inv.to_cplx(), - par.reborrow(), - ); + with_width(par.reborrow(), width, |par| { + matmul::triangular::matmul( + y2.rb_mut(), + BlockStructure::Rectangular, + Accum::Replace, + A22.rb(), + BlockStructure::TriangularLower, + x2.rb(), + BlockStructure::Rectangular, + tau_inv.to_cplx(), + par, + ); + }); + with_width(par.reborrow(), width, |par| { + matmul::triangular::matmul( + y2.rb_mut(), + BlockStructure::Rectangular, + Accum::Add, + A22.rb().adjoint(), + BlockStructure::StrictTriangularUpper, + x2.rb(), + BlockStructure::Rectangular, + tau_inv.to_cplx(), + par, + ); + }); } z!(y2.rb_mut(), A21.rb()).for_each( |uz!(y, a): Zip!(&mut T, &T)| { diff --git a/faer/src/linalg/householder.rs b/faer/src/linalg/householder.rs index 94b812b9..91359c21 100644 --- a/faer/src/linalg/householder.rs +++ b/faer/src/linalg/householder.rs @@ -161,21 +161,31 @@ pub fn upgrade_householder_factor( // small against the basis) and give every thread at least four rows (`4 * threads`). The // third is the cost rule: `deep_rows * n * n` multiply-adds at the GEMM rate must pay for one // section. - let threads = crate::utils::thread::parallelism_degree(&par); + // The width is the widest that pays, with at least four rows per participant; a narrower + // width than the token's reruns this call on a sub-team of that width, where the same + // decision then takes the whole sub-team. let deep_rows = essentials.nrows().unbound().saturating_sub(n); #[cfg(feature = "par-trace")] crate::utils::cost::trace::shape("upgrade_above_floor", deep_rows, n, n); - let above_floor = deep_rows >= 2 * n - && deep_rows >= 4 * threads - && crate::utils::cost::par_pays::( + let width = if deep_rows >= 2 * n { + crate::utils::cost::par_width_upto::( &par, deep_rows.saturating_mul(n).saturating_mul(n), crate::utils::cost::rate::gemm_mac(), 1, - ); + deep_rows / 4, + ) + } else { + 1 + }; #[cfg(feature = "par-trace")] crate::utils::cost::trace::shape("", 0, 0, 0); - let mut par = if above_floor { par } else { Par::Seq }; + if width >= 2 && width < par.degree() { + return crate::utils::thread::with_width(par, width, |par| { + upgrade_householder_factor(householder_factor, essentials, block_size, prev_block_size, par) + }); + } + let mut par = if width >= 2 { par } else { Par::Seq }; let mut householder_factor = householder_factor; let essentials = essentials; assert!(householder_factor.nrows() == householder_factor.ncols()); @@ -232,20 +242,23 @@ pub fn upgrade_householder_factor( // threads for `m = 4096, n = 16`, against 0.040 when its rows are split. The partials here // are `n x n` each, tiny against the basis, and only the strict upper triangle is stored // (the diagonal is the caller's `tau`), so the split costs almost nothing. - let threads = crate::utils::thread::parallelism_degree(&par); let rows = basis_bot.nrows().unbound(); - // shape conditions (`2 * n`, `4 * threads`) as above, then the cost rule + // shape conditions (`2 * n`, four rows per participant) as above, then the widest width + // that pays #[cfg(feature = "par-trace")] crate::utils::cost::trace::shape("upgrade_deep", rows, n, n); - let deep = threads > 1 - && rows >= 2 * n - && rows >= 4 * threads - && crate::utils::cost::par_pays::( + let threads = if rows >= 2 * n { + crate::utils::cost::par_width_upto::( &par, rows.saturating_mul(n).saturating_mul(n), crate::utils::cost::rate::gemm_mac(), 1, - ); + rows / 4, + ) + } else { + 1 + }; + let deep = threads > 1; #[cfg(feature = "par-trace")] crate::utils::cost::trace::shape("", 0, 0, 0); if n <= NARROW_BASIS { @@ -270,17 +283,18 @@ pub fn upgrade_householder_factor( #[cfg(feature = "rayon")] { use rayon::prelude::*; - let chunk = rows.div_ceil(threads); - let n_chunks = rows.div_ceil(chunk); + // exactly `threads` balanced row ranges (at least four rows each), so the split runs + // at the width the gate chose + let n_chunks = threads; let mut partials = alloc::vec![T::zero(); n_chunks * n * n]; let source = basis_bot; - crate::utils::thread::for_each( + crate::utils::thread::with_width(par, threads, |par| crate::utils::thread::for_each( par, n_chunks, partials.par_chunks_mut(n * n).enumerate(), |(index, part)| { - let lo = index * chunk; - let hi = rows.min(lo + chunk); + let (lo, len) = crate::utils::thread::par_split_indices(rows, index, n_chunks); + let hi = lo + len; let mut part = MatMut::from_column_major_slice_mut(part, n, n); matmul_with_conj( @@ -294,7 +308,7 @@ pub fn upgrade_householder_factor( Par::Seq, ); }, - ); + )); for index in 0..n_chunks { let part = MatRef::from_column_major_slice( &partials[index * n * n..(index + 1) * n * n], @@ -545,7 +559,7 @@ fn apply_block_householder_on_the_left_in_place_generic< N.unbound(), matrix.ncols().unbound(), ); - let big_rank_one = crate::utils::cost::par_pays::( + let big_rank_one = crate::utils::cost::par_width::( &par, householder_basis .nrows() @@ -553,7 +567,7 @@ fn apply_block_householder_on_the_left_in_place_generic< .saturating_mul(matrix.ncols().unbound()), crate::utils::cost::rate::householder_rank_one(), 1, - ); + ) > 1; let use_simd = N.unbound() == 1 && T::SIMD_CAPABILITIES.is_simd() && !(crate::utils::thread::parallelism_degree(&par) > 1 && big_rank_one); @@ -681,42 +695,52 @@ fn apply_block_householder_on_the_left_in_place_generic< // region inside each -- `512 x 512` `qr_in_place` f64 at eight threads: 5.27 ms with // the cap, 4.74 ms with eight tasks, and `1024 x 1024`: 20.5 -> 16.9 ms. let mut n_tasks = Ord::min(threads, K.unbound()); + let mut split_width = threads; // The column split is one section over `n_tasks` participants; its work is the reduction // plus the two products of the update, `2 * M * N * K` multiply-adds, plus the streaming // of the `M * K` target // (`cost::block_apply_estimate`). With the multiply-adds alone this estimate was 0.6 to // 0.8 of the measured serial time of an `N = 8` block (`QR 128`). if n_tasks > 1 { - let width = core::num::NonZeroUsize::new(n_tasks).unwrap(); let est = crate::utils::cost::block_apply_estimate::( M.unbound(), N.unbound(), K.unbound(), ); + // the split: the widest number of tasks (at most `K`) that pays for its one section; + // the region: the widest width that pays, shared out among the tasks + let tasks = crate::utils::cost::par_width_time_upto(&par, est, 1, n_tasks); + let width = crate::utils::cost::par_width_time(&par, est, 1); #[cfg(feature = "par-trace")] #[cfg(feature = "rayon")] { crate::utils::cost::trace::shape("n_tasks", M.unbound(), N.unbound(), K.unbound()); + let w = core::num::NonZeroUsize::new(width.max(2)).unwrap(); crate::utils::cost::trace_pub( M.unbound() * N.unbound() * K.unbound(), 1, est, - par.shape().dispatch_cost(width), - par.shape().pays(est, 1, width), + par.shape().dispatch_cost(w), + width, + width > 1, core::panic::Location::caller(), ); } - let pays = par.shape().pays(est, 1, width); #[cfg(feature = "par-trace")] - let pays = match crate::utils::cost::trace::N_TASKS_MODE.load(core::sync::atomic::Ordering::Relaxed) { - 1 => (M.unbound() * K.unbound()).saturating_mul(4 * M.unbound()) >= gemm::get_threading_threshold(), - 2 => true, - 3 => false, - _ => pays, + let (tasks, width) = match crate::utils::cost::trace::N_TASKS_MODE.load(core::sync::atomic::Ordering::Relaxed) { + 1 => { + if (M.unbound() * K.unbound()).saturating_mul(4 * M.unbound()) >= gemm::get_threading_threshold() { + (n_tasks, threads) + } else { + (1, 1) + } + }, + 2 => (n_tasks, threads), + 3 => (1, 1), + _ => (tasks, width), }; - if !pays { - n_tasks = 1; - } + n_tasks = if width > 1 { tasks } else { 1 }; + split_width = width; } let func = |(mut tmp, mut matrix): ( MatMut<'_, T, Dim<'N>>, @@ -796,17 +820,20 @@ fn apply_block_householder_on_the_left_in_place_generic< #[cfg(feature = "rayon")] { use rayon::prelude::*; - crate::utils::thread::for_each_with_par(par, - n_tasks, - tmp.rb_mut() - .par_col_partition_mut(n_tasks) - .zip_eq(matrix.rb_mut().par_col_partition_mut(n_tasks)), - func, - ); + crate::utils::thread::with_width(par, split_width, |par| { + crate::utils::thread::for_each_with_par( + par, + n_tasks, + tmp.rb_mut() + .par_col_partition_mut(n_tasks) + .zip_eq(matrix.rb_mut().par_col_partition_mut(n_tasks)), + func, + ) + }); } #[cfg(not(feature = "rayon"))] { - let _ = (n_tasks, func); + let _ = (n_tasks, split_width, func); unreachable!("no parallel route in this build") } } diff --git a/faer/src/linalg/lu/full_pivoting/factor.rs b/faer/src/linalg/lu/full_pivoting/factor.rs index 1ffa9d41..5e5d1d02 100644 --- a/faer/src/linalg/lu/full_pivoting/factor.rs +++ b/faer/src/linalg/lu/full_pivoting/factor.rs @@ -2,6 +2,7 @@ use crate::internal_prelude::*; use crate::utils::cost::{ParThreshold, rate}; use crate::perm::{swap_cols_idx, swap_rows_idx}; use crate::utils::thread::par_split_indices; +use crate::utils::thread::with_width; use faer_traits::{Real, RealReg}; use linalg::matmul::matmul; use pulp::Simd; @@ -352,6 +353,10 @@ fn lu_in_place_unblocked( return 0; } let mut par = par; + // `Never` keeps the whole call sequential (the steps decide their own width otherwise) + if matches!(params.par_threshold, ParThreshold::Never) { + par = Par::Seq; + } let mut A = A; let (mut max_row, mut max_col, mut max_score) = best_in_matrix(A.rb()); for k in 0..Ord::min(m, n) { @@ -390,68 +395,68 @@ fn lu_in_place_unblocked( } // One section per step (the fused rank-one update and pivot search over column // blocks); the unit is an element of the trailing block. - if !params.par_threshold.region_pays::( - par.shape(), + let width = params.par_threshold.region_width::( + &par, (m - k - 1) * (n - k - 1), rate::lu_full_pivot(), 1, - ) { - par = Par::Seq; - } + ); let (_, A01, A10, mut A11) = A.rb_mut().split_at_mut(k + 1, k + 1); let lhs = A10.col(k); let rhs = A01.row(k); - match par { - Par::Seq => { - (max_row, max_col, max_score) = - rank_one_update_and_best_in_matrix( - A11.rb_mut(), - lhs, - rhs, - simd_align(k + 1), + with_width(par.reborrow(), width, |par| { + match par { + Par::Seq => { + (max_row, max_col, max_score) = + rank_one_update_and_best_in_matrix( + A11.rb_mut(), + lhs, + rhs, + simd_align(k + 1), + ); + }, + #[cfg(feature = "rayon")] + _ => { + 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::>(); + let full_cols = A11.ncols(); + crate::utils::thread::for_each(par, + nthreads, + best.par_iter_mut() + .zip_eq(A11.rb_mut().par_col_partition_mut(nthreads)) + .zip_eq(rhs.par_partition(nthreads)) + .enumerate(), + |(idx, (((max_row, max_col, max_score), A11), rhs))| { + (*max_row, *max_col, *max_score) = { + let (a, mut b, c) = + rank_one_update_and_best_in_matrix( + A11, + lhs, + rhs, + simd_align(k + 1), + ); + b += par_split_indices(full_cols, idx, nthreads).0; + (a, b, c) + }; + }, ); - }, - #[cfg(feature = "rayon")] - _ => { - 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::>(); - let full_cols = A11.ncols(); - crate::utils::thread::for_each(par.reborrow(), - nthreads, - best.par_iter_mut() - .zip_eq(A11.rb_mut().par_col_partition_mut(nthreads)) - .zip_eq(rhs.par_partition(nthreads)) - .enumerate(), - |(idx, (((max_row, max_col, max_score), A11), rhs))| { - (*max_row, *max_col, *max_score) = { - let (a, mut b, c) = - rank_one_update_and_best_in_matrix( - A11, - lhs, - rhs, - simd_align(k + 1), - ); - b += par_split_indices(full_cols, idx, nthreads).0; - (a, b, c) - }; - }, - ); - 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; + 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; + } } - } - }, - } + }, + } + }); max_row += k + 1; max_col += k + 1; } diff --git a/faer/src/linalg/lu/partial_pivoting/factor.rs b/faer/src/linalg/lu/partial_pivoting/factor.rs index ba25b873..7eee7887 100644 --- a/faer/src/linalg/lu/partial_pivoting/factor.rs +++ b/faer/src/linalg/lu/partial_pivoting/factor.rs @@ -95,30 +95,32 @@ const SWAP_STAGE_SECTIONS: u32 = 4; /// Sections of the swap stage when one side holds every column (the first panel of a wide /// matrix, the last of a tall one): one split of the columns over the whole budget. const SWAP_ONE_SIDE_SECTIONS: u32 = 1; -use crate::utils::cost::{gemm_ns_per_mac, rate, serial_time}; +use crate::utils::cost::{gemm_ns_per_mac, par_width_time, rate, serial_time}; +use crate::utils::thread::with_width; pub use crate::utils::cost::ParThreshold; -/// Whether the trailing update (`b x cols` solve, `rows x cols x b` product) of a level runs -/// parallel. [`ParThreshold::Derived`] derives the answer from the serial time and the -/// dispatch cost of `shape`; [`Fixed`](ParThreshold::Fixed) is an explicit floor in `f64` -/// elements of the updated block, counted in flops (a complex element carries four times the -/// arithmetic). -fn trailing_update_pays( - shape: ParShape, +/// The width the trailing update (`b x cols` solve, `rows x cols x b` product) of a level runs +/// at (1: sequential). [`ParThreshold::Derived`] is the widest width of `par` that pays for the +/// serial time; [`Fixed`](ParThreshold::Fixed) is an explicit floor in `f64` elements of the +/// updated block, counted in flops (a complex element carries four times the arithmetic), and +/// takes the whole degree above it. +fn trailing_update_width( + par: &Par, par_threshold: ParThreshold, rows: usize, cols: usize, b: usize, -) -> bool { +) -> usize { let floor = match par_threshold { - ParThreshold::Never => return false, + ParThreshold::Never => return 1, ParThreshold::Fixed(n) => n, ParThreshold::Derived => { let macs = cols.saturating_mul( b.saturating_mul(b) / 2 + rows.saturating_mul(b), ); - return shape.pays_full( + return par_width_time( + par, serial_time(macs, gemm_ns_per_mac::()), TRAILING_UPDATE_SECTIONS, ); @@ -127,24 +129,27 @@ fn trailing_update_pays( let work = rows .saturating_mul(cols) .saturating_mul(if T::IS_REAL { 2 } else { 8 }); - work > floor.saturating_mul(2) + if work > floor.saturating_mul(2) { par.degree() } else { 1 } } -/// Whether the swap stage (`n` row swaps over `cols` columns of an `m` row matrix) runs -/// parallel; `par_threshold` as in [`trailing_update_pays`], where an explicit value counts -/// `m * cols` elements. -fn swap_stage_pays( - shape: ParShape, +/// The width the swap stage (`n` row swaps over `cols` columns of an `m` row matrix) runs at; +/// `par_threshold` as in [`trailing_update_width`], where an explicit value counts `m * cols` +/// elements. +fn swap_stage_width( + par: &Par, par_threshold: ParThreshold, m: usize, n: usize, cols: usize, one_side: bool, -) -> bool { +) -> usize { match par_threshold { - ParThreshold::Never => false, - ParThreshold::Fixed(floor) => m.saturating_mul(cols) > floor, - ParThreshold::Derived => shape.pays_full( + ParThreshold::Never => 1, + ParThreshold::Fixed(floor) => { + if m.saturating_mul(cols) > floor { par.degree() } else { 1 } + }, + ParThreshold::Derived => par_width_time( + par, serial_time(n.saturating_mul(cols), rate::lu_swap().of::()), if one_side { SWAP_ONE_SIDE_SECTIONS } else { SWAP_STAGE_SECTIONS }, ), @@ -195,33 +200,25 @@ pub(crate) fn lu_in_place_recursion( // cold caches. Below it take `Par::Seq`; the recursion keeps the caller's // `par`, so every level decides for its own trailing update. let (rows1, cols1, b) = (A11.nrows(), A11.ncols(), A10.ncols()); - let mut par_here = if trailing_update_pays::( - par.shape(), - params.par_threshold, - rows1, - cols1, - b, - ) { - par.reborrow() - } else { - Par::Seq - }; + let width = trailing_update_width::(&par, params.par_threshold, rows1, cols1, b); + with_width(par.reborrow(), width, |mut par_here| { #[cfg(feature = "rayon")] - { - linalg::triangular_solve::solve_unit_lower_triangular_in_place( - A00.rb(), - A01.rb_mut(), - par_here.reborrow(), + { + linalg::triangular_solve::solve_unit_lower_triangular_in_place( + A00.rb(), + A01.rb_mut(), + par_here.reborrow(), + ); + } + linalg::matmul::matmul( + A11.rb_mut(), + Accum::Add, + A10.rb(), + A01.rb(), + -one::(), + par_here, ); - } - linalg::matmul::matmul( - A11.rb_mut(), - Accum::Add, - A10.rb(), - A01.rb(), - -one::(), - par_here, - ); + }); n_trans += lu_in_place_recursion( A.rb_mut().get_mut(block_size..m, ..), block_size, @@ -248,70 +245,68 @@ pub(crate) fn lu_in_place_recursion( }; let (A_left, A_right) = A.rb_mut().split_at_col_mut(start); let A_right = A_right.get_mut(.., end - start..ncols - start); - let par = if swap_stage_pays::( - par.shape(), + let width = swap_stage_width::( + &par, params.par_threshold, m, n, ncols - n, start == 0 || end == ncols, - ) { - par - } else { - Par::Seq - }; - match par { - Par::Seq => { - swap(A_left); - swap(A_right); - }, - #[cfg(feature = "rayon")] - _ => { - let nthreads = core::num::NonZeroUsize::new(par.degree()).unwrap(); - let nthreads = nthreads.get(); - let len = (A_left.ncols() + A_right.ncols()) as f64; - let left_threads = Ord::min( - (nthreads as f64 * (A_left.ncols() as f64 / len)) as usize, - nthreads, - ); - let (left_cols, right_cols) = (A_left.ncols(), A_right.ncols()); - if left_cols == 0 || right_cols == 0 { - // one side holds every column: it takes the whole budget - let (A, _) = if left_cols == 0 { (A_right, A_left) } else { (A_left, A_right) }; - crate::utils::thread::for_each( - par, + ); + with_width(par.reborrow(), width, |par| { + match par { + Par::Seq => { + swap(A_left); + swap(A_right); + }, + #[cfg(feature = "rayon")] + _ => { + let nthreads = core::num::NonZeroUsize::new(par.degree()).unwrap(); + let nthreads = nthreads.get(); + let len = (A_left.ncols() + A_right.ncols()) as f64; + let left_threads = Ord::min( + (nthreads as f64 * (A_left.ncols() as f64 / len)) as usize, nthreads, - A.par_col_partition_mut(nthreads), - |A| swap(A), ); - } else { - // the two sides run on the two halves of the split, each on the - // parallelism it is handed - crate::utils::thread::join_sized( - par, - left_threads, - |par| { - let n = par.degree(); - crate::utils::thread::for_each( - par, - n, - A_left.par_col_partition_mut(n), - |A| swap(A), - ) - }, - |par| { - let n = par.degree(); - crate::utils::thread::for_each( - par, - n, - A_right.par_col_partition_mut(n), - |A| swap(A), - ) - }, - ); - } - }, - } + let (left_cols, right_cols) = (A_left.ncols(), A_right.ncols()); + if left_cols == 0 || right_cols == 0 { + // one side holds every column: it takes the whole budget + let (A, _) = if left_cols == 0 { (A_right, A_left) } else { (A_left, A_right) }; + crate::utils::thread::for_each( + par, + nthreads, + A.par_col_partition_mut(nthreads), + |A| swap(A), + ); + } else { + // the two sides run on the two halves of the split, each on the + // parallelism it is handed + crate::utils::thread::join_sized( + par, + left_threads, + |par| { + let n = par.degree(); + crate::utils::thread::for_each( + par, + n, + A_left.par_col_partition_mut(n), + |A| swap(A), + ) + }, + |par| { + let n = par.degree(); + crate::utils::thread::for_each( + par, + n, + A_right.par_col_partition_mut(n), + |A| swap(A), + ) + }, + ); + } + }, + } + }); n_trans } /// $LU$ factorization tuning parameters @@ -467,36 +462,43 @@ mod tests { #[test] fn region_rules_follow_the_serial_time() { use core::num::NonZeroUsize; - let seq = ParShape::Seq; - // sequential never pays, whatever the size - assert!(!trailing_update_pays::(seq, ParThreshold::Derived, 4096, 4096, 256)); - assert!(!swap_stage_pays::(seq, ParThreshold::Derived, 4096, 256, 4096, false)); - { - let par = ParShape::team(6); - // tiny regions stay serial, huge ones go parallel - assert!(!trailing_update_pays::(par, ParThreshold::Derived, 16, 16, 8)); - assert!(trailing_update_pays::(par, ParThreshold::Derived, 2048, 2048, 128)); - assert!(!swap_stage_pays::(par, ParThreshold::Derived, 64, 16, 48, false)); - assert!(swap_stage_pays::(par, ParThreshold::Derived, 4096, 256, 4096, false)); - // the swap stage of a 64 x 1024 matrix (64 swaps over 960 columns, every column on - // one side: one section) follows the rule on the tables of the host class, for both - // dtypes; the same region counted as the two sided split (four sections) is what kept - // it serial. The verdict itself depends on the class, so the test states the rule. - for w in [4, 8] { - let par = ParShape::team(w); - let rule = |ns: f64| pays_at(serial_time(64 * 960, ns), SWAP_ONE_SIDE_SECTIONS, par.dispatch_cost(NonZeroUsize::new(w).unwrap())); - assert_eq!(swap_stage_pays::(par, ParThreshold::Derived, 64, 64, 960, true), rule(rate::lu_swap().of::()), "f64 {w}"); - assert_eq!(swap_stage_pays::(par, ParThreshold::Derived, 64, 64, 960, true), rule(rate::lu_swap().of::()), "c64 {w}"); - } + let seq = Par::Seq; + // sequential never goes wider than one, whatever the size + assert_eq!(trailing_update_width::(&seq, ParThreshold::Derived, 4096, 4096, 256), 1); + assert_eq!(swap_stage_width::(&seq, ParThreshold::Derived, 4096, 256, 4096, false), 1); + crate::team::with_team(6, |par| { + // tiny regions stay serial, huge ones take the whole team + assert_eq!(trailing_update_width::(&par, ParThreshold::Derived, 16, 16, 8), 1); + assert_eq!(trailing_update_width::(&par, ParThreshold::Derived, 2048, 2048, 128), 6); + assert_eq!(swap_stage_width::(&par, ParThreshold::Derived, 64, 16, 48, false), 1); + assert_eq!(swap_stage_width::(&par, ParThreshold::Derived, 4096, 256, 4096, false), 6); // a complex trailing update reaches the floor at fewer elements than a real one let (r, c) = (160, 160); - let real = trailing_update_pays::(par, ParThreshold::Derived, r, c, 32); - let cplx = trailing_update_pays::(par, ParThreshold::Derived, r, c, 32); + let real = trailing_update_width::(&par, ParThreshold::Derived, r, c, 32); + let cplx = trailing_update_width::(&par, ParThreshold::Derived, r, c, 32); assert!(cplx >= real); - // an explicit floor restores the fixed rule and `usize::MAX` disables - assert!(trailing_update_pays::(par, ParThreshold::Fixed(1), 64, 64, 8)); - assert!(!trailing_update_pays::(par, ParThreshold::Never, 4096, 4096, 256)); - assert!(!swap_stage_pays::(par, ParThreshold::Never, 4096, 256, 4096, false)); + // an explicit floor restores the fixed rule (the whole team) and `usize::MAX` disables + assert_eq!(trailing_update_width::(&par, ParThreshold::Fixed(1), 64, 64, 8), 6); + assert_eq!(trailing_update_width::(&par, ParThreshold::Never, 4096, 4096, 256), 1); + assert_eq!(swap_stage_width::(&par, ParThreshold::Never, 4096, 256, 4096, false), 1); + }); + // the swap stage of a 64 x 1024 matrix (64 swaps over 960 columns, every column on one + // side: one section) takes the widest width at which the rule holds on the tables of the + // host class, for both dtypes; the same region counted as the two sided split (four + // sections) is what kept it serial. The verdict depends on the class, so the test states + // the rule. + for w in [4, 8] { + crate::team::with_team(w, |par| { + let shape = par.shape(); + let rule = |ns: f64| { + (2..=w) + .rev() + .find(|&k| pays_at(serial_time(64 * 960, ns), SWAP_ONE_SIDE_SECTIONS, shape.dispatch_cost(NonZeroUsize::new(k).unwrap()))) + .unwrap_or(1) + }; + assert_eq!(swap_stage_width::(&par, ParThreshold::Derived, 64, 64, 960, true), rule(rate::lu_swap().of::()), "f64 {w}"); + assert_eq!(swap_stage_width::(&par, ParThreshold::Derived, 64, 64, 960, true), rule(rate::lu_swap().of::()), "c64 {w}"); + }); } } diff --git a/faer/src/linalg/matmul/internal/mod.rs b/faer/src/linalg/matmul/internal/mod.rs index 72de4ec5..7fa61b35 100644 --- a/faer/src/linalg/matmul/internal/mod.rs +++ b/faer/src/linalg/matmul/internal/mod.rs @@ -78,17 +78,8 @@ pub fn spicy_matmul( let team_gate_pending = matches!(&par, Par::Team(_)) && par.degree() >= 2; #[cfg(not(feature = "threadteam"))] let team_gate_pending = false; - let mut par = if team_gate_pending - || crate::utils::cost::par_pays::( - &par, - units, - crate::utils::cost::rate::gemm_mac(), - 1, - ) { - par - } else { - Par::Seq - }; + // a token that is not a team is sequential: nothing to decide before the kernel + let mut par = par; if let Some(row_idx) = row_idx { for &i in row_idx { assert!(i.zx() < C.nrows()); @@ -284,15 +275,23 @@ pub fn spicy_matmul( return; } } - // no kernel for a team here: the one-section gate of the generic path - if team_gate_pending - && !crate::utils::cost::par_pays::( + // no kernel for a team here: the one-section gate of the generic path, at the widest width + // that pays; a narrower width reruns the call on a sub-team of that width, where the same + // decision then takes the whole sub-team + if team_gate_pending { + let width = crate::utils::cost::par_width::( &par, units, crate::utils::cost::rate::gemm_mac(), 1, - ) { - par = Par::Seq; + ); + if width < par.degree() { + return crate::utils::thread::with_width(par, width, |par| { + spicy_matmul( + C, C_block, row_idx, col_idx, beta, A, conj_A, B, conj_B, D, alpha, par, stack, + ) + }); + } } let (mut out, stack) = unsafe { temp_mat_uninit::( diff --git a/faer/src/linalg/matmul/mod.rs b/faer/src/linalg/matmul/mod.rs index b6f636d3..5213db2b 100644 --- a/faer/src/linalg/matmul/mod.rs +++ b/faer/src/linalg/matmul/mod.rs @@ -773,102 +773,72 @@ mod matvec_rowmajor { // One section (a split of the rows or columns); the unit is a matrix element. The // product reads each element once, so the rate is the streaming one. let size = *lhs.nrows() * *lhs.ncols(); - let par = if crate::utils::cost::par_pays::( + let width = crate::utils::cost::par_width::( &par, size, crate::utils::cost::rate::matvec(), 1, - ) { - par - } else { - Par::Seq - }; - match par { - Par::Seq => { - pub struct Impl<'a, 'M, 'K, T: ComplexField> { - dst: ColMut<'a, T, Dim<'M>>, - beta: Accum, - lhs: MatRef<'a, T, Dim<'M>, Dim<'K>, isize, ContiguousFwd>, - conj_lhs: Conj, - rhs: ColRef<'a, T, Dim<'K>, ContiguousFwd>, - conj_rhs: Conj, - alpha: &'a T, - } - impl<'a, 'M, 'K, T: ComplexField> pulp::WithSimd for Impl<'a, 'M, 'K, T> { - type Output = (); + ); + crate::utils::thread::with_width(par, width, |par| { + match par { + Par::Seq => { + pub struct Impl<'a, 'M, 'K, T: ComplexField> { + dst: ColMut<'a, T, Dim<'M>>, + beta: Accum, + lhs: MatRef<'a, T, Dim<'M>, Dim<'K>, isize, ContiguousFwd>, + conj_lhs: Conj, + rhs: ColRef<'a, T, Dim<'K>, ContiguousFwd>, + conj_rhs: Conj, + alpha: &'a T, + } + impl<'a, 'M, 'K, T: ComplexField> pulp::WithSimd for Impl<'a, 'M, 'K, T> { + type Output = (); - #[inline(always)] - fn with_simd(self, simd: S) -> Self::Output { - let Self { - dst, - beta, - lhs, - conj_lhs, - rhs, - conj_rhs, - alpha, - } = self; - let simd = T::simd_ctx(simd); - let mut dst = dst; - let K = lhs.ncols(); - let simd = SimdCtx::new(simd, K); - for i in lhs.nrows().indices() { - let dst = &mut dst[i]; - let lhs = lhs.row(i); - let rhs = rhs; - let mut tmp = if conj_lhs == conj_rhs { - dot::inner_prod_no_conj_simd::( - simd, - lhs.transpose(), - rhs, - ) - } else { - dot::inner_prod_conj_lhs_simd::( - simd, - lhs.transpose(), - rhs, - ) - }; - if conj_rhs == Conj::Yes { - tmp = tmp.conj(); - } - tmp = alpha * tmp; - if let Accum::Add = beta { - tmp = &*dst + tmp; + #[inline(always)] + fn with_simd(self, simd: S) -> Self::Output { + let Self { + dst, + beta, + lhs, + conj_lhs, + rhs, + conj_rhs, + alpha, + } = self; + let simd = T::simd_ctx(simd); + let mut dst = dst; + let K = lhs.ncols(); + let simd = SimdCtx::new(simd, K); + for i in lhs.nrows().indices() { + let dst = &mut dst[i]; + let lhs = lhs.row(i); + let rhs = rhs; + let mut tmp = if conj_lhs == conj_rhs { + dot::inner_prod_no_conj_simd::( + simd, + lhs.transpose(), + rhs, + ) + } else { + dot::inner_prod_conj_lhs_simd::( + simd, + lhs.transpose(), + rhs, + ) + }; + if conj_rhs == Conj::Yes { + tmp = tmp.conj(); + } + tmp = alpha * tmp; + if let Accum::Add = beta { + tmp = &*dst + tmp; + } + *dst = tmp; } - *dst = tmp; } } - } - dispatch!( - Impl { - dst, - beta, - lhs, - conj_lhs, - rhs, - conj_rhs, - alpha, - }, - Impl, - T - ); - }, - #[cfg(feature = "rayon")] - _ => { - let nthreads = core::num::NonZeroUsize::new(par.degree()).unwrap(); - let nthreads = nthreads.get(); - use rayon::prelude::*; - crate::utils::thread::for_each(par, - nthreads, - dst.par_partition_mut(nthreads) - .zip_eq(lhs.par_row_partition(nthreads)), - |(dst, lhs)| { - make_guard!(M); - let nrows = dst.nrows().bind(M); - let dst = dst.as_row_shape_mut(nrows); - let lhs = lhs.as_row_shape(nrows); - matvec( + dispatch!( + Impl { dst, beta, lhs, @@ -876,12 +846,40 @@ mod matvec_rowmajor { rhs, conj_rhs, alpha, - Par::Seq, - ); - }, - ) - }, - } + }, + Impl, + T + ); + }, + #[cfg(feature = "rayon")] + _ => { + let nthreads = core::num::NonZeroUsize::new(par.degree()).unwrap(); + let nthreads = nthreads.get(); + use rayon::prelude::*; + crate::utils::thread::for_each(par, + nthreads, + dst.par_partition_mut(nthreads) + .zip_eq(lhs.par_row_partition(nthreads)), + |(dst, lhs)| { + make_guard!(M); + let nrows = dst.nrows().bind(M); + let dst = dst.as_row_shape_mut(nrows); + let lhs = lhs.as_row_shape(nrows); + matvec( + dst, + beta, + lhs, + conj_lhs, + rhs, + conj_rhs, + alpha, + Par::Seq, + ); + }, + ) + }, + } + }); } } mod matvec_colmajor { @@ -904,157 +902,155 @@ mod matvec_colmajor { // One section (a split of the rows or columns); the unit is a matrix element. The // product reads each element once, so the rate is the streaming one. let size = *lhs.nrows() * *lhs.ncols(); - let par = if crate::utils::cost::par_pays::( + let width = crate::utils::cost::par_width::( &par, size, crate::utils::cost::rate::matvec(), 1, - ) { - par - } else { - Par::Seq - }; - match par { - Par::Seq => { - pub struct Impl<'a, 'M, 'K, T: ComplexField> { - dst: ColMut<'a, T, Dim<'M>, ContiguousFwd>, - beta: Accum, - lhs: MatRef<'a, T, Dim<'M>, Dim<'K>, ContiguousFwd, isize>, - conj_lhs: Conj, - rhs: ColRef<'a, T, Dim<'K>>, - conj_rhs: Conj, - alpha: &'a T, - } - impl<'a, 'M, 'K, T: ComplexField> pulp::WithSimd for Impl<'a, 'M, 'K, T> { - type Output = (); + ); + crate::utils::thread::with_width(par, width, |par| { + match par { + Par::Seq => { + pub struct Impl<'a, 'M, 'K, T: ComplexField> { + dst: ColMut<'a, T, Dim<'M>, ContiguousFwd>, + beta: Accum, + lhs: MatRef<'a, T, Dim<'M>, Dim<'K>, ContiguousFwd, isize>, + conj_lhs: Conj, + rhs: ColRef<'a, T, Dim<'K>>, + conj_rhs: Conj, + alpha: &'a T, + } + impl<'a, 'M, 'K, T: ComplexField> pulp::WithSimd for Impl<'a, 'M, 'K, T> { + type Output = (); - #[inline(always)] - fn with_simd(self, simd: S) -> Self::Output { - let Self { - dst, - beta, - lhs, - conj_lhs, - rhs, - conj_rhs, - alpha, - } = self; - let simd = T::simd_ctx(simd); - let M = lhs.nrows(); - let simd = SimdCtx::::new(simd, M); - let indices = simd.indices(); - let mut dst = dst; - match beta { - Accum::Add => {}, - Accum::Replace => { - simd_iter!(for i in [indices] { - simd.write(dst.rb_mut(), i, simd.zero()); - }); - }, - } - for j in lhs.ncols().indices() { - let mut dst = dst.rb_mut(); - let lhs = lhs.col(j); - let rhs = &rhs[j]; - let rhs = conj_rhs.apply_rt(rhs); - let rhs = rhs * alpha; - let vrhs = simd.splat(&rhs); - if conj_lhs == Conj::Yes { - simd_iter!(for i in [indices] { - let y = simd.read(dst.rb(), i); - let x = simd.read(lhs, i); - simd.write( - dst.rb_mut(), - i, - simd.conj_mul_add(x, vrhs, y), - ); - }); - } else { - simd_iter!(for i in [indices] { - let y = simd.read(dst.rb(), i); - let x = simd.read(lhs, i); - simd.write( - dst.rb_mut(), - i, - simd.mul_add(x, vrhs, y), - ); - }); + #[inline(always)] + fn with_simd(self, simd: S) -> Self::Output { + let Self { + dst, + beta, + lhs, + conj_lhs, + rhs, + conj_rhs, + alpha, + } = self; + let simd = T::simd_ctx(simd); + let M = lhs.nrows(); + let simd = SimdCtx::::new(simd, M); + let indices = simd.indices(); + let mut dst = dst; + match beta { + Accum::Add => {}, + Accum::Replace => { + simd_iter!(for i in [indices] { + simd.write(dst.rb_mut(), i, simd.zero()); + }); + }, + } + for j in lhs.ncols().indices() { + let mut dst = dst.rb_mut(); + let lhs = lhs.col(j); + let rhs = &rhs[j]; + let rhs = conj_rhs.apply_rt(rhs); + let rhs = rhs * alpha; + let vrhs = simd.splat(&rhs); + if conj_lhs == Conj::Yes { + simd_iter!(for i in [indices] { + let y = simd.read(dst.rb(), i); + let x = simd.read(lhs, i); + simd.write( + dst.rb_mut(), + i, + simd.conj_mul_add(x, vrhs, y), + ); + }); + } else { + simd_iter!(for i in [indices] { + let y = simd.read(dst.rb(), i); + let x = simd.read(lhs, i); + simd.write( + dst.rb_mut(), + i, + simd.mul_add(x, vrhs, y), + ); + }); + } } } } - } - dispatch!( - Impl { - dst, - lhs, - conj_lhs, - rhs, - conj_rhs, - beta, - alpha, - }, - Impl, - T - ) - }, - #[cfg(feature = "rayon")] - _ => { - let nthreads = core::num::NonZeroUsize::new(par.degree()).unwrap(); - use rayon::prelude::*; - let nthreads = nthreads.get(); - let mut mem = MemBuffer::new(temp_mat_scratch::( - dst.nrows().unbound(), - nthreads, - )); - let stack = MemStack::new(&mut mem); - alloca!('stack: { - let tmp = - unsafe { mat![uninit::, dst.nrows(), nthreads] }; - }); - let mut tmp = tmp.try_as_col_major_mut().unwrap(); - let mut dst = dst; - make_guard!(Z); - let Z = 0usize.bind(Z); - let z = IdxInc::new_checked(0, lhs.ncols()); - crate::utils::thread::for_each(par, - nthreads, - tmp.rb_mut() - .par_col_iter_mut() - .zip_eq(lhs.par_col_partition(nthreads)) - .zip_eq(rhs.par_partition(nthreads)), - |((dst, lhs), rhs)| { - make_guard!(K); - let K = lhs.ncols().bind(K); - let lhs = lhs.as_col_shape(K); - let rhs = rhs.as_row_shape(K); - matvec( + dispatch!( + Impl { dst, - Accum::Replace, lhs, conj_lhs, rhs, conj_rhs, + beta, alpha, - Par::Seq, - ); - }, - ); - matvec( - dst.rb_mut(), - beta, - lhs.subcols(z, Z), - conj_lhs, - rhs.subrows(z, Z), - conj_rhs, - &zero(), - Par::Seq, - ); - for j in 0..nthreads { - zip!(dst.rb_mut(), tmp.rb().col(j)) - .for_each(|unzip!(dst, src)| *dst = &*dst + &*src) - } - }, - } + }, + Impl, + T + ) + }, + #[cfg(feature = "rayon")] + _ => { + let nthreads = core::num::NonZeroUsize::new(par.degree()).unwrap(); + use rayon::prelude::*; + let nthreads = nthreads.get(); + let mut mem = MemBuffer::new(temp_mat_scratch::( + dst.nrows().unbound(), + nthreads, + )); + let stack = MemStack::new(&mut mem); + alloca!('stack: { + let tmp = + unsafe { mat![uninit::, dst.nrows(), nthreads] }; + }); + let mut tmp = tmp.try_as_col_major_mut().unwrap(); + let mut dst = dst; + make_guard!(Z); + let Z = 0usize.bind(Z); + let z = IdxInc::new_checked(0, lhs.ncols()); + crate::utils::thread::for_each(par, + nthreads, + tmp.rb_mut() + .par_col_iter_mut() + .zip_eq(lhs.par_col_partition(nthreads)) + .zip_eq(rhs.par_partition(nthreads)), + |((dst, lhs), rhs)| { + make_guard!(K); + let K = lhs.ncols().bind(K); + let lhs = lhs.as_col_shape(K); + let rhs = rhs.as_row_shape(K); + matvec( + dst, + Accum::Replace, + lhs, + conj_lhs, + rhs, + conj_rhs, + alpha, + Par::Seq, + ); + }, + ); + matvec( + dst.rb_mut(), + beta, + lhs.subcols(z, Z), + conj_lhs, + rhs.subrows(z, Z), + conj_rhs, + &zero(), + Par::Seq, + ); + for j in 0..nthreads { + zip!(dst.rb_mut(), tmp.rb().col(j)) + .for_each(|unzip!(dst, src)| *dst = &*dst + &*src) + } + }, + } + }); } } mod rank_update { @@ -1379,7 +1375,7 @@ fn matmul_imp<'M, 'N, 'K, T: ComplexField>( ); return; } else { - let par = match par { + let _ = match par { #[cfg(feature = "threadteam")] Par::Team(team) => { // the team path runs both kernels below on the team. On @@ -1463,24 +1459,9 @@ fn matmul_imp<'M, 'N, 'K, T: ComplexField>( rhs.col_stride(), conj_rhs == Conj::Yes, &raw const alpha as *const (), - // The kernel below (`private_gemm_x86`) threads whenever it is - // handed more than one thread and has no work floor of its own, - // so the decision is faer's: the serial time of the product - // (its MACs at the dtype's rate, `utils::cost::ns_per_mac`) - // against `PAYS` times the cost of one section of this shape - // (`ParShape::dispatch_cost`). One section: the kernel cuts - // the product into one grid of tiles. - if par.shape().pays_full( - $crate::utils::cost::serial_time( - (*M).saturating_mul(*N).saturating_mul(*K), - $crate::utils::cost::gemm_ns_per_mac::(), - ), - 1, - ) { - par.degree() - } else { - 1 - }, + // one thread: a team returned above (`gemm_team`), so `par` + // is sequential here + 1, // a team is the only parallel route: the `gemm` crate's own // parallelism is never handed a thread (the team runs the // kernel on its members instead, or the kernel stays serial) diff --git a/faer/src/linalg/qr/col_pivoting/factor.rs b/faer/src/linalg/qr/col_pivoting/factor.rs index 4808e5c2..d868bdbd 100644 --- a/faer/src/linalg/qr/col_pivoting/factor.rs +++ b/faer/src/linalg/qr/col_pivoting/factor.rs @@ -2,6 +2,7 @@ pub use super::super::no_pivoting::factor::recommended_block_size; use crate::assert; use crate::internal_prelude::*; use crate::utils::cost::{ParThreshold, rate}; +use crate::utils::thread::with_width; use crate::perm::swap_cols_idx; use linalg::householder::{self, HouseholderInfo}; use pulp::Simd; @@ -128,6 +129,10 @@ fn qr_in_place_unblocked<'out, I: Index, T: ComplexField>( let mut A = A; let mut H = H; let mut par = par; + // `Never` keeps the whole call sequential (the steps decide their own width otherwise) + if matches!(params.par_threshold, ParThreshold::Never) { + par = Par::Seq; + } assert!(size == Ord::min(m, n)); for j in 0..n { col_perm[j] = I::truncate(j); @@ -162,6 +167,9 @@ fn qr_in_place_unblocked<'out, I: Index, T: ComplexField>( } best = best * scale_bwd; let mut best_threshold = best * threshold; + // the width the previous step's decision chose: it also covers the product that the + // next step opens with when the update was not delayed + let mut width = par.degree(); 'unscale: { for k in 0..size { let mut new_best = zero::(); @@ -180,14 +188,16 @@ fn qr_in_place_unblocked<'out, I: Index, T: ComplexField>( let (_, _, A10, mut A11) = A.rb_mut().split_at_mut(k, k); let dot = dot.rb().get(k..); let A10 = A10.rb().col(k - 1); - linalg::matmul::matmul( - A11.rb_mut(), - Accum::Add, - A10, - dot, - one(), - par.reborrow(), - ); + with_width(par.reborrow(), width, |par| { + linalg::matmul::matmul( + A11.rb_mut(), + Accum::Add, + A10, + dot, + one(), + par, + ); + }); best = zero(); for j in k..n { let val = A11.rb().col(j - k).norm_l2(); @@ -247,67 +257,69 @@ fn qr_in_place_unblocked<'out, I: Index, T: ComplexField>( // One section per step: the fused update of the trailing block and the // norms (or, without the delayed update, the product, which the same // decision covers); the unit is an element of the trailing block. - if !params.par_threshold.region_pays::( - par.shape(), + width = params.par_threshold.region_width::( + &par, (m - k - 1) * (n - k - 1), rate::qr_col_pivot(), 1, - ) { - par = Par::Seq; - } + ); if delayed_update { let A10 = A10.col(k - 1).get(1..); - match par { - Par::Seq => { - update_mat_and_dot_simd( - norm.rb_mut(), - dot.rb_mut(), - B01.rb_mut(), - B11.rb_mut().try_as_col_major_mut().unwrap(), - A10.try_as_col_major().unwrap(), - B10.rb().try_as_col_major().unwrap(), - l.copy(), - tau_inv.copy(), - simd_align(k + 1), - ); - }, - #[cfg(feature = "rayon")] - _ => { - let nthreads = core::num::NonZeroUsize::new(par.degree()).unwrap(); - let nthreads = nthreads.get(); - use rayon::prelude::*; - crate::utils::thread::for_each(par.reborrow(), - nthreads, - norm.par_partition_mut(nthreads) - .zip(dot.par_partition_mut(nthreads)) - .zip(B01.par_partition_mut(nthreads)) - .zip(B11.par_col_partition_mut(nthreads)), - |(((norm, dot), B01), B11)| { - update_mat_and_dot_simd( - norm, - dot, - B01, - B11.try_as_col_major_mut().unwrap(), - A10.try_as_col_major().unwrap(), - B10.rb().try_as_col_major().unwrap(), - l.copy(), - tau_inv.copy(), - simd_align(k + 1), - ); - }, - ); - }, - } + with_width(par.reborrow(), width, |par| { + match par { + Par::Seq => { + update_mat_and_dot_simd( + norm.rb_mut(), + dot.rb_mut(), + B01.rb_mut(), + B11.rb_mut().try_as_col_major_mut().unwrap(), + A10.try_as_col_major().unwrap(), + B10.rb().try_as_col_major().unwrap(), + l.copy(), + tau_inv.copy(), + simd_align(k + 1), + ); + }, + #[cfg(feature = "rayon")] + _ => { + let nthreads = core::num::NonZeroUsize::new(par.degree()).unwrap(); + let nthreads = nthreads.get(); + use rayon::prelude::*; + crate::utils::thread::for_each(par, + nthreads, + norm.par_partition_mut(nthreads) + .zip(dot.par_partition_mut(nthreads)) + .zip(B01.par_partition_mut(nthreads)) + .zip(B11.par_col_partition_mut(nthreads)), + |(((norm, dot), B01), B11)| { + update_mat_and_dot_simd( + norm, + dot, + B01, + B11.try_as_col_major_mut().unwrap(), + A10.try_as_col_major().unwrap(), + B10.rb().try_as_col_major().unwrap(), + l.copy(), + tau_inv.copy(), + simd_align(k + 1), + ); + }, + ); + }, + } + }); } else { dot.copy_from(B01.rb()); - linalg::matmul::matmul( - dot.rb_mut(), - Accum::Add, - B10.rb().adjoint(), - B11.rb(), - one(), - par.reborrow(), - ); + with_width(par.reborrow(), width, |par| { + linalg::matmul::matmul( + dot.rb_mut(), + Accum::Add, + B10.rb().adjoint(), + B11.rb(), + one(), + par, + ); + }); zip!(B01.rb_mut(), dot.rb_mut(), norm.rb_mut()).for_each( |unzip!(a, dot, norm)| { *dot = -dot.mul_real(tau_inv); @@ -377,6 +389,10 @@ pub fn qr_in_place<'out, I: Index, T: ComplexField>( // copy's, so the factor, the `tau` matrix and the permutations are the // column-major ones. let mut A = A; + // `Never` keeps the whole call sequential, the compact-factor upgrade included + if matches!(params.config.par_threshold, ParThreshold::Never) { + par = Par::Seq; + } if A.row_stride().unsigned_abs() != 1 && A.nrows() > 1 && A.ncols() > 1 diff --git a/faer/src/linalg/svd/bidiag.rs b/faer/src/linalg/svd/bidiag.rs index 5bcf6081..5485a3bb 100644 --- a/faer/src/linalg/svd/bidiag.rs +++ b/faer/src/linalg/svd/bidiag.rs @@ -2,6 +2,7 @@ use crate::assert; use crate::internal_prelude::*; use linalg::householder::*; use crate::utils::cost::{ParThreshold, rate}; +use crate::utils::thread::with_width; use linalg::matmul::{dot, matmul}; /// width of one blocked panel: `BLOCK_SIZE` columns and rows are reduced with level-2 work @@ -99,6 +100,10 @@ pub fn bidiag_in_place( let mut Hl = H_left; let mut Hr = H_right; let mut par = par; + // `Never` keeps the whole call sequential, the compact-factor upgrades included + if matches!(params.par_threshold, ParThreshold::Never) { + par = Par::Seq; + } // blocked phase: whole panels while the remaining dimension clears the floor. its regions are // two vector products per step and one update per panel, and they run fastest when the caller // holds a team around the call (see the crate docs): spinning members that keep their column @@ -196,27 +201,23 @@ fn bidiag_blocked_phase( let mut w = w.as_mat_mut(); let mut p = 0; while size - p > BLOCK_FLOOR { - // the decision is made on the full remaining area, never latched by a panel - let par = if par_threshold.region_pays::( - par.shape(), - (m - p) * (n - p), - rate::bidiag(), - 2, - ) { - par.reborrow() - } else { - Par::Seq - }; - bidiag_blocked_panel( - A.rb_mut().get_mut(p.., p..), - Hl.rb_mut().get_mut(p..p + BLOCK_SIZE), - Hr.rb_mut().get_mut(p..p + BLOCK_SIZE), - X.rb_mut().get_mut(p.., ..), - Y.rb_mut().get_mut(p.., ..), - de.rb_mut(), - w.rb_mut().col_mut(0).get_mut(p..), - par, - ); + // the decision is made on the full remaining area, never latched by a panel. It is the + // floor of the panel only: the panel's own products choose their widths, so a paying + // panel keeps the whole `par` rather than capping them at this estimate's width + let width = par_threshold.region_width::(&par, (m - p) * (n - p), rate::bidiag(), 2); + let width = if width > 1 { par.degree() } else { 1 }; + with_width(par.reborrow(), width, |par| { + bidiag_blocked_panel( + A.rb_mut().get_mut(p.., p..), + Hl.rb_mut().get_mut(p..p + BLOCK_SIZE), + Hr.rb_mut().get_mut(p..p + BLOCK_SIZE), + X.rb_mut().get_mut(p.., ..), + Y.rb_mut().get_mut(p.., ..), + de.rb_mut(), + w.rb_mut().col_mut(0).get_mut(p..), + par, + ); + }); p += BLOCK_SIZE; } p @@ -512,6 +513,10 @@ fn bidiag_unblocked( let mut Hl = Hl; let mut Hr = Hr; let mut par = par; + // `Never` keeps the whole call sequential (the steps decide their own width otherwise) + if matches!(par_threshold, ParThreshold::Never) { + par = Par::Seq; + } let mut w = w0.rb_mut(); let mut wp = w1.rb_mut(); for k in 0..size { @@ -548,62 +553,60 @@ fn bidiag_unblocked( make_householder_in_place(a11, A21.rb_mut()); let ref tl_inv = tl.recip(); Hl[k] = tl.to_cplx(); - if !par_threshold.region_pays::( - par.shape(), - (m - k - 1) * (n - k - 1), - rate::bidiag(), - 2, - ) { - par = Par::Seq; - } + // the step's two products run at the widest width that pays, decided per step + let width = par_threshold.region_width::(&par, (m - k - 1) * (n - k - 1), rate::bidiag(), 2); if k > 0 { let k1 = k - 1; let up = A20.rb().col(k1); - match par.reborrow() { - Par::Seq => bidiag_fused_op( - A22.rb_mut(), - A21.rb(), - up.rb(), - z2.rb(), - y2.rb_mut(), - vp, - simd_align(row_offset + k + 1), - ), - #[cfg(feature = "rayon")] - _ => { - use rayon::prelude::*; - let nthreads = - core::num::NonZeroUsize::new(par.degree()).unwrap().get(); - crate::utils::thread::for_each( - par.reborrow(), - nthreads, - A22.rb_mut() - .par_col_partition_mut(nthreads) - .zip_eq(y2.rb_mut().par_partition_mut(nthreads)) - .zip_eq(vp.par_partition(nthreads)), - |((A22, y2), vp)| { - bidiag_fused_op( - A22, - A21.rb(), - up.rb(), - z2.rb(), - y2, - vp.rb(), - simd_align(row_offset + k + 1), - ); - }, - ); - }, - } + with_width(par.reborrow(), width, |mut par| { + match par.reborrow() { + Par::Seq => bidiag_fused_op( + A22.rb_mut(), + A21.rb(), + up.rb(), + z2.rb(), + y2.rb_mut(), + vp, + simd_align(row_offset + k + 1), + ), + #[cfg(feature = "rayon")] + _ => { + use rayon::prelude::*; + let nthreads = + core::num::NonZeroUsize::new(par.degree()).unwrap().get(); + crate::utils::thread::for_each( + par.reborrow(), + nthreads, + A22.rb_mut() + .par_col_partition_mut(nthreads) + .zip_eq(y2.rb_mut().par_partition_mut(nthreads)) + .zip_eq(vp.par_partition(nthreads)), + |((A22, y2), vp)| { + bidiag_fused_op( + A22, + A21.rb(), + up.rb(), + z2.rb(), + y2, + vp.rb(), + simd_align(row_offset + k + 1), + ); + }, + ); + }, + } + }); } else { - matmul( - y2.rb_mut(), - Accum::Replace, - A21.rb().adjoint(), - A22.rb(), - one(), - par.reborrow(), - ); + with_width(par.reborrow(), width, |par| { + matmul( + y2.rb_mut(), + Accum::Replace, + A21.rb().adjoint(), + A22.rb(), + one(), + par, + ); + }); } z!(y2.rb_mut(), wk.rb_mut()).for_each(|uz!(y, a)| { *y = (&*y + &*a).mul_real(tl_inv); @@ -623,14 +626,16 @@ fn bidiag_unblocked( }; z!(wk.rb_mut()).for_each(|uz!(a)| *a = a.mul_real(norm_inv)); } - matmul( - z2.rb_mut(), - Accum::Replace, - A22.rb(), - wk.rb().adjoint(), - one(), - par.reborrow(), - ); + with_width(par.reborrow(), width, |par| { + matmul( + z2.rb_mut(), + Accum::Replace, + A22.rb(), + wk.rb().adjoint(), + one(), + par, + ); + }); if k + 1 == size { A12.copy_from(wk.rb()); break; diff --git a/faer/src/linalg/triangular_solve.rs b/faer/src/linalg/triangular_solve.rs index fbeda8ec..bcc7c53e 100644 --- a/faer/src/linalg/triangular_solve.rs +++ b/faer/src/linalg/triangular_solve.rs @@ -467,29 +467,30 @@ fn solve_unit_lower_triangular_in_place_imp<'N, 'K, T: ComplexField>( // rule for a sequential run. #[cfg(feature = "rayon")] { - let jobs = Ord::min(par.degree(), k.div_ceil(64)); - if jobs >= 2 - && crate::utils::cost::par_pays::( - &par, - n.saturating_mul(n).saturating_mul(k) / 2, - crate::utils::cost::rate::gemm_mac(), - 1, - ) { - crate::utils::thread::for_each_with_par( - par, - jobs, - rhs.as_dyn_cols_mut().par_col_partition_mut(jobs), - |chunk, par| { - make_guard!(CHUNK); - let width = chunk.ncols().bind(CHUNK); - solve_unit_lower_triangular_in_place_imp( - tril, - conj_lhs, - chunk.as_col_shape_mut(width), - par, - ) - }, - ); + // the split: the widest number of jobs, at most `ceil(k / 64)` (the inherited cap), that + // pays for its section; the region: the widest width that pays, shared out among the jobs + let units = n.saturating_mul(n).saturating_mul(k) / 2; + let rate = crate::utils::cost::rate::gemm_mac(); + let jobs = crate::utils::cost::par_width_upto::(&par, units, rate, 1, k.div_ceil(64)); + let width = crate::utils::cost::par_width::(&par, units, rate, 1); + if jobs >= 2 { + crate::utils::thread::with_width(par, width, |par| { + crate::utils::thread::for_each_with_par( + par, + jobs, + rhs.as_dyn_cols_mut().par_col_partition_mut(jobs), + |chunk, par| { + make_guard!(CHUNK); + let width = chunk.ncols().bind(CHUNK); + solve_unit_lower_triangular_in_place_imp( + tril, + conj_lhs, + chunk.as_col_shape_mut(width), + par, + ) + }, + ) + }); return; } } @@ -589,29 +590,30 @@ fn solve_lower_triangular_in_place_imp<'N, 'K, T: ComplexField>( if k > 64 && n <= 128 { #[cfg(feature = "rayon")] { - let jobs = Ord::min(par.degree(), k.div_ceil(64)); - if jobs >= 2 - && crate::utils::cost::par_pays::( - &par, - n.saturating_mul(n).saturating_mul(k) / 2, - crate::utils::cost::rate::gemm_mac(), - 1, - ) { - crate::utils::thread::for_each_with_par( - par, - jobs, - rhs.as_dyn_cols_mut().par_col_partition_mut(jobs), - |chunk, par| { - make_guard!(CHUNK); - let width = chunk.ncols().bind(CHUNK); - solve_lower_triangular_in_place_imp( - tril, - conj_lhs, - chunk.as_col_shape_mut(width), - par, - ) - }, - ); + // the split: the widest number of jobs, at most `ceil(k / 64)` (the inherited cap), that + // pays for its section; the region: the widest width that pays, shared out among the jobs + let units = n.saturating_mul(n).saturating_mul(k) / 2; + let rate = crate::utils::cost::rate::gemm_mac(); + let jobs = crate::utils::cost::par_width_upto::(&par, units, rate, 1, k.div_ceil(64)); + let width = crate::utils::cost::par_width::(&par, units, rate, 1); + if jobs >= 2 { + crate::utils::thread::with_width(par, width, |par| { + crate::utils::thread::for_each_with_par( + par, + jobs, + rhs.as_dyn_cols_mut().par_col_partition_mut(jobs), + |chunk, par| { + make_guard!(CHUNK); + let width = chunk.ncols().bind(CHUNK); + solve_lower_triangular_in_place_imp( + tril, + conj_lhs, + chunk.as_col_shape_mut(width), + par, + ) + }, + ) + }); return; } } diff --git a/faer/src/sparse/linalg/matmul.rs b/faer/src/sparse/linalg/matmul.rs index d271bf40..98d1bfc3 100644 --- a/faer/src/sparse/linalg/matmul.rs +++ b/faer/src/sparse/linalg/matmul.rs @@ -112,15 +112,15 @@ pub fn sparse_sparse_matmul_numeric< unsafe impl Sync for SyncWrapper {} unsafe impl Send for SyncWrapper {} let c_values = SyncWrapper(&*UnsafeCell::from_mut(c_values)); - // one section over the column ranges; the work is the symbolic flop count, which is already - // known. The workspace keeps one column per thread of `par` either way. - let parallel = crate::utils::cost::par_pays::( + // one section over the column ranges, at the widest width that pays; the work is the + // symbolic flop count, which is already known. The workspace keeps one column per thread of + // `par` either way. + let nthreads = crate::utils::cost::par_width::( &par, total_flop_count as usize, crate::utils::cost::rate::sparse_matmul(), 1, ); - let nthreads = if parallel { par.degree() } else { 1 }; let job = &|tid: usize| { assert!(tid < nthreads); fn partition_fn( @@ -170,26 +170,22 @@ pub fn sparse_sparse_matmul_numeric< } } }; - match par { - Par::Seq => { - job(0); - }, - #[cfg(feature = "rayon")] - _ if !parallel => { - job(0); - }, - #[cfg(feature = "rayon")] - _ => { - let nthreads = core::num::NonZeroUsize::new(par.degree()).unwrap(); - use rayon::prelude::*; - crate::utils::thread::for_each( - par, - nthreads.get(), - (0..nthreads.get()).into_par_iter(), - |tid| job(tid), - ); - }, + if nthreads < 2 { + job(0); + return; } + #[cfg(feature = "rayon")] + crate::utils::thread::with_width(par, nthreads, |par| { + use rayon::prelude::*; + crate::utils::thread::for_each( + par, + nthreads, + (0..nthreads).into_par_iter(), + |tid| job(tid), + ); + }); + #[cfg(not(feature = "rayon"))] + unreachable!("no parallel route in this build"); } /// performs a numeric matrix multiplication of a sparse matrix `lhs` by a diff --git a/faer/src/utils/cost.rs b/faer/src/utils/cost.rs index 084f7ed4..de442b11 100644 --- a/faer/src/utils/cost.rs +++ b/faer/src/utils/cost.rs @@ -207,43 +207,31 @@ impl From for usize { } impl ParThreshold { - /// Whether a region of `units` of work in the unit of `rate` that opens `sections` - /// sections runs in parallel on `shape`, at the full degree of the shape. + /// The width a region of `units` of work in the unit of `rate` that opens `sections` + /// sections runs at on `par` (1: sequential); the caller runs it on a sub-team of that width + /// ([`with_width`](crate::utils::thread::with_width)). /// - /// [`Derived`](Self::Derived) is [`ParShape::pays_full`] of `units * rate`; an explicit - /// [`Fixed`](Self::Fixed) floor keeps the previous rule, `units >= floor`. + /// [`Derived`](Self::Derived) is [`par_width`] of `units * rate`, the widest width that pays; + /// an explicit [`Fixed`](Self::Fixed) floor keeps its rule, the whole degree when + /// `units >= floor`; [`Never`](Self::Never) is 1. #[cfg_attr(feature = "par-trace", track_caller)] - pub(crate) fn region_pays( + pub(crate) fn region_width( self, - shape: ParShape, + par: &crate::Par, units: usize, rate: Rate, sections: u32, - ) -> bool { - // a sequential shape, or a token of one participant (a `Team` of size 1), never pays: - // skip the estimate (a QR 128 asks hundreds of times) - if shape.degree() < 2 && !matches!(self, Self::Fixed(_)) { - return false; - } - let r = self.region_pays_impl::(shape, units, rate, sections); - #[cfg(feature = "par-trace")] - trace_decision(shape, shape.degree(), units, rate.of::(), sections, r, core::panic::Location::caller()); - r - } - - fn region_pays_impl( - self, - shape: ParShape, - units: usize, - rate: Rate, - sections: u32, - ) -> bool { + ) -> usize { match self { - Self::Never => false, - Self::Fixed(floor) => units >= floor, - Self::Derived => { - shape.pays_full(serial_time(units, rate.of::()), sections) + Self::Never => 1, + Self::Fixed(floor) => { + if units >= floor { + par.degree() + } else { + 1 + } }, + Self::Derived => par_width::(par, units, rate, sections), } } } @@ -520,8 +508,10 @@ pub mod trace { pub sections: u32, /// estimated serial time in ns pub estimate_ns: f64, - /// dispatch cost of one section at the shape's width, in ns + /// dispatch cost of one section at `width`, in ns pub dispatch_ns: f64, + /// the width the gate chose, or the width it asked about (1: sequential) + pub width: usize, /// the decision pub parallel: bool, } @@ -554,6 +544,7 @@ pub mod trace { sections: u32, estimate_ns: f64, dispatch_ns: f64, + width: usize, parallel: bool, ) { if !ON.load(Ordering::Relaxed) { @@ -565,7 +556,7 @@ pub mod trace { } else { label }; - LOG.lock().unwrap().push(Record { site, dims, units, sections, estimate_ns, dispatch_ns, parallel }); + LOG.lock().unwrap().push(Record { site, dims, units, sections, estimate_ns, dispatch_ns, width, parallel }); } } @@ -616,24 +607,72 @@ pub(crate) fn par_can_pay_time_team(team: &threadteam::TeamExec<'_>, estimate: D && pays_at(estimate, sections, team.dispatch_costs().min_cost(threadteam::DispatchKind::Team)) } -/// Whether a region of `units` of work at `rate` that opens `sections` sections pays on `par`, -/// at the full degree of `par` ([`ParShape::pays_full`]). +/// The widest width in `2..=degree` at which a region of serial time `estimate` that opens +/// `sections` sections pays on `par`, or 1 (sequential) when none does. +/// +/// Every width is considered: the search runs from the full width down and stops at the first +/// that pays, which is the widest. A narrower width that does not pay says nothing about a wider +/// one, because a measured row can fall with the width (`aarch64-apple`'s team row is 3.0, 2.3, +/// 1.4 us at widths 2, 4, 6). +pub(crate) fn par_width_time(par: &crate::Par, estimate: Duration, sections: u32) -> usize { + par_width_time_upto(par, estimate, sections, usize::MAX) +} + +/// [`par_width_time`] among the widths up to `max` (a region that cannot use more participants, +/// such as a split into at most `max` tasks). +pub(crate) fn par_width_time_upto(par: &crate::Par, estimate: Duration, sections: u32, max: usize) -> usize { + match par { + #[cfg(feature = "threadteam")] + crate::Par::Team(team) => par_width_team(team, estimate, sections, max), + _ => { + let _ = (estimate, sections, max); + 1 + }, + } +} + +/// [`par_width_time_upto`] for a team token. +#[cfg(feature = "threadteam")] +pub(crate) fn par_width_team(team: &threadteam::TeamExec<'_>, estimate: Duration, sections: u32, max: usize) -> usize { + // the cheap rejection first: a region that does not pay at the cheapest width of the table + // pays at none (a QR 128 asks hundreds of times) + if max < 2 || !par_can_pay_time_team(team, estimate, sections) { + return 1; + } + (2..=team.size().get().min(max)) + .rev() + .find(|&w| pays_at(estimate, sections, team.dispatch_cost(NonZeroUsize::new(w).unwrap()))) + .unwrap_or(1) +} + +/// [`par_width_time`] of `units` of work at `rate`: the width the region runs at (1: +/// sequential). #[cfg_attr(feature = "par-trace", track_caller)] -pub(crate) fn par_pays( +pub(crate) fn par_width( par: &crate::Par, units: usize, rate: Rate, sections: u32, -) -> bool { - // a sequential token, or one of a single participant (a `Team` of size 1: it behaves exactly - // like `Par::Seq`), never pays: skip the estimate (a QR 128 asks 500 times) - if par.degree() < 2 { - return false; +) -> usize { + par_width_upto::(par, units, rate, sections, usize::MAX) +} + +/// [`par_width`] among the widths up to `max`. +#[cfg_attr(feature = "par-trace", track_caller)] +pub(crate) fn par_width_upto( + par: &crate::Par, + units: usize, + rate: Rate, + sections: u32, + max: usize, +) -> usize { + if par.degree() < 2 || max < 2 { + return 1; } - let r = par_pays_time(par, serial_time(units, rate.of::()), sections, par.degree()); + let w = par_width_time_upto(par, serial_time(units, rate.of::()), sections, max); #[cfg(feature = "par-trace")] - trace_decision(par.shape(), par.shape().degree(), units, rate.of::(), sections, r, core::panic::Location::caller()); - r + trace_decision(par.shape(), w, units, rate.of::(), sections, w > 1, core::panic::Location::caller()); + w } #[cfg(feature = "par-trace")] @@ -649,7 +688,7 @@ fn trace_decision( let d = NonZeroUsize::new(width) .map(|w| shape.dispatch_cost(w).as_secs_f64() * 1e9) .unwrap_or(f64::INFINITY); - trace::record(loc, units, sections, units as f64 * ns, d, r); + trace::record(loc, units, sections, units as f64 * ns, d, width, r); } /// Whether a join of two halves pays: one section of width 2 (the two halves are the @@ -847,7 +886,7 @@ mod tests { sections: u32, ) -> usize { let pays = |u: usize| { - ParThreshold::Derived.region_pays::(shape, u, r, sections) + shape.pays_full(serial_time(u, r.of::()), sections) }; let mut hi = 1usize; while !pays(hi) { @@ -1020,12 +1059,14 @@ mod tests { #[test] fn fixed_and_never_keep_their_meaning() { - let s = ParShape::team(6); let r = rate::lu_full_pivot(); - assert!(!ParThreshold::Never.region_pays::(s, usize::MAX, r, 1)); - assert!(ParThreshold::Fixed(100).region_pays::(ParShape::Seq, 100, r, 1)); - assert!(!ParThreshold::Fixed(100).region_pays::(s, 99, r, 1)); - assert!(ParThreshold::Fixed(100).region_pays::(s, 100, r, 1)); + // a fixed floor takes the whole degree at the floor, nothing below it; never is serial + assert_eq!(ParThreshold::Fixed(100).region_width::(&crate::Par::Seq, 100, r, 1), 1); + crate::team::with_team(6, |par| { + assert_eq!(ParThreshold::Never.region_width::(&par, usize::MAX, r, 1), 1); + assert_eq!(ParThreshold::Fixed(100).region_width::(&par, 99, r, 1), 1); + assert_eq!(ParThreshold::Fixed(100).region_width::(&par, 100, r, 1), 6); + }); // the integer form of the C API assert_eq!(ParThreshold::from(0), ParThreshold::Derived); assert_eq!(ParThreshold::from(usize::MAX), ParThreshold::Never); @@ -1035,19 +1076,11 @@ mod tests { #[test] fn derived_gates_never_pay_on_a_sequential_shape() { + let seq = crate::Par::Seq; for (_, r, sections) in gates(rates()) { - assert!(!ParThreshold::Derived.region_pays::( - ParShape::Seq, - usize::MAX / 2, - r, - sections - )); - assert!(!ParThreshold::Derived.region_pays::( - ParShape::Seq, - usize::MAX / 2, - r, - sections - )); + assert_eq!(ParThreshold::Derived.region_width::(&seq, usize::MAX / 2, r, sections), 1); + assert_eq!(ParThreshold::Derived.region_width::(&seq, usize::MAX / 2, r, sections), 1); + assert!(!ParShape::Seq.pays_full(serial_time(usize::MAX / 2, r.of::()), sections)); } } @@ -1073,7 +1106,7 @@ mod tests { } // representative sizes at width 6 on the M5 table (D = 1.4 us, 4.2 us per section) let pays = |u: usize, r: Rate, f: u32| { - ParThreshold::Derived.region_pays::(team, u, r, f) + team.pays_full(serial_time(u, r.of::()), f) }; // full pivoting LU: a 100 x 100 trailing block stays serial, 200 x 200 is parallel assert!(!pays(100 * 100, M5.lu_full_pivot, 1)); @@ -1091,15 +1124,17 @@ mod tests { } #[test] - fn par_pays_reads_the_degree_of_the_token() { + fn par_width_reads_the_degree_of_the_token() { let seq = crate::Par::Seq; - assert!(!par_pays::(&seq, usize::MAX / 2, rate::matvec(), 1)); + assert_eq!(par_width::(&seq, usize::MAX / 2, rate::matvec(), 1), 1); crate::team::with_team(6, |par| { - assert!(par_pays::(&par, 10_000_000, rate::matvec(), 1)); - assert!(!par_pays::(&par, 1_000, rate::matvec(), 1)); + assert_eq!(par_width::(&par, 10_000_000, rate::matvec(), 1), 6); + assert_eq!(par_width::(&par, 1_000, rate::matvec(), 1), 1); + // a cap below the degree bounds the search + assert_eq!(par_width_upto::(&par, 10_000_000, rate::matvec(), 1, 3), 3); }); crate::team::with_team(1, |par| { - assert!(!par_pays::(&par, 10_000_000, rate::matvec(), 1)); + assert_eq!(par_width::(&par, 10_000_000, rate::matvec(), 1), 1); }); } @@ -1416,8 +1451,8 @@ mod tests { assert!(!pays(no)); } // a sparse product of 1e6 multiply-adds is above the floor, 1e3 below - assert!(ParThreshold::Derived.region_pays::(team, 1_000_000, M5.sparse_matmul, 1)); - assert!(!ParThreshold::Derived.region_pays::(team, 1_000, M5.sparse_matmul, 1)); + assert!(team.pays_full(serial_time(1_000_000, M5.sparse_matmul.of::()), 1)); + assert!(!team.pays_full(serial_time(1_000, M5.sparse_matmul.of::()), 1)); } } @@ -1428,10 +1463,11 @@ pub fn trace_pub( sections: u32, estimate: Duration, dispatch: Duration, + width: usize, r: bool, loc: &'static core::panic::Location<'static>, ) { - trace::record(loc, units, sections, estimate.as_secs_f64() * 1e9, dispatch.as_secs_f64() * 1e9, r); + trace::record(loc, units, sections, estimate.as_secs_f64() * 1e9, dispatch.as_secs_f64() * 1e9, width, r); } #[cfg(all(test, feature = "sparse-linalg"))] diff --git a/faer/src/utils/mod.rs b/faer/src/utils/mod.rs index 71cf4baa..9ab7c8b6 100644 --- a/faer/src/utils/mod.rs +++ b/faer/src/utils/mod.rs @@ -68,6 +68,45 @@ pub mod thread { parallelism, ) } + /// Runs `f` on a part of `par` of `width` participants: `par` itself when `width` is its + /// degree or more, `Par::Seq` below two, otherwise a sub-team of `width` (a split into one + /// child, which is led by the caller and dispatches nothing; the members it leaves out stay + /// idle). A region decides its width with [`cost::par_width`](crate::utils::cost) and runs + /// here, so only that region is narrowed and the caller keeps its whole `par`. + pub fn with_width(par: Par, width: usize, f: impl FnOnce(Par) -> R) -> R { + if width >= par.degree() { + return f(par); + } + if width < 2 { + return f(Par::Seq); + } + match par { + Par::Seq => f(Par::Seq), + #[cfg(feature = "threadteam")] + Par::Team(team) => { + super::team::note_width(width); + struct Once(core::cell::UnsafeCell<(Option, Option)>); + // SAFETY: a split into the single child `[width]` runs its callback once, on the + // calling thread (child 0 is led by the caller), so the cell is never shared. + unsafe impl Sync for Once {} + impl Once { + fn get(&self) -> *mut (Option, Option) { + self.0.get() + } + } + let once = Once(core::cell::UnsafeCell::new((Some(f), None))); + team.split(&[width], |_, child| { + // SAFETY: see `Once`. + let (f, out) = unsafe { &mut *once.get() }; + *out = Some((f.take().expect("one child"))(Par::Team(child))); + }) + .expect("2 <= width < size is a valid split"); + once.0.into_inner().1.expect("the child ran") + }, + #[cfg(not(feature = "threadteam"))] + Par::__Lifetime(x, _) => match x {}, + } + } /// unsafe [`Send`] and [`Sync`] pointer type pub struct Ptr(pub *mut T); unsafe impl Send for Ptr {} diff --git a/faer/src/utils/team.rs b/faer/src/utils/team.rs index 5b74664b..27751dc0 100644 --- a/faer/src/utils/team.rs +++ b/faer/src/utils/team.rs @@ -120,13 +120,17 @@ pub mod stats { pub(super) static RUN: [AtomicUsize; W] = [Z; W]; /// splits and forks by number of parts (index 1 = ran inline on the caller) pub(super) static SPLIT: [AtomicUsize; W] = [Z; W]; + /// bins of [`WIDTH`]: one per width up to 256 (the last bin collects anything wider) + const WW: usize = 257; + /// regions narrowed to a sub-team, by width (`thread::with_width` and the x86 kernel's fork) + pub(super) static WIDTH: [AtomicUsize; WW] = [Z; WW]; pub(super) static GEMM_BELOW_FLOOR: AtomicUsize = AtomicUsize::new(0); pub(super) static GEMM_ONE_TILE: AtomicUsize = AtomicUsize::new(0); pub(super) static GEMM_TILED: AtomicUsize = AtomicUsize::new(0); pub(super) static GEMM_TEAM_KERNEL: AtomicUsize = AtomicUsize::new(0); - pub(super) fn hist(h: &[AtomicUsize; W], n: usize) { + pub(super) fn hist(h: &[AtomicUsize; N], n: usize) { if ON.load(Relaxed) { - h[n.min(W - 1)].fetch_add(1, Relaxed); + h[n.min(N - 1)].fetch_add(1, Relaxed); } } pub(super) fn bump(c: &AtomicUsize) { @@ -141,6 +145,8 @@ pub mod stats { pub run: Vec, /// splits and forks, index = parts pub split: Vec, + /// regions narrowed to a sub-team, index = its width (up to 256) + pub width: Vec, /// GEMM calls decided sequential by the work floor pub gemm_below_floor: usize, /// GEMM calls whose tile grid had a single tile @@ -157,10 +163,13 @@ pub mod stats { } /// Returns and resets the counters. pub fn take() -> Stats { - let t = |h: &[AtomicUsize; W]| h.iter().map(|c| c.swap(0, Relaxed)).collect(); + fn t(h: &[AtomicUsize; N]) -> Vec { + h.iter().map(|c| c.swap(0, Relaxed)).collect() + } Stats { run: t(&RUN), split: t(&SPLIT), + width: t(&WIDTH), gemm_below_floor: GEMM_BELOW_FLOOR.swap(0, Relaxed), gemm_one_tile: GEMM_ONE_TILE.swap(0, Relaxed), gemm_tiled: GEMM_TILED.swap(0, Relaxed), @@ -169,6 +178,11 @@ pub mod stats { } } +/// Records a region narrowed to a sub-team of `width` (probes only). +pub(crate) fn note_width(width: usize) { + stats::hist(&stats::WIDTH, width); +} + /// Records a split of the team into `parts` concurrent children (probes only). pub(crate) fn note_split(parts: usize) { stats::hist(&stats::SPLIT, parts); @@ -375,27 +389,15 @@ pub(crate) mod kernel { /// The width at which `team` runs a product of serial time `estimate` through the kernel, or /// `None` (run sequentially on the caller) when no width pays. The kernel forks twice per - /// call (packing, then compute), so the gate is `pays(estimate, 2, width)` at the width - /// `width_for` chooses. + /// call (packing, then compute): the widest width at which two sections pay. pub(crate) fn kernel_width(team: &TeamExec<'_>, estimate: Duration) -> Option { - // the cheap rejection first: a product that does not pay at the cheapest width of the table - // pays nowhere (whatever the shape of the table), so `width_for` need not search (a QR 128 - // asks for every product of its 59 block applies) - if !crate::utils::cost::par_can_pay_time_team(team, estimate, 2) { - stats::bump(&stats::GEMM_BELOW_FLOOR); - return None; - } - // the widest width that pays, on the cost the shapes read, searched from the full width - // down: a measured row can fall with the width (aarch64-apple's team row is 3.0, 2.3, - // 1.4 us at widths 2, 4, 6), so a narrow width that does not pay says nothing about a wide one - let size = team.size().get(); - let best = (2..=size).rev().find(|&w| crate::utils::cost::pays_at(estimate, 2, team.dispatch_cost(core::num::NonZeroUsize::new(w).unwrap()))); - match best { - Some(width) => Some(width), - None => { + // the widest width that pays, the rule of every region (`cost::par_width_time`) + match crate::utils::cost::par_width_team(team, estimate, 2, usize::MAX) { + 1 => { stats::bump(&stats::GEMM_BELOW_FLOOR); None }, + width => Some(width), } } @@ -405,6 +407,9 @@ pub(crate) mod kernel { stats::bump(&stats::GEMM_TEAM_KERNEL); let size = team.size().get(); let width = width.clamp(1, size); + if width < size { + super::note_width(width); + } if width == size { return f(&TeamFork { team: RefCell::new(team), width }); } @@ -579,18 +584,12 @@ pub(crate) unsafe fn gemm_tiled( gemm::Parallelism::None, ) }; - let nt = team.size().get(); // The floor is time, not a count: the serial time of the product (its MACs at the dtype's - // rate) against `PAYS` times the cost of one section of this team. A complex MAC is 4.5 - // real ones, so the floor of a complex product is that much lower in MACs. + // rate) against `PAYS` times the cost of one section, at the widest width that pays. A + // complex MAC is 4.5 real ones, so the floor of a complex product is that much lower in MACs. let estimate = crate::utils::cost::serial_time(m.saturating_mul(n).saturating_mul(k), ns_per_mac); - if nt <= 1 - || !crate::utils::cost::pays_at( - estimate, - 1, - team.dispatch_cost(core::num::NonZeroUsize::new(nt).unwrap()), - ) - { + let nt = crate::utils::cost::par_width_team(&team, estimate, 1, usize::MAX); + if nt <= 1 { stats::bump(&stats::GEMM_BELOW_FLOOR); return seq(dst, lhs, rhs, m, n); } @@ -599,6 +598,17 @@ pub(crate) unsafe fn gemm_tiled( stats::bump(&stats::GEMM_ONE_TILE); return seq(dst, lhs, rhs, m, n); } + // the section runs on `pm * pn` participants, which the grid can leave below `nt`: that width + // must pay as well (a row need not be monotone) + if pm * pn < nt + && !crate::utils::cost::pays_at(estimate, 1, team.dispatch_cost(core::num::NonZeroUsize::new(pm * pn).unwrap())) + { + stats::bump(&stats::GEMM_BELOW_FLOOR); + return seq(dst, lhs, rhs, m, n); + } + if pm * pn < team.size().get() { + note_width(pm * pn); + } stats::bump(&stats::GEMM_TILED); let (dst, lhs, rhs) = ( SendPtr(dst), @@ -630,6 +640,7 @@ mod tests { use crate::{c64, Accum, Mat, MatRef, Par}; use core::sync::atomic::{AtomicUsize, Ordering}; use threadteam::Pool; + use core::time::Duration as StdDuration; fn team_of(width: usize) -> threadteam::Team { Pool::builder() @@ -763,6 +774,296 @@ mod tests { .expect("a team token"); } + /// A team of `width` whose dispatch cost row is `anchors` (widths in increasing order). + fn with_row(width: usize, anchors: &[(usize, StdDuration)], f: impl FnOnce(Par<'_>) -> R) -> R { + use threadteam::{DispatchCosts, DispatchKind, WidthCosts}; + let _turn = team_turn(); + let costs = DispatchCosts::for_class(threadteam::HostClass::Generic) + .with(DispatchKind::Team, WidthCosts::from_anchors(anchors).expect("a valid row")); + let team = Pool::builder() + .threads(width) + .acquire_timeout(std::time::Duration::from_secs(30)) + .dispatch_costs(costs) + .build() + .unwrap() + .into_team(width) + .expect("an idle pool"); + let mut exec = team.exec(width); + exec.team(|team| f(Par::Team(team))).ok().expect("a team token") + } + + /// The width of a region is the widest that pays, every width considered: on a row that + /// rises past 4 the region stops at 4, on a row that falls it takes a width past one that + /// does not pay, and a cap below the degree searches inside the cap. + #[test] + fn par_width_takes_the_widest_width_that_pays() { + use crate::utils::cost::{par_width_time, par_width_time_upto}; + let us = StdDuration::from_micros; + // 1 us at 2 and 4, 100 us at 8: 6 us of work pays for one section at 4 (3 * 2 us) but not + // at 5 to 8 (interpolated 26.5 us and up) + with_row(8, &[(2, us(1)), (4, us(2)), (8, us(100))], |par| { + assert_eq!(par_width_time(&par, us(6), 1), 4); + // 5 us pays at 3 (interpolated 1.5 us) but not at 4 + assert_eq!(par_width_time(&par, us(5), 1), 3); + assert_eq!(par_width_time(&par, us(2), 1), 1); + assert_eq!(par_width_time(&par, us(2000), 1), 8); + }); + // a falling row: 9 us at 4, 1 us at 6 and 8; 3 us pays at 8 only, 9 us at 2 inside a cap of 5 + with_row(8, &[(2, us(3)), (4, us(9)), (6, us(1)), (8, us(1))], |par| { + assert_eq!(par_width_time(&par, us(3), 1), 8); + assert_eq!(par_width_time_upto(&par, us(3), 1, 5), 1); + assert_eq!(par_width_time_upto(&par, us(9), 1, 5), 2); + }); + assert_eq!(par_width_time(&Par::Seq, us(1_000_000), 1), 1); + } + + /// `with_width` hands the region a sub-team of exactly that width, the whole token at its + /// degree or more, and `Seq` below two, and the parent keeps its width afterwards. + #[test] + fn with_width_runs_on_a_sub_team_of_that_width() { + use crate::utils::thread::with_width; + with_par(8, |mut par| { + stats::enable(true); + let _ = stats::take(); + for (w, want) in [(3, 3), (8, 8), (9, 8), (1, 1), (0, 1)] { + let seen = AtomicUsize::new(0); + let degree = with_width(par.reborrow(), w, |par| { + let degree = par.degree(); + if let Par::Team(team) = par { + team.run(|_, _| { + seen.fetch_add(1, Ordering::Relaxed); + }); + } + degree + }); + assert_eq!(degree, want, "width {w}"); + if want > 1 { + assert_eq!(seen.load(Ordering::Relaxed), want, "width {w}"); + } + } + assert_eq!(par.degree(), 8); + stats::enable(false); + let s = stats::take(); + // only the narrowing to 3 is a sub-team; 8 and 9 are the token itself + assert_eq!(s.width[3], 1, "{s:?}"); + assert_eq!(s.width.iter().sum::(), 1, "{s:?}"); + }); + } + + /// Every converted family on a team of 8 whose row makes width 4 the widest that pays: the + /// regions narrow to 4 (counted) and the results equal the sequential ones. + #[test] + fn the_families_narrow_their_regions_and_match_sequential() { + use crate::dyn_stack::{MemBuffer, MemStack}; + use crate::linalg; + let ns = StdDuration::from_nanos; + // cheap up to 4, prohibitive from 5 on + let row = [(2, ns(10)), (4, ns(10)), (5, StdDuration::from_secs(1))]; + let tol = |a: MatRef<'_, f64>, b: MatRef<'_, f64>, what: &str| { + let scale = (0..a.ncols()).flat_map(|j| (0..a.nrows()).map(move |i| (i, j))).fold(1.0f64, |m, (i, j)| m.max(a[(i, j)].abs())); + let d = max_diff(a, b); + assert!(d <= 1e-10 * scale, "{what}: {d} against {scale}"); + }; + type Run = fn(Par<'_>) -> Vec>; + let families: [(&str, Run); 14] = [ + ("full pivoting lu", |par| { + let n = 200; + let mut a = rand_mat(n, n, 11); + let (mut p, mut pi, mut q, mut qi) = (vec![0usize; n], vec![0usize; n], vec![0usize; n], vec![0usize; n]); + let req = linalg::lu::full_pivoting::factor::lu_in_place_scratch::(n, n, par.shape(), Default::default()); + linalg::lu::full_pivoting::factor::lu_in_place(a.as_mut(), &mut p, &mut pi, &mut q, &mut qi, par, MemStack::new(&mut MemBuffer::new(req)), Default::default()); + vec![a] + }), + ("partial pivoting lu", |par| { + let n = 400; + let mut a = rand_mat(n, n, 12); + let (mut p, mut pi) = (vec![0usize; n], vec![0usize; n]); + let req = linalg::lu::partial_pivoting::factor::lu_in_place_scratch::(n, n, par.shape(), Default::default()); + linalg::lu::partial_pivoting::factor::lu_in_place(a.as_mut(), &mut p, &mut pi, par, MemStack::new(&mut MemBuffer::new(req)), Default::default()); + vec![a] + }), + ("column pivoting qr", |par| { + let n = 200; + let mut a = rand_mat(n, n, 13); + let bs = linalg::qr::no_pivoting::factor::recommended_block_size::(n, n); + let mut h = Mat::::zeros(bs, n); + let (mut q, mut qi) = (vec![0usize; n], vec![0usize; n]); + let req = linalg::qr::col_pivoting::factor::qr_in_place_scratch::(n, n, bs, par.shape(), Default::default()); + linalg::qr::col_pivoting::factor::qr_in_place(a.as_mut(), h.as_mut(), &mut q, &mut qi, par, MemStack::new(&mut MemBuffer::new(req)), Default::default()); + vec![a, h] + }), + ("hessenberg", |par| { + let n = 160; + let mut a = rand_mat(n, n, 14); + let bs = 8; + let mut h = Mat::::zeros(bs, n - 1); + let req = linalg::evd::hessenberg::hessenberg_in_place_scratch::(n, bs, par.shape(), Default::default()); + linalg::evd::hessenberg::hessenberg_in_place(a.as_mut(), h.as_mut(), par, MemStack::new(&mut MemBuffer::new(req)), Default::default()); + vec![a, h] + }), + ("tridiagonalization", |par| { + let n = 200; + let x = rand_mat(n, n, 15); + let mut a = &x + x.transpose(); + let bs = 8; + let mut h = Mat::::zeros(bs, n - 1); + let req = linalg::evd::tridiag::tridiag_in_place_scratch::(n, par.shape(), Default::default()); + linalg::evd::tridiag::tridiag_in_place(a.as_mut(), h.as_mut(), par, MemStack::new(&mut MemBuffer::new(req)), Default::default()); + vec![a, h] + }), + ("bidiagonalization", |par| { + let (m, n) = (300, 200); + let mut a = rand_mat(m, n, 16); + let bs = 8; + let mut hl = Mat::::zeros(bs, n); + let mut hr = Mat::::zeros(bs, n - 1); + let req = linalg::svd::bidiag::bidiag_in_place_scratch::(m, n, par.shape(), Default::default()); + linalg::svd::bidiag::bidiag_in_place(a.as_mut(), hl.as_mut(), hr.as_mut(), par, MemStack::new(&mut MemBuffer::new(req)), Default::default()); + vec![a, hl, hr] + }), + ("bunch kaufman, full pivoting", |par| { + use linalg::cholesky::lblt::factor::{LbltParams, PivotingStrategy}; + let n = 200; + let x = rand_mat(n, n, 17); + let mut a = &x + x.transpose(); + let mut sub = crate::diag::Diag::::zeros(n); + let (mut p, mut pi) = (vec![0usize; n], vec![0usize; n]); + let params = LbltParams { + pivoting: PivotingStrategy::Full, + ..>::auto() + }; + let req = linalg::cholesky::lblt::factor::cholesky_in_place_scratch::(n, par.shape(), params.into()); + linalg::cholesky::lblt::factor::cholesky_in_place(a.as_mut(), sub.as_mut(), &mut p, &mut pi, par, MemStack::new(&mut MemBuffer::new(req)), params.into()); + vec![a, sub.column_vector().as_mat().to_owned()] + }), + ("triangular solve", |par| { + let (n, k) = (100, 512); + let l = rand_mat(n, n, 18) + Mat::::identity(n, n) * (n as f64); + let mut x = rand_mat(n, k, 19); + linalg::triangular_solve::solve_lower_triangular_in_place(l.as_ref(), x.as_mut(), par); + vec![x] + }), + ("matvec, row split", |par| { + let a = rand_mat(1000, 1000, 20); + let v = rand_mat(1000, 1, 21); + let mut y = Mat::::zeros(1000, 1); + matmul(y.as_mut(), Accum::Replace, a.as_ref(), v.as_ref(), 1.0, par); + vec![y] + }), + ("matvec, column split", |par| { + let a = rand_mat(1000, 1000, 22); + let v = rand_mat(1000, 1, 23); + let mut y = Mat::::zeros(1, 1000); + matmul(y.as_mut(), Accum::Replace, v.transpose(), a.as_ref(), 1.0, par); + vec![y] + }), + ("block householder apply", |par| { + let (m, b, k) = (2000, 8, 256); + let basis = rand_mat(m, b, 24); + let factor = Mat::::from_fn(b, b, |i, j| if i == j { 1.5 } else if i < j { 0.1 } else { 0.0 }); + let mut target = rand_mat(m, k, 25); + let req = linalg::householder::apply_block_householder_on_the_left_in_place_scratch::(m, b, k); + linalg::householder::apply_block_householder_on_the_left_in_place_with_conj( + basis.as_ref(), factor.as_ref(), crate::Conj::No, target.as_mut(), par, MemStack::new(&mut MemBuffer::new(req)), + ); + vec![target] + }), + ("householder factor upgrade", |par| { + let (rows, b) = (4000, 16); + let essentials = rand_mat(rows, b, 26); + let mut factor = Mat::::from_fn(b, b, |i, j| if i == j { 1.25 } else { 0.0 }); + linalg::householder::upgrade_householder_factor(factor.as_mut(), essentials.as_ref(), b, 1, par); + vec![factor] + }), + ("sparse product", |par| { + use crate::sparse::{SparseColMat, Triplet}; + let n = 3000; + let triplets: Vec<_> = (0..n) + .flat_map(|j| (0..4).map(move |t| Triplet::new((j * 7 + t * 131) % n, j, ((j + t) as f64).sin()))) + .collect(); + let a = SparseColMat::::try_new_from_triplets(n, n, &triplets).unwrap(); + let c = crate::sparse::linalg::matmul::sparse_sparse_matmul(a.as_ref(), a.as_ref(), 1.0, par).unwrap(); + vec![c.to_dense()] + }), + ("tiled product (the non-x86 route)", |par| { + let (m, n, k) = (301, 263, 190); + let a = rand_mat(m, k, 27); + let b = rand_mat(k, n, 28); + let mut c = Mat::::zeros(m, n); + match par { + Par::Team(team) => unsafe { + gemm_tiled::( + team, m, n, k, 1.0, c.as_ptr_mut(), c.col_stride(), c.row_stride(), false, + a.as_ptr(), a.col_stride(), a.row_stride(), b.as_ptr(), b.col_stride(), b.row_stride(), + 0.0, 1.0, false, false, + ) + }, + Par::Seq => matmul(c.as_mut(), Accum::Replace, a.as_ref(), b.as_ref(), 1.0, Par::Seq), + } + vec![c] + }), + ]; + for (name, run) in families { + let want = run(Par::Seq); + let (got, s) = with_row(8, &row, |par| { + stats::enable(true); + let _ = stats::take(); + let got = run(par); + stats::enable(false); + (got, stats::take()) + }); + // narrowed to 4 somewhere, and never wider than 4 in a split + assert!(s.width[4] > 0, "{name}: no region narrowed to 4 ({s:?})"); + // and the family's own region (a split or a flat section, not only the GEMM kernel's + // fork) ran on 4 + assert!(s.split[4] + s.run[4] > 0, "{name}: no split or section of 4 ({s:?})"); + assert!(s.split.iter().skip(5).all(|&c| c == 0), "{name}: a split wider than 4 ({s:?})"); + for (k, (a, b)) in want.iter().zip(got.iter()).enumerate() { + tol(a.as_ref(), b.as_ref(), &alloc::format!("{name} output {k}")); + } + } + } + + /// `Never` keeps the whole call sequential, the compact-factor upgrade after the steps included, + /// even where the upgrade's own gate would pay (a 33 x 9 column pivoting QR at block 9 and a + /// bidiagonalization, on a row that makes every width up to 4 cheap). + #[test] + fn never_reaches_the_compact_factor_upgrades() { + use crate::dyn_stack::{MemBuffer, MemStack}; + use crate::linalg; + use crate::utils::cost::ParThreshold; + let ns = StdDuration::from_nanos; + let row = [(2, ns(10)), (4, ns(10)), (5, StdDuration::from_secs(1))]; + let s = with_row(8, &row, |mut par| { + stats::enable(true); + let _ = stats::take(); + { + use linalg::qr::col_pivoting::factor::*; + let (m, n) = (33, 9); + let mut a = rand_mat(m, n, 31); + let mut h = Mat::::zeros(n, n); + let (mut q, mut qi) = (vec![0usize; n], vec![0usize; n]); + let params = ColPivQrParams { par_threshold: ParThreshold::Never, ..>::auto() }; + let req = qr_in_place_scratch::(m, n, n, par.shape(), params.into()); + qr_in_place(a.as_mut(), h.as_mut(), &mut q, &mut qi, par.reborrow(), MemStack::new(&mut MemBuffer::new(req)), params.into()); + } + { + use linalg::svd::bidiag::*; + let (m, n) = (300, 200); + let mut a = rand_mat(m, n, 32); + let mut hl = Mat::::zeros(8, n); + let mut hr = Mat::::zeros(8, n - 1); + let params = BidiagParams { par_threshold: ParThreshold::Never, ..>::auto() }; + let req = bidiag_in_place_scratch::(m, n, par.shape(), params.into()); + bidiag_in_place(a.as_mut(), hl.as_mut(), hr.as_mut(), par, MemStack::new(&mut MemBuffer::new(req)), params.into()); + } + stats::enable(false); + stats::take() + }); + let parallel = s.run.iter().skip(2).sum::() + s.split.iter().skip(2).sum::() + s.gemm_tiled + s.gemm_team_kernel; + assert_eq!(parallel, 0, "{s:?}"); + } + /// `par_pays_time` is `ParShape::pays` without the shape. #[test] fn par_pays_time_is_the_shape_gate() { @@ -1124,7 +1425,7 @@ mod tests { } #[test] fn a_team_of_one_behaves_exactly_like_seq() { - use crate::utils::cost::{par_pays, rate, ParThreshold}; + use crate::utils::cost::{par_width, rate, ParThreshold}; // it reads the process-wide counters let _turn = team_turn(); threadteam::Exec::serial() @@ -1132,9 +1433,9 @@ mod tests { let par = Par::Team(team); assert_eq!(par.degree(), 1); // no gate pays on a single participant, however large the work - assert!(!par_pays::(&par, 1 << 40, rate::gemm_mac(), 1)); + assert_eq!(par_width::(&par, 1 << 40, rate::gemm_mac(), 1), 1); assert!(!crate::utils::cost::join_pays::(&par, 1 << 40, rate::gemm_mac(), 1)); - assert!(!ParThreshold::Derived.region_pays::(par.shape(), 1 << 40, rate::gemm_mac(), 1)); + assert_eq!(ParThreshold::Derived.region_width::(&par, 1 << 40, rate::gemm_mac(), 1), 1); // and the product equals the sequential one, with no team section at all let a = rand_mat(200, 150, 1); let b = rand_mat(150, 180, 2);