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 @@ -28,6 +28,7 @@ upstream `faer`.
| EVD eigenvector normalization + overflow-safe back substitution | `evd_real`/`evd_cplx` return eigenvectors scaled to unit 2-norm (a conjugate pair is scaled as one complex vector, matching `xGEEV`), and the back substitution is rescaled when an entry would exceed `sqrt(max)` instead of overflowing to `NaN`. A finite input that still overflows, or a Schur iteration that does not converge, returns `NoConvergence`. | KEEP | tenferro-rs#2131 | Focused test `eig_2131_tests` (f32/f64 Jordan block, conjugate pair with left and right vectors, complex input, sequential + recursive solver). Normalization runs at the `evd_real`/`evd_cplx` boundary, so `operator/eigen` keeps its previous Ritz-vector scaling while still getting the overflow-safe finite Schur eigenvectors. |
| 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. |
| 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
13 changes: 13 additions & 0 deletions faer/src/linalg/householder.rs
Original file line number Diff line number Diff line change
Expand Up @@ -146,6 +146,19 @@ pub fn upgrade_householder_factor<T: ComplexField>(
return;
}
let n = essentials.ncols();
// Below, this function threads its `V^H V` products, and its own `deep` predicate states when
// the reduction has the shape to pay for that: deep against its `n x n` output, over a work
// floor, with at least one row per thread. Below that floor the same products still open a
// region each -- and a blocked factorisation calls this function once per block, at a block
// size the recursion halves every level, so a `128 x 128` QR pays the fork-join cost of
// `n`-wide reductions ~50 times for work the floor itself calls too small to split. The
// decision here is the floor, asked once: no region at all when the products below would not
// split either.
let threads = crate::utils::thread::parallelism_degree(par);
let deep_rows = essentials.nrows().unbound().saturating_sub(n);
let above_floor =
deep_rows >= 2 * n && deep_rows >= 4 * threads && deep_rows * n * n >= 8 * 1024;
let par = if above_floor { par } else { Par::Seq };
let mut householder_factor = householder_factor;
let essentials = essentials;
assert!(householder_factor.nrows() == householder_factor.ncols());
Expand Down
Loading