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 @@ -30,6 +30,7 @@ upstream `faer`.
| 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`. |
| 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. |
| Row-split gate for the block Householder apply's top block | The row split slices `tmp = V^H C`'s rows, but the top `n` rows are applied through the unit-triangular head's own `n x n` factor: a chunk narrower than `n` would need that factor column-sliced, and the trapezoid that leaves is a structure `triangular::matmul` cannot state. The split's predicate now also requires `ceil(M / threads) >= n`, so the head is always applied whole; a shape whose chunk is narrower than its panel takes the products below. Every shape the split's own measurements used (`M = 4096` against `K = 240`, `M = 2048 x 4`) has `chunk >= n` and is unchanged. Test: `threaded_panel_narrower_than_the_row_chunk_matches_the_sequential_one`. | KEEP | tensor4all/faer#12 | Fixes the `lhs_ncols == rhs_nrows` assertion that aborted `qr_in_place` at 8 threads (64 against 32) and blocked the campaign's `householder_apply` page. The row split is fork-only, so this carries no upstream feedback. |
| 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. |
| Typed uninitialized overwrite destination | A `MatUninitMut<T>`-style destination that never reads the prior output, accepted by `matmul_with_conj`. | TODO | strided-rs#198 | strided-rs currently returns a Faer-only `Unsupported`. |
Expand Down
71 changes: 71 additions & 0 deletions docs/worklogs/2026-10-14-row-split-top-block-gate.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,71 @@
# The block Householder row split's top block

PR: [tensor4all/faer#13](https://github.com/tensor4all/faer/pull/13). Ledger row:
`Row-split gate for the block Householder apply's top block`.
Consumer: [tensor4all/faer#12](https://github.com/tensor4all/faer/issues/12), the
campaign blocker the `householder_apply` page could not be recorded past.

## Decisions

- **Take the row split only while one chunk swallows the top `n` rows.** The split sums
`tmp = V^H C` from `n_chunks` partials over the target's rows. The first `n` rows of `V`
are its unit-lower triangular head, and their contribution is applied through that head's
own `n x n` factor (`essentials_top.transpose()`, `UnitTriangularUpper`); the rows below
go through the rectangular `essentials_bot`, which slices correctly. A chunk whose rows are
a strict *subrange* of `[0, n)` would need the head factor's columns `lo..top_hi` -- a
`top_hi x (top_hi - lo)` unit-upper trapezoid, a shape `triangular::matmul` cannot state
because its precondition wants a square triangular operand. The predicate asked only for
depth (`M >= 2K`), a thread floor and a work floor, so any shape with
`ceil(M / threads) < n` handed the full `n x n` factor to a shorter row chunk and tripped
`lhs_ncols == rhs_nrows` inside `qr_in_place`. One more conjunct,
`M.div_ceil(threads) >= n`, is the whole change; no product is computed differently.
- **The owner is the library, not the harness.** The abort is reachable with no harness at
all: `qr_in_place` on a 2048x2048 `f64` at 8 threads with its *own*
`recommended_block_size` (64) aborts, and so does 4096x4096 at 4 and 8 threads, before any
`householder_apply` call runs. The harness's own fixture is legal: it passes
`block_size = min(m, n)` to `qr_in_place_scratch`/`qr_in_place`, whose only contract on that
argument is `> 0` and `Q_coeff.ncols() == min(A.nrows(), A.ncols())`. The campaign's
256x256 case is exactly this shape, and the assertion fires *inside* `qr_in_place`, not in
the apply that follows it. The row split is fork-only (`fa22a69`, `c061812`, `c8fea58`);
upstream's generic apply has only the column split, whose operands cannot mismatch. So this
is a fork defect and no upstream revision was measured for it; there is nothing to report
upstream.
- **Rejected: slice the head factor's columns instead of gating.** The correct partial is
`(V_top^H)[:, lo..top_hi] * C[lo..top_hi, :]`; expressing it needs the `lo x (top_hi - lo)`
strict-upper block and the `(top_hi - lo) x (top_hi - lo)` unit upper block as two products
into disjoint row ranges of the partial, plus an explicit zero fill of the rows the head
leaves empty. That is 15-20 lines and a new invariant ("the tail of every partial stays
zero") to buy parallel work on the top `n` rows -- a range the ledger's own numbers do not
show the split paying for. Rejected for now; the gate is one predicate.
- **Not rejected: the row split itself.** Every shape the split was measured on has
`chunk >= n` and is untouched by the gate: the tall class is a 4096-row target with a 16-wide
basis (`chunk = 512`, `K = 240`) and the recorded tests use 2048x4 and 4096x16, so none of
the split's recorded wins is removed. That is a shape argument, not a re-measurement: no
timing was taken for this change.

## Verification conclusions and constraints

- Establishes: a threaded no-pivot QR is correct again at every block size, and the failing
shape class computes the same map as the sequential call.
- The focused test is
`linalg::householder::tests::threaded_panel_narrower_than_the_row_chunk_matches_the_sequential_one`
(a 256x128 basis with a 64-column target at 4 threads, `chunk = 64 < n = 128`). On the
parent it fails with the same assertion as the issue (`lhs_ncols = 128`, `rhs_nrows = 64`);
with the gate it passes. An isolated `qr_in_place` sweep confirms the broader fix:
256x256, 2048x2048 and 4096x4096 at 1/2/4/8 threads and block sizes 64/128/256 no longer
abort, where 2048x2048 at 4 and 8 threads aborted on the parent with the *recommended* block
size.
- Reduces the campaign blocker: the harness's `verify` over the whole declared `householder`
shape grid (4, 8, 16, 32, 64, 256, 1024x64, 4096x256; `f64`/`c64`; both orders; 1, 4 and 8
threads; fork and upstream rows) reports "all comparisons within tolerance" and exits 0.
That is the check `record_run.py` makes before publishing, so a page for a population
containing `householder_apply` is recordable again.
- Limits: the gate is a *restriction*, not a better top block. A deep reduction with a chunk
narrower than its panel now runs its top rows on the column path; that shape class has not
been measured, and no win is claimed or lost for it. The whole-fork gate was
`cargo test -p t4a-faer --lib` (204 passed, 0 failed, 2 ignored); it carries no pre-existing
failure. In particular `sparse::linalg::qr::tests::test_numeric_simplicial_qr_1_transpose`,
which the issue names as failing on 22f23cf, passes here on both 22f23cf and 9de95f7 -- its
calls are all `Par::Seq`, so it cannot reach the split at all. The second reproducer in the
issue was therefore misattributed; the isolated `qr_in_place` sweep above is what
establishes the library ownership.
20 changes: 19 additions & 1 deletion faer/src/linalg/householder.rs
Original file line number Diff line number Diff line change
Expand Up @@ -624,11 +624,18 @@ fn apply_block_householder_on_the_left_in_place_generic<
// trailing update, `M = 512`, `K = 496`) loses for the partials it pays for. The thread and
// work floors keep small problems on the products below.
let threads = crate::utils::thread::parallelism_degree(par);
// The split below slices the reduction's rows, but the top `n` rows carry the unit
// diagonal and are applied through their own `n x n` factor -- a chunk whose rows
// are a *subrange* of those `n` would need that factor column-sliced, and the
// triangular kernel cannot express the trapezoid that leaves. So the split is taken
// only while one chunk swallows the whole top block. Below that the products further
// down stay exactly as they were.
let row_split = cfg!(feature = "rayon")
&& threads > 1
&& M.unbound() >= 2 * K.unbound()
&& M.unbound() >= 4 * threads
&& M.unbound() * K.unbound() >= 8 * 1024;
&& M.unbound() * K.unbound() >= 8 * 1024
&& M.unbound().div_ceil(threads) >= N.unbound();
if row_split {
#[cfg(feature = "rayon")]
{
Expand Down Expand Up @@ -1271,6 +1278,17 @@ mod tests {
check_panel(128, 8, 16);
}

/// The row split slices the reduction's rows, but the top `n` rows are applied through one
/// `n x n` factor: the split is only valid while a single chunk carries that whole top block.
/// A panel wider than the row chunk does not -- and it used to trip the triangular shape
/// assertion (64 factor columns against a 32-row chunk) from `qr_in_place`, which is the abort
/// [tensor4all/faer#12](https://github.com/tensor4all/faer/issues/12) reports. Both sides of
/// the threshold must compute the same map.
#[test]
fn threaded_panel_narrower_than_the_row_chunk_matches_the_sequential_one() {
check_panel(256, 128, 64);
}

#[test]
fn threaded_rank_one_matches_the_sequential_one() {
check_f64(2048, 8);
Expand Down
Loading