Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions PATCHSET.md
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@ upstream `faer`.
| Unpivoted QR column skip | The skip lived only in `qr_in_place_unblocked`; `qr_in_place_blocked` has no threshold of its own and recurses into that core, so removing it there covers the tall and wide paths at once. The reflector is applied and the row advances for every nonzero column, matching LAPACK; the tall/wide `svd` path goes through it too. | KEEP | tlinalg-rs#28, tensor4all-rs#836 | Adopted. Tests: `test_near_dependent_columns_are_not_skipped` (tall `m x 3`, m=100/1000/8192) and `test_near_dependent_rows_are_not_skipped` (its transpose, wide `3 x n`, n=100/1000/8192), both at `<=1e-13`; measured 1e-16 to 5e-15. The wide fixture was already accurate before the change: at m=3 the removed threshold `16 (m-k) eps ||column||` is ~1e-14 and never tripped it (a wide skip's relative reconstruction error is ~`16 eps sqrt(m) sqrt(m/n)`), so the wide path needed no separate treatment — the same core change covers it. tlinalg's reported wide `3.33e-7` is a reconstruction-indexing bug in that test (entry `(i,j)` written to flat `j+3i` while the factor is column-major `i+3j`); faer's factors for that matrix are accurate to ~4e-16. Existing columns are unchanged; formerly skipped nonzero columns now pay reflector/update flops (the threshold norm bookkeeping is removed). No existing benchmark suite was run. Release smoke, 3-run mean, before→after: 100000×8 2.984→2.986 ms; 1000×1000 19.793→19.921 ms. |
| Row-split reduction in the compact-WY factor's deep half | `upgrade_householder_factor`'s flat branch (`prev_block_size < 8`) builds `T = diag(tau) + strict_upper(V^H V)`, and its second product reduces over every row below the basis' diagonal with an `n x n` output -- the same shape as the block apply's reduction, and serial for the same reason: 0.271 ms at eight threads at `m = 4096, n = 16`, against 0.040 with the rows split (6.7x; 2.96x at `m = 512`, 4.59x at `m = 1024`). The partials here are `n x n` each, tiny against the basis, and only the strict upper triangle is stored, so the split's overhead is negligible and no depth-against-width gate is needed -- the thread and work floors are enough. End to end, 8T pool: `householder_factor` 4096x256 b1 c64 15.72 -> 13.38 ms (f64 5.63 -> 4.84), `qr` c64 42.76 -> 37.00, and the 512x512 and 256x256 cells improve too (4.49 -> 4.14, 0.98 -> 0.92) with every 1T row unchanged. Test: `threaded_upgrade_matches_the_sequential_one` (4096x16, 2048x8, 4096x1, and 64x8 below the floors). | KEEP | tlinalg-rs#46 | Costs one `n * n * threads` buffer per call; the same scratch note as the row below. Own code: no OpenBLAS or MKL source. |
| No parallel region below the compact factor's own floor | `upgrade_householder_factor` threads its `V^H V` products, and its `deep` predicate already states when the reduction has the shape to pay for that (deeper than `2 * n`, at least one row per thread, over the `8 * 1024` work floor). Below that floor the products still opened a region each, and a blocked QR calls this function once per block -- at a block size the recursion halves every level -- so a `128 x 128` QR paid the fork-join cost of `n`-wide reductions about fifty times for work that same floor calls too small to split: the threaded call measured 4.2-4.6x *slower* than the sequential one at `n = 128` (1.05 vs 0.23 ms), 2.1x at 256, 1.3x at 384, and only broke even at 512 (0.96). Asking the floor once, before opening anything, gives 0.99 at 128, 0.79 at 256, 0.60 at 384 and 0.44 at 512 (f64), 0.72 / 0.50 / 0.30 (c64), and no region is opened below the floor at all. The tall shapes the ledger's other rows serve are untouched (`4096x256` f64 6.35 -> 6.33 ms, c64 16.67 -> 16.67), because their reductions are above the floor; `1024x64` goes from 1.4x slower than sequential to parity. No new constant: the predicate is the one the products below already use. | KEEP | tlinalg-rs#46 | Measured with an isolated driver on the call tlinalg's `qr` family makes (`qr_in_place` f64/c64, 128-768 square, 4096x256 / 1024x64 / 512x512), CPUs 4-11, best of 7. Why, and the candidates that lost: `docs/worklogs/2026-10-14-compact-factor-region-floor.md`. |
| Bidiagonalisation parallel hint | The lead that `bidiag_in_place_scratch`'s `_ = par; _ = params;` and an unread `BidiagParams::par_threshold` make the bidiagonalisation sequential, however the caller calls it. The scratch function's hint is not needed (the workspace is shape-sized, not thread-sized); the worker `bidiag_in_place` reads the threshold and forwards `par` into the trailing update. That floor makes `128 x 128` fully sequential and closes the threaded window at `256 x 256` at k=34; it is load-bearing: removing it costs 2.50x (f64) / 1.77x (c64) on `bidiag_in_place` at 128x128, 1.44x at 256x256 and 1.08-1.12x at 512x512, and end to end 1.81x (`svd_values`) / 1.52x (`svd_thin`) at 128x128. Keeping the caller's hint for the trailing compact-factor upgrades changed nothing (every row within noise, 0.94-1.03). | DROPPED | — | Isolated driver at f64/c64, 128/256/512 square plus 1024x64 and 4096x256, CPUs 0-3 (four cores), `RAYON_NUM_THREADS=8`, best of 7, interleaved A/B. Why, and the candidates that lost with their numbers: `docs/worklogs/2026-10-14-bidiag-parallel-hint.md`. |
| Row-split reduction in the block Householder apply | `tmp = V^H C` has an `N x K` output and reduces over every row of the target. With `N` small and the target long, the products inside the generic path were one serial reduction inside one task, so the column split (capped at `min(par, K, 4)`) could not use the depth: 4096x16 basis against a 4096x240 target measured 2.03 ms at 8T where OpenBLAS's `?gemm` on the same shape is 0.53-0.71. Splitting the rows is a split along the depth -- it adds `threads` partials of `N x K` to write, read and add -- so it is taken only while the reduction is *deeper than it is wide*: a sweep over `M` against `K` at eight threads puts the crossover at `M ~ 2 K` (four times deeper gains 3-4x; as wide as it is deep, which is the 512x512 trailing update at `M = 512`, `K = 496`, loses for the partials it pays for). With that gate, 8T pool, `householder_factor` 4096x256 b1 c64 17.44 -> 15.72 ms (f64 5.83 -> 5.54), `qr` c64 46.75 -> 42.68, `svd_thin` c64 77.11 -> 68.16, and the 512x512 and 256x256 cells stay at baseline; every 1T row is unchanged. Tests: `threaded_panel_matches_the_sequential_one` (2048x4, 4096x16, and two shapes below the thresholds) plus the two existing cases. | KEEP | tlinalg-rs#46 | Costs one `N * K * threads` buffer per call on the new path; the public scratch function is unchanged, so upstreaming this should carry it there. Own code: no OpenBLAS or MKL source is used or adapted. |
| Threaded rank-one Householder apply | The `N == 1` arm of `apply_block_householder_on_the_left_in_place_generic` is a sequential rank-one update per target column and it never read `par`. A tall blocked QR's panel is made almost entirely of those updates, so the panel was the serial part of it: at 4096x256 c64 the panel is 5.96 ms at 1T and 6.31 ms at 8T, against OpenBLAS's `?gemv`/`?geru` panel at 5.97 -> 2.58 ms on the same eight cores. When the caller asks for threads and the update is at least 8192 elements, the arm now falls through to the generic path, which partitions the target's columns and threads its two products. Through tlinalg's compact factor (`householder_factor`), pool row, 4T / 8T: 4096x256 f64 6.73 -> 5.77 and 6.67 -> 5.83 ms, 4096x256 c64 21.42 -> 18.75 and 19.50 -> 17.44 ms (11-14%), with 512x512, 128x128, `qr`, `svd_thin` and every 1T row unchanged within noise. Tests: `threaded_rank_one_matches_the_sequential_one` (f64 and c64; 2048x8, 16384x1, 64x8) and `rank_one_negates_its_own_reflector_column`. | KEEP | tlinalg-rs#46 | Own code: no OpenBLAS or MKL source is used or adapted, their timings are only the target. |
| Tall compact-WY trailing update scaling | The blocked compact QR applies its panel through `apply_block_householder_transpose_on_the_left_in_place_with_conj`, whose generic path partitions the target's columns with `n_tasks = min(parallelism_degree(par), K, 4)` and `rayon(par / n_tasks)` inside each. On a tall DRAM-bound panel (4096x16 basis, 4096x240 target, eight physical Zen 5c cores sharing 8 MiB of L3, 16.8 MB target) the update is 8.72 ms at `Par::Seq` and 2.45 ms at `Par::rayon(8)` = 3.56x, where OpenBLAS reaches 5.5x on the same cores for the same factorization; splitting the columns by hand saturates at 3.87x, so the ceiling is the kernel, not the decomposition. Candidates: partition rows for the `C -= V tmp` half, which is row-parallel and has the long dimension (rows >> columns, while the column partition caps at 4), and the DRAM-bound path of the inner matmuls. | TODO | tlinalg-rs#46 | Reproducer and the phase split: tlinalg-benchmark `docs/worklog.md`, 2026-10-12. Drives no code change in tlinalg until the kernel moves. |
Expand Down
91 changes: 91 additions & 0 deletions docs/worklogs/2026-10-14-bidiag-parallel-hint.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,91 @@
# The bidiagonalisation's parallel hint is not dropped

PR: [tensor4all/faer#11](https://github.com/tensor4all/faer/pull/11). Ledger
row: `Bidiagonalisation parallel hint`.

## Decisions

- **No code change; the experiment is dropped.** The lead closed here was left as
a hypothesis in `2026-10-14-compact-factor-region-floor.md`: that
`bidiag_in_place_scratch` discards its hint (`_ = par; _ = params;`) and that
`BidiagParams::par_threshold` is never read, so the bidiagonalisation is
sequential however the caller calls it. The first half is true of the
*scratch-size* function only, and correctly so — the workspace is one `n`-vector
per side and does not depend on the thread count. The worker,
`bidiag_in_place`, reads `params.par_threshold` on every column and forwards
`par` into the trailing update (the fused op's column split and both trailing
matmuls). There is no missing `par`. Measured `bidiag_in_place` at `Par::Seq`
and `Par::rayon(8)`, f64 and c64: the two are identical wherever the gate is
closed and differ only where it is open.
- **What actually makes `128 x 128` sequential.** The gate is
`(m - k - 1) * (n - k - 1) < params.par_threshold` and it *latches*: once the
trailing area drops below `192 * 256 = 49152`, `par` is `Par::Seq` for the rest
of the factorisation. At `128 x 128` the first column already has
`127 * 127 = 16129`, so the whole call is sequential by construction. That is
the same shape-of-the-work decision #9 makes for the compact factor ("no region
below the floor"), applied to the trailing update: below the floor the update
would fork and join for work the same predicate calls too small to split.
- **Rejected, with the numbers that rejected it: remove the gate** (forward the
caller's `par` unconditionally). It regresses every measured shape. Isolation
driver on `bidiag_in_place`, CPUs 0-3, `RAYON_NUM_THREADS=8`, best of 7,
interleaved A/B, min of 3 pairs, ms at `Par::rayon(8)`, gate → no gate:
`128 x 128` f64 **0.319 → 0.797** (2.50x), c64 **0.617 → 1.093** (1.77x);
`256 x 256` f64 2.328 → 3.345 (1.44x), c64 4.646 → 5.852 (1.26x);
`512 x 512` f64 15.05 → 16.29 (1.08x), c64 22.38 → 25.14 (1.12x). End to end
(`linalg::svd::svd` as tlinalg's provider calls it) at `128 x 128`: `svd_values`
f64 0.625 → 1.134 (1.81x), `svd_thin` f64 0.932 → 1.421 (1.52x), c64 1.419 →
1.879 (1.32x). The fork-and-join cost at these sizes is real; the gate is
load-bearing.
- **Rejected: keep the caller's hint for the two trailing
`upgrade_householder_factor` calls.** The latch reaches those calls too: for
`4096 x 256` it drops to `Par::Seq` at `k >= 243`, so both compact-factor
upgrades run sequential although #9's own `deep` floor would thread them.
Preserving the caller's hint there changed nothing — no row moved more than 6% in a
2-pair interleaved A/B, f64 and c64, from `128 x 128` to `4096 x 256` — so the
starved upgrades are not a measurable loss on these cores.
- **Upstream.** `faer/src/linalg/svd/` is byte-identical to the upstream mirror
(`git diff` against the GitHub mirror clone `sarah-ek/faer-rs` at `ddbd2ae` is
empty), so the numbers above are upstream's behaviour for this file and not a
fork divergence. The source of truth is Codeberg. No upstream binary was
measured separately (**none measured**).

## Verification conclusions and constraints

- Establishes: the campaign's flat `svd_values` at 128 and 256 is not a dropped
parallel hint. The bidiagonalisation is sequential because the work is below the
floor the same code already uses, and threading it does not pay on the measured
cores: at `128 x 128` `Par::Seq` and `Par::rayon(8)` are identical (0.319 / 0.319
ms f64, 0.617 / 0.617 c64), at `256 x 256` the threaded call is 4.7% slower
(2.223 → 2.328 ms f64), and only `512 x 512` pays (f64 16.06 → 15.05, c64 32.73
→ 22.38). End to end at the campaign's shapes, 1T / 8T, f64 / c64 (faer API;
`svd_values` with the QR params the provider uses): `128 x 128` 0.640 / 0.625 and
0.930 / 0.928, `256 x 256` 3.484 / 3.642 and 6.039 / 6.009. `svd_thin` scales
where the vector applications do, not where the bidiagonalisation does:
`128 x 128` 1.040 → 0.932 (f64) / 1.792 → 1.419 (c64), `256 x 256` 5.811 → 4.214
/ 11.58 → 7.549.
- Limits: **four cores.** Every number here is pinned to CPUs 0-3 (the campaign
holds 4-11 and their SMT siblings), so `Par::rayon(8)` runs eight workers on
four cores and the eight-physical-core crossover is not measured. A shape that
loses here could be a small win on eight idle cores; the campaign page's own
`svd_values` rows (1T vs 8T) say it is not for 128 and 256, but that is the
provider path, not this driver.
- Limits: these are faer-API numbers, not the tlinalg page's rows. At `128 x 128`
f64 my `svd_values` call is 0.64 ms while the page's `faer-1lane` row reads
1.41 ms, so ~0.8 ms sits in the provider (per-item scratch reuse and the output
push in `tlinalg/src/svd.rs`) rather than in `bidiag_in_place`. The page's
`svd_thin` row also carries a reconstruction check and a possible QR retry that
this driver does not, so this record does not attribute the page's gap to MKL
to faer's SVD core.
- Gate: `cargo test -p t4a-faer --lib` on this revision passes — 203 passed,
0 failed, 2 ignored. `sparse::linalg::qr::tests::
test_numeric_simplicial_qr_1_transpose` passes as well (five reruns); the
failure recorded in `2026-10-14-compact-factor-region-floor.md` did not
reproduce on `9de95f7`.
- Remains unverified, and these are hypotheses only:
- The eight-core crossover for the bidiag gate. The four-core A/B cannot see it.
- A blocked bidiagonalisation (LAPACK's `dgebrd` shape: a panel plus a `GEMM`
trailing update) is the candidate for the algorithm-level single-thread
difference, not threading. This record measures no blocked variant, and at
these shapes the faer-API call is already at or below the campaign page's MKL
row, so a faer-versus-MKL gap at `128 x 128` is a provider-level question
first.
16 changes: 13 additions & 3 deletions docs/worklogs/2026-10-14-compact-factor-region-floor.md
Original file line number Diff line number Diff line change
Expand Up @@ -68,9 +68,19 @@ Consumer: [tensor4all/tlinalg-rs#56](https://github.com/tensor4all/tlinalg-rs/pu
- `svd/bidiag.rs`: `bidiag_in_place_scratch` discards its hint (`_ = par; _ =
params;`) and `BidiagParams::par_threshold` is never read, so the bidiagonalisation
is sequential however the caller calls. Consistent with `svd_values` being flat from
1T to 8T (128x128 b1 f64 1.412 → 1.433 ms against MKL's 0.704 → 0.678).
- `evd/tridiag.rs` uses no parallel degree at all, and `eigh` is flat for the same
shape (1.26 → 1.24 ms against MKL's 1.08 → 0.586).
1T to 8T (128x128 b1 f64 1.412 → 1.433 ms against MKL's 0.704 → 0.678). Measured
and closed: `bidiag_in_place` (the worker) does read the threshold, and its floor
is load-bearing — see `2026-10-14-bidiag-parallel-hint.md`.
- `evd/tridiag.rs`: **corrected on re-reading, and the correction is the same lesson as the
bidiagonalisation lead above.** The scratch function ignoring its hint is the workspace sizing
it correctly (shape-sized, not thread-sized), and `tridiag_in_place` *does* read
`params.par_threshold` — `rem * rem / 2 < par_threshold` latches to `Par::Seq` below the
floor, which at 128x128 (`127² / 2 = 8064 < 49152`) closes it for the whole factorisation.
The flat `eigh` reading (1.26 → 1.24 ms against MKL's 1.08 → 0.586) is therefore *this
library's own floor doing its job*, not a dropped hint: MKL threads that size and faer
declines to. Whether the floor sits where it should is a measured question about the
crossover — the `M ~ 2K` episode in the ledger is the warning about moving such a constant
from the wrong cell — and it has not been measured. Not a lead to chase by reading.
- `lu/partial_pivoting/factor.rs` is *not* missing its hint (panel solve, trailing
matmul and recursion all take `par`, and its `par_threshold` is read), but
`lu_factor` 512x512 b1 f64 scales only 1.85x from 1T to 8T where MKL reaches 4.9x.
Expand Down
Loading