perf(lu): floor the blocked partial-pivoting recursion's trailing update

PATCHSET.md

@@ -33,6 +33,7 @@ upstream `faer`.
 | 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. |
+| Trailing-update size floor in the blocked partial-pivoting LU | `lu_in_place_recursion` halves its panel every level and ran, per level, a panel solve and a trailing product at the caller's `par` whatever their size. At `512 x 512` that is ~60 rayon regions, most of them a few thousand flops. Both products now take `Par::Seq` while `A11.nrows() * A11.ncols()` is not above `params.par_threshold` (`128 * 128`), the same field the swap stage already reads; the recursion still forwards the caller's `par`, so each level decides for its own trailing update. The floor is a sweep over 256x256, 512x512 and 1024x1024, not a single cell. Parent -> child, best of 4, min over 14 interleaved reps, f64, 4 threads (CPUs 0-3): 128x128 0.0980 -> 0.0582 (1.68x), 192x192 0.1597 -> 0.1120 (1.43x), 256x256 0.3108 -> 0.2550 (1.22x), 384x384 0.6046 -> 0.4806 (1.26x), 512x512 1.2799 -> 1.1673 (1.10x), 1024x1024 5.2400 -> 5.1188 (1.02x); 512x512 at 1T unchanged (1.8818 -> 1.8724). The vendor gap is not closed: the sequential leaf kernel and the level-0 products remain, and the swap stage's own region was deliberately left alone. Sixteen cores to eight cores of the campaign's page are not reproduced; four physical cores only. | KEEP | tlinalg-rs#46 | Own code: no vendor source is used or adapted. |
 | 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`. |
 | Native SVD / eigh direct output | At nbatch=1, stop writing U and the eigenvectors into an intermediate `Mat` and copying them out; write into the caller's output region instead. | TODO | tlinalg-rs#25 | Preserve output shape, initialization, and failure cleanup. |

docs/worklogs/2026-10-14-lu-trailing-update-floor.md

@@ -0,0 +1,94 @@
+# The blocked partial-pivoting LU's trailing-update size floor
+
+Ledger row: `Trailing-update size floor in the blocked partial-pivoting LU`.
+Consumer: the campaign's `lu_factor` row for `512 x 512` `f64`, which reads 2.83 ms
+at one thread and 1.53 ms at eight (1.85x) where MKL reaches 3.00 -> 0.61 (4.9x).
+
+## Decisions
+
+- **The level's two products are gated on the trailing update's size, nothing else
+  changes.** `lu_in_place_recursion` splits its panel in half at every level and runs,
+  per level, a panel solve (`solve_unit_lower_triangular_in_place`) and a trailing
+  product (`A11 -= A10 A01`). Both were handed the caller's `par` unchanged at every
+  depth, so the recursion opened a fork and a join around a trailing update that the
+  next level halves again: at `512 x 512` that is roughly sixty regions, most of them
+  a few thousand flops. They now take `Par::Seq` while `A11.nrows() * A11.ncols()` is
+  not above `params.par_threshold`. The recursion still forwards the caller's `par` to
+  its children, so each level decides for its own trailing update, and the panel
+  subtree that a level already ran is untouched by that level's decision.
+- **The floor is the existing `par_threshold`, not a new constant.** The measured
+  crossover is `128 * 128`, which is `PartialPivLuParams::par_threshold`'s default, and
+  the swap stage below already reads that field as the size at which parallelism stops
+  paying. One knob, and its doc comment ("threshold at which size parallelism should be
+  disabled") now describes both of its uses. A separate constant would need its own
+  sweep and its own justification for a value the sweep picks anyway.
+- **Measured by sweep, not by the cell that regressed.** A trailing-update floor was
+  swept over 0, 1024, 2048, 4096, 8192, 16384, 24576, 32768 and 65536 at four threads,
+  on `256 x 256` (where the region overhead is the largest share of the call),
+  `512 x 512` and `1024 x 1024`: 16384 is the best of the nine on all three, 8192 and
+  24576 are within 3%, and 65536 loses on all three. The sweep is what fixes the value;
+  `256 x 256` is not the shape the change was chosen on.
+- **Rejected: gating the swap stage's region too.** The swap stage is 5-10% of the call
+  and it too measures worse under the caller's `par` than under `Par::Seq`: in both
+  instrumented runs its bucket is larger at four threads than at one (0.1075 -> 0.1679
+  and 0.1586 -> 0.2249 ms). It has its own predicate,
+  `m * (ncols - n) > par_threshold`, which already refuses the top-level swap at
+  `512 x 512`. Tightening it is a different predicate with its own crossover, and in the
+  interleaved comparison below the same swap bucket moves to 0.1382 without a line of its
+  code changing, so this record could not separate a swap gate from that drift. Left
+  alone rather than moved on an unmeasured reading.
+- **Not rejected: the recursion itself.** Nothing about the panel/blocking structure
+  changed. The vendor gap at `512 x 512` is not closed by this: the remaining serial
+  part is the leaf kernel (`lu_in_place_unblocked`, 12-19% of the call and sequential by
+  construction) plus the level-0 matmul and solve, and the change removes only the
+  region overhead below the floor.
+
+## Verification conclusions and constraints
+
+- Establishes: the `512 x 512` scaling gap is partly this library's own scheduling, and
+  not only the vendor kernels. Instrumenting one `lu_in_place` (phase counters, that
+  build only) splits the time at four threads into trailing product and panel solve
+  (both bucketed by output size), the swap stage and the sequential leaf. The products
+  with an output of at least 16384 elements scale, the 4096..16383 bucket and below do
+  not: with the floor open vs forced at 16384, per-call mean over eight calls, the
+  smallest such mean over ten interleaved repetitions, `512 x 512` `f64`, four threads
+  (ms; each column is its own minimum, so the columns do not sum to the wall column) --
+  - trailing product, >= 16384: 0.2336 -> 0.2276; 4096..16383: 0.4648 -> 0.3982;
+    < 4096: 0.0480 -> 0.0121
+  - panel solve, >= 16384: 0.2487 -> 0.2390; 4096..16383: 0.0872 -> 0.0745;
+    < 4096: 0.0302 -> 0.0290
+  - swap stage: 0.1601 -> 0.1382; sequential leaf: 0.2974 -> 0.2814; total:
+    1.5960 -> 1.4134.
+  The two buckets the change touches account for roughly 0.12 of the 0.18 ms; the rest is
+  drift on buckets the change cannot reach, which is why the same comparison is also run
+  on the shipped pair of builds below.
+- Establishes, on the shipped code, two separately built binaries of the parent and the
+  child interleaved: best of four calls, smallest over fourteen interleaved repetitions,
+  `f64`, four threads (ms, parent -> child) -- `128 x 128` 0.0980 -> 0.0582 (1.68x),
+  `192 x 192` 0.1597 -> 0.1120 (1.43x), `256 x 256` 0.3108 -> 0.2550 (1.22x),
+  `384 x 384` 0.6046 -> 0.4806 (1.26x), `512 x 512` 1.2799 -> 1.1673 (1.10x),
+  `1024 x 1024` 5.2400 -> 5.1188 (1.02x). At one thread the two builds are the same
+  code path and measure the same (`512 x 512`: 1.8818 -> 1.8724); at two threads
+  `512 x 512` is 1.4943 -> 1.4308. The win is entirely in the deep levels, so it is
+  largest where the deep levels are the largest share of the call.
+- Limits: **four physical cores**, CPUs 0-3 (pinned throughout; the campaign holds
+  4-11 and their siblings), so the eight-physical-core row of the campaign's page is
+  not reproduced here and the four-thread crossover is not the eight-thread one. The
+  floor is a constant that a machine with more cores would want lower, and this
+  measurement cannot say where.
+- Limits: with an eight-thread rayon pool asked for eight workers on those same four
+  cores, `512 x 512` `f64` takes 57 ms on the parent -- fifty times slower, not a
+  marginal effect. It is reproducible and it is in the parent build, so it is not this
+  change; it is recorded so that the 8T column above is not read as an eight-core
+  result.
+- The focused test is
+  `linalg::lu::partial_pivoting::factor::tests::threaded_lu_matches_the_sequential_one`
+  (96, 128, 192 below `128 * 128` at the top level, 256 and 512 above it, sequential and
+  four-thread factors compared entry-wise, permutations and transposition counts
+  compared exactly). The gate schedules the same arithmetic, so the test cannot fail on
+  a wrong constant; it fixes that both sides of the floor agree.
+- Gate: `cargo test -p t4a-faer --lib` on this revision passes, 205 passed, 0 failed,
+  2 ignored. That includes
+  `sparse::linalg::qr::tests::test_numeric_simplicial_qr_1_transpose`, which an earlier
+  record blamed on the parent revision: it passes on this parent too, so the blame was
+  misattributed, as `2026-10-14-bidiag-parallel-hint.md` had already found.

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

@@ -100,11 +100,23 @@ pub(crate) fn lu_in_place_recursion<I: Index, T: ComplexField>(
 			A.rb_mut().split_at_mut(block_size, block_size);
 		let A00 = A00.rb();
 		let A10 = A10.rb();
+		// The panel solve and the trailing product below are this level's whole
+		// parallel work, and the block size halves every level of the recursion, so
+		// past a few levels they are a fork and a join around a handful of thousands
+		// of flops. Below the size floor the swap stage already uses, take
+		// `Par::Seq`; the recursion keeps the caller's `par`, so every level decides
+		// for its own trailing update. The crossover comes from a sweep over 256x256,
+		// 512x512 and 1024x1024 threaded, not from one shape.
+		let par_here = if A11.nrows() * A11.ncols() > params.par_threshold {
+			par
+		} else {
+			Par::Seq
+		};
 		{
 			linalg::triangular_solve::solve_unit_lower_triangular_in_place(
 				A00.rb(),
 				A01.rb_mut(),
-				par,
+				par_here,
 			);
 		}
 		linalg::matmul::matmul(
@@ -113,7 +125,7 @@ pub(crate) fn lu_in_place_recursion<I: Index, T: ComplexField>(
 			A10.rb(),
 			A01.rb(),
 			-one::<T>(),
-			par,
+			par_here,
 		);
 		n_trans += lu_in_place_recursion(
 			A.rb_mut().get_mut(block_size..m, ..),
@@ -402,4 +414,47 @@ mod tests {
 			assert!(p.inverse() * L * U ~ A);
 		}
 	}
+	fn factor_lu(A: &Mat<f64>, par: Par) -> (Mat<f64>, Vec<usize>, usize) {
+		let n = A.nrows();
+		let mut LU = A.cloned();
+		let perm = &mut *vec![0usize; n];
+		let perm_inv = &mut *vec![0usize; n];
+		let params: Spec<PartialPivLuParams, f64> = default();
+		let stack =
+			&mut MemBuffer::new(lu_in_place_scratch::<usize, f64>(n, n, par, params));
+		let stack = MemStack::new(stack);
+		let (info, _) = lu_in_place(LU.as_mut(), perm, perm_inv, par, stack, params);
+		(LU, perm.to_vec(), info.transposition_count)
+	}
+	/// The two products at a recursion level take `Par::Seq` while the trailing update is smaller
+	/// than `par_threshold`, instead of the caller's hint. That is a choice of how to schedule the
+	/// same arithmetic, so both sides of the floor have to return the same factorization: `96`,
+	/// `128` and `192` are below `128 * 128` at the top level, `256` and `512` above it, and every
+	/// level in between halves down through it.
+	#[test]
+	fn threaded_lu_matches_the_sequential_one() {
+		let rng = &mut StdRng::seed_from_u64(0);
+		let pool = rayon::ThreadPoolBuilder::new()
+			.num_threads(4)
+			.build()
+			.unwrap();
+		let approx_eq = CwiseMat(ApproxEq {
+			abs_tol: 1e-12,
+			rel_tol: 1e-12,
+		});
+		for n in [96, 128, 192, 256, 512] {
+			let A = CwiseMatDistribution {
+				nrows: n,
+				ncols: n,
+				dist: StandardNormal,
+			}
+			.rand::<Mat<f64>>(rng);
+			let (seq, perm_seq, trans_seq) = factor_lu(&A, Par::Seq);
+			let (par, perm_par, trans_par) =
+				pool.install(|| factor_lu(&A, Par::rayon(4)));
+			assert!(seq.as_ref() ~ par.as_ref());
+			assert!(perm_seq == perm_par);
+			assert!(trans_seq == trans_par);
+		}
+	}
 }