From 033e346980752b5e2838122f5c5d82c749400dff Mon Sep 17 00:00:00 2001 From: Hiroshi Shinaoka Date: Mon, 9 Nov 2026 18:03:13 +0900 Subject: [PATCH 1/6] perf(lu): factor small square row-major matrices through the full-pivot dual (#314) Full-pivoting `lu_in_place` took the shared row-major copy from 64 * 64 elements up, ahead of its own transpose dual. On square matrices below 128 the dual runs within 1.05x of a column-major input (f64 64x64 one thread: 1.01x on Apple M5 Max, 1.03x on Zen 5) while the copy costs 1.24x / 1.32x, so those shapes now skip the copy (`row_major_takes_the_copy`). From 128 up the dual's row scaling walks a strided row every step (1024 one thread: 0.22 -> 3.26 ms for that stage alone, the update and pivot search equal), and on tall shapes the long swap turns strided and the trailing block wide (1.4-2.1x), so every other shape keeps the copy. Tests: `faer/tests/full_piv_lu_row_major_no_copy.rs` bounds the largest single allocation below m * n * 8 on 64/96/127 square and sees the copy on 128x128 and 256x32; the full-pivoting arm of the row-major adaptation test asks the new guard for bitwise versus rounding agreement and covers 127x127. Co-Authored-By: Claude Opus 5.5 --- PATCHSET.md | 2 +- .../2026-11-09-full-piv-lu-row-major-dual.md | 69 +++++++++++ faer/src/linalg/lu/full_pivoting/factor.rs | 33 +++++- faer/src/linalg/mod.rs | 29 ++++- faer/tests/full_piv_lu_row_major_no_copy.rs | 112 ++++++++++++++++++ 5 files changed, 236 insertions(+), 9 deletions(-) create mode 100644 docs/worklogs/2026-11-09-full-piv-lu-row-major-dual.md create mode 100644 faer/tests/full_piv_lu_row_major_no_copy.rs diff --git a/PATCHSET.md b/PATCHSET.md index 562ddf0a..43b4a817 100644 --- a/PATCHSET.md +++ b/PATCHSET.md @@ -121,7 +121,7 @@ gate has changed and two kinds of adoption are now in the ledger. | Empty complex tridiagonal EVD with eigenvectors | `tridiagonal_self_adjoint_evd` had no empty-input return, unlike `self_adjoint_evd`. For a complex type with eigenvectors requested, its phase loop wrote `scale[0]` into an empty buffer and panicked; real types and the values-only path did not reach that write. It now returns `Ok(())` for `n == 0`. The function is upstream faer code (sarah-quiñones); nothing was reported upstream. | KEEP | tensor4all/faer#31 | Correctness only. Found by the #31 LAPACK validation harness (`lapack_tridiagonal_evd`, 36 panicking calls per complex type). Test `empty_tridiagonal_with_vectors` (f64/c64/c32) fails before the change and passes after. Upstream report: [`upstream-reports/faer-021-tridiag-evd-empty-complex`](upstream-reports/faer-021-tridiag-evd-empty-complex/). | | Row-major QR layout adaptation | `qr_in_place` factors a matrix whose columns are not the contiguous direction through a column-major copy: the panel is not layout-aware (`make_householder_in_place` reads the pivot column and scales it, and the rank-one update beside it walked the block with `z!`), so a row-major matrix cost one cache line per element -- at `4096x256` f64 the panel was 25x its column-major cost and the whole deficit, and the factorization did not scale (47.8 / 42.9 / 43.0 ms at 1T / 4T / 8T). `qr_in_place` now copies into a column-major buffer, factors that, and copies the factor back, which is what LAPACK's row-major interface pays internally; the arithmetic is the copy's, so the factor is the column-major factor of the same matrix bit for bit. Tiled copy; with threads it is split along the destination's major dimension in whole 16-element tiles (columns into the buffer, rows on the write-back, so a thread writes whole lines of its own: tensor4all/faer#92), and the work buffer is allocated uninitialised and initialised by the threads that copy into it, so its pages are first touched in parallel rather than by a serial `vec![zero; m * n]` on the caller (`linalg::column_major_copy`, shared by the four row-major adaptations). One `m * n` heap buffer per call, because `qr_in_place_scratch` cannot know the strides. The copy has no cost gate yet (tensor4all/faer#92 part 2, after #91's dispatch measurements). Paired, best of 21, CPUs 4 / 4-7 / 4-11, fork `bfe4041` -> `044f947`, C order ms (before -> after, F order after): `4096x256` f64 47.83 / 42.85 / 42.96 -> 19.44 / 10.22 / 8.39 (15.77 / 8.37 / 7.23), `1024x64` f64 1.090 / 2.687 / 2.919 -> 0.587 / 0.613 / 0.678, `512x512` f64 9.44 / 9.56 / 12.73 -> 5.24 / 3.11 / 3.24, `256x256` f64 1.455 / 1.596 / 2.315 -> 0.921 / 0.833 / 0.904. Against the campaign's MKL C-order column (27.73 / 8.66 / 8.20) `4096x256` f64 goes from 1.72x / 4.95x / 5.24x behind to 0.70x / 1.18x / 1.02x. Test: `test_row_major_factor_is_the_column_major_one` -- the factor, `tau` and rank bit for bit at f64 and c64, tall/wide/square, sequential and threaded, backwards strides, and a zero leading column. A one-row matrix is left on the direct path (it eliminates at most one column). The `perf-warn` that told row-major callers to transpose is gone. | KEEP | tlinalg-rs#46 | Heap allocation inside a public call, the same boundary the fork's two row splits already cross. Own code: no OpenBLAS or MKL source is used or adapted. Exactly the adapted shapes (tall or square) are the column-major factorisation bit for bit; the guards' exclusions are not. Why, and what lost: `docs/worklogs/2026-10-15-qr-row-major-layout.md`. | | Row-major col-pivoting QR | `qr::col_pivoting::factor::qr_in_place` computes the norm of every one of the `n` columns before it starts and re-reduces them per step, and its incremental-norm SIMD path is gated on `A.row_stride() == 1` and unwraps a column-major view, so a row-major matrix pays one cache line per element on every shape. It now factors a column-major copy and writes the factor back, wide included. Paired, best of 21, CPUs 4 / 4-7 / 4-11, `f64`, C order before -> after (F order after): `256x256` 58.96 / 59.63 / 59.61 -> 1.52 / 1.40 / 1.19 (1.45 / 1.33 / 1.54), `1024x64` 22.23 / 22.29 / 22.45 -> 0.726 / 0.520 / 0.509 (0.472 / 0.522 / 0.536); `c64` `256x256` one thread 64.40 -> 3.29, `1024x64` 24.69 -> 1.64. So the C/F ratio at `256x256` one thread is 38.8x -> 1.05 and the C order is *faster* than the F order at eight threads. Nothing in the campaign records a column-pivoted QR, so this is the record's own driver. Test: `linalg::tests::test_row_major_adaptation_matches_the_column_major_one`. | KEEP | tlinalg-rs#46 | Heap allocation inside a public call, as for the no-pivoting row. Why, and what lost: `docs/worklogs/2026-10-15-row-major-adaptation-lu-colpiv.md`. | -| Row-major LU, partial pivoting | `lu::partial_pivoting::factor::lu_in_place` walks each pivot column twice per step (search, then scaling) and swaps elements *inside* a column, which is a cache line per element when the columns are not contiguous. It now factors a column-major copy and writes the factor and its permutations back. Paired, best of 21, `f64`, C order before -> after (F order after): `512x512` 12.14 / 7.02 / 6.34 -> 3.37 / 2.01 / 1.89 (3.04 / 1.88 / 1.81) at 1T / 4T / 8T, `1024x1024` 59.80 / 27.62 / 22.93 -> 22.91 / 10.04 / 8.21 (19.89 / 8.19 / 6.65), `256x256` one thread 1.530 -> 0.486; `c64` `512x512` one thread 18.08 -> 9.50. Test: `linalg::tests::test_row_major_adaptation_matches_the_column_major_one` -- factor, permutations and transposition count bit for bit. The campaign's `lu` cell takes the composite, so `to_owned` below is what moves it. | KEEP | tlinalg-rs#46 | The full-pivoting flavour is deliberately **not** adapted: it already factors a row-major input through `mat.transpose_mut()`, and the copy measured a no-op there (512x512 f64 one thread 13.758 against 13.712 C order). Why: `docs/worklogs/2026-10-15-row-major-adaptation-lu-colpiv.md`. | +| Row-major LU, partial pivoting | `lu::partial_pivoting::factor::lu_in_place` walks each pivot column twice per step (search, then scaling) and swaps elements *inside* a column, which is a cache line per element when the columns are not contiguous. It now factors a column-major copy and writes the factor and its permutations back. Paired, best of 21, `f64`, C order before -> after (F order after): `512x512` 12.14 / 7.02 / 6.34 -> 3.37 / 2.01 / 1.89 (3.04 / 1.88 / 1.81) at 1T / 4T / 8T, `1024x1024` 59.80 / 27.62 / 22.93 -> 22.91 / 10.04 / 8.21 (19.89 / 8.19 / 6.65), `256x256` one thread 1.530 -> 0.486; `c64` `512x512` one thread 18.08 -> 9.50. Test: `linalg::tests::test_row_major_adaptation_matches_the_column_major_one` -- factor, permutations and transposition count bit for bit. The campaign's `lu` cell takes the composite, so `to_owned` below is what moves it. | KEEP | tlinalg-rs#46 | The full-pivoting flavour is deliberately **not** adapted: it already factors a row-major input through `mat.transpose_mut()`, and the copy measured a no-op there (512x512 f64 one thread 13.758 against 13.712 C order); `10ee1c9` adapted it after all, behind the same `row_major_copy_pays` floor. Since tensor4all/faer#314 the square row-major matrices below 128 skip that copy and go through the transpose dual, which runs within 1.05x of a column-major input there (64x64 f64 one thread: 1.24x through the copy, 1.01x through the dual on Apple M5 Max; 1.32x and 1.03x on Zen 5); from 128 up and on every tall shape the copy stays, because the dual's strided row scaling (and, tall, its strided long swap and wide trailing block) costs more than the copy. Guard `full_pivoting::factor::row_major_takes_the_copy`; tests `faer/tests/full_piv_lu_row_major_no_copy.rs` (no `m * n` allocation on the dual shapes, the copy still seen on 128x128 and 256x32) and the full-pivoting arm of `test_row_major_adaptation_matches_the_column_major_one` (127x127 added). Why: `docs/worklogs/2026-11-09-full-piv-lu-row-major-dual.md`. Why: `docs/worklogs/2026-10-15-row-major-adaptation-lu-colpiv.md`. | | Row-major `to_owned` | `MatRef::to_owned` builds its result with `Mat::from_fn`, which fills column by column, so a row-major source is read with a stride of `ncols` and every element costs its own cache line: `512x512` f64 1.01 ms against 0.057, `1024x1024` 6.36 against 0.36, `4096x256` 5.08 against 0.37, `256x4096` 4.30 against 0.36. It now stages a tile when the source's columns are not contiguous and the source is already canonical. After: 0.275 / 2.007 / 2.045 / 1.774 ms. This is what the campaign's `lu` and `qr` cells pay -- a probe at `n = 512` f64 one thread C order reads `lu` 6.1417 -> 4.156 and 4.5329 -> 3.033 at eight threads, `qr` 12.5143 -> 11.185 -- and `full_piv_lu` 12.95 -> 12.63, `cholesky`, `eigh` and `svd_values` unchanged, because those three copy through their own routes rather than `to_owned`. | KEEP | tlinalg-rs#46 | The tile adds a `Mat::zeros` pass and is serial (`to_owned` has no `par`). Gated on `T::IS_CANONICAL`: `canonical()` reinterprets memory and does not apply `Conj::apply`, so a conjugated source (`A.adjoint().to_owned()`, which the sparse Cholesky tests use) keeps the elementwise path -- the ungated attempt broke five of them. The attempted generic `cloned` tile is **dropped**: initializing every destination entry by cloning the first source element amplifies heap-owned payloads by `m * n` before overwrite. `cloned` retains the original elementwise implementation; `test_cloned_does_not_replicate_the_first_heap_value` checks that each source is cloned once. The `Option` tile still moves each cloned numeric value into `to_owned`/factorization destinations without a second clone. Why: `docs/worklogs/2026-10-15-row-major-adaptation-lu-colpiv.md`. | | Row-major `copy_from_triangular_lower` | `MatMut::copy_from_triangular_lower` (what `Llt::new` and `Ldlt::new` copy their input with) ran the elementwise zip, which for a row-major source into a column-major destination reads with a stride of the source's row length: `1024x1024` f64 one thread 4.16 ms against 0.30 ms from a column-major source. It now copies each 16-column band as its diagonal tile plus the rectangle under it through `copy_across_layouts_block`, the tile `to_owned` uses: 0.65 ms. The campaign's `cholesky` cell, 1024 f64 8T, three processes per arm against MKL in the same process: row-major 5.445 -> 2.236 ms (MKL 3.33 / 3.39), column-major 1.674 -> 1.649 (MKL 2.65 / 2.67). Tests: `linalg::tests::tiled_triangular_copy_writes_what_the_elementwise_one_does` and `linalg::tests::llt_and_partial_piv_lu_factor_both_orders_bit_for_bit` (the factor of the same matrix through the two orders, bit for bit). | KEEP | faer-benchmark slide re-take (1024, 8T) | Serial, like `to_owned`: the method has no `par`. A conjugated source keeps the elementwise path. `copy_from` (the full-matrix copy) still runs elementwise: 8.6 ms row-major against 0.47 ms at `1024x1024`; it is not on the `cholesky` or `lu` path. Why: `docs/worklogs/2026-11-06-row-major-triangular-copy.md`. | | 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. | diff --git a/docs/worklogs/2026-11-09-full-piv-lu-row-major-dual.md b/docs/worklogs/2026-11-09-full-piv-lu-row-major-dual.md new file mode 100644 index 00000000..0e0f30da --- /dev/null +++ b/docs/worklogs/2026-11-09-full-piv-lu-row-major-dual.md @@ -0,0 +1,69 @@ +# Full-pivoting LU: small square row-major inputs skip the copy (tensor4all/faer#314) + +## Decisions + +- `full_pivoting::factor::lu_in_place` used to take the shared row-major copy + (`linalg::row_major_copy_pays`, at least `64 * 64` elements, tall or square) before + its own transpose dual (`lu_in_place_unblocked(.., transpose = true)`) was even tried. + It now skips the copy for **square matrices below 128** and keeps it everywhere else + (`row_major_takes_the_copy`). +- The rule is the issue's acceptance criterion applied shape by shape: the dual has to + run within 1.05x of a column-major input. Square `64..=127` passes on both machines + measured (Apple M5 Max, this change: f64 C/A 1.009 / 1.039 at 64 and 1.007 / 1.005 + at 128, one / four threads; in an ungated scan, 64-120 square at 0.98-1.04 in f64 and c64; Zen 5, the + issue's earlier comment: n = 64 and 96 at 0.94-1.04 over f64/c64/f32 at 1T and 8T). + Square `128` is mixed on Zen 5 (f64 1T 1.114) and from 256 up the dual leaves the band + on at least one machine (M5 four threads 1.087 at 512, 1.089 at 1024; Zen 5 + 1.08-1.27), so the boundary sits below 128. On M5 one thread the square dual stays + at 1.03-1.04 up to 1024, so a single-thread-only rule could go further; not taken. +- Off the square the band ends early (`80 x 64` f64 1.06x, `128 x 64` 1.11-1.14x on + M5, the same ungated scan), so near-square shapes keep the copy. +- No change to the dual kernel in this slice, and no change to `row_major_copy_pays` + or to the other families that share it. + +## Where the dual loses (stage breakdown, M5 Max, f64, one thread) + +Per-stage timers inside `lu_in_place_unblocked` (local patch, not committed): + +- **Square, n >= 256**: the rank-one update with the fused pivot search costs the same + through both layouts (1024: 104.4 ms column-major, 103.7 ms dual). The whole gap is + the dual's row scaling `A[(k, j)] *= inv`, which walks a strided row every step, + one cache line per element, O(n^2) in total: 0.22 -> 3.26 ms at 1024, + 0.06 -> 0.60 ms at 512, 0.016 -> 0.069 ms at 256. The pivot search is fused into the + update and is not where the time goes. +- These stage timers ran before the process gate below existed (load 5-7, no competing + build checked); they compare stages inside one call, not cells across processes. +- **Tall**: three things at once. The scaling again (4096x256: 0.39 -> 5.04 ms), the + long swap turning strided (1.25 -> 7.63 ms), and the trailing block turning wide + (update 36.1 -> 40.2 ms), the last one being the issue comment's "wide trailing block" + effect. +- So a dual that scaled the pivot row contiguously (fused into the update's column + walk, or staged through a contiguous buffer) would plausibly close the square gap; + that is a kernel change for a later slice and is **not measured**. + +## Verification conclusions and constraints + +- Tests: `faer/tests/full_piv_lu_row_major_no_copy.rs` records the largest single + allocation of the calling thread during the call; on 64, 96 and 127 square it stays + below `m * n * 8`, and on 128x128 and 256x32 (copy kept) it reaches it, so the + counter is shown to see the copy. It is an integration test because the lib's test + binary already has a `#[global_allocator]` (`sparse::linalg::amd`). With the + square exclusion removed from the guard the test fails at `64x64` (an allocation of + 32768 bytes, the whole matrix). The full-pivoting arm of + `test_row_major_adaptation_matches_the_column_major_one` now asks the new guard + whether the two orders must agree bit for bit (copy) or up to rounding (dual), and + covers 127x127. +- Measurement: tensor4all/faer#314 comment (fbench `01a1cb8`, fork row only, fork + `b16d6b8` plus a local switch that bypasses the copy, f64, `--reps 5`, three + repetitions per cell, the arms run back to back within each repetition and the cells + cycled per repetition rather than grouped by arm). +- Gate: a 1-minute load below 4 is not reachable on this host (its idle baseline is + about 6 from WindowServer and sysmond), so a cell was kept only if no competing + `rustc`, busy `cargo`, or other `target/` or bench binary was running before and + after it. Two of 44 cells were discarded and rerun. Load averages recorded over the + kept cells (1 / 5 / 15 min): 7.58-49.01 / 29.23-54.65 / 26.63-32.71, decaying from + two other agents' builds just before; the per-cell values are in the issue comment. +- macOS cannot pin threads, so the 4T cells mix core types: the two column-major + controls of a 4T cell differ by up to 15% (1024x64), against at most 1% at 1T. The 4T + ratios on large shapes are the least reliable numbers here. +- Upstream faer was not measured. diff --git a/faer/src/linalg/lu/full_pivoting/factor.rs b/faer/src/linalg/lu/full_pivoting/factor.rs index 02b4ed49..40914947 100644 --- a/faer/src/linalg/lu/full_pivoting/factor.rs +++ b/faer/src/linalg/lu/full_pivoting/factor.rs @@ -571,6 +571,29 @@ pub fn lu_in_place_scratch( .array(2) .and(StackReq::new::<(usize, usize, T::Real)>(par.degree())) } +/// Whether `lu_in_place` factors a row-major `nrows x ncols` matrix through a +/// column-major copy rather than in place through the transpose dual +/// (`lu_in_place_unblocked(.., transpose = true)`). +/// +/// The copy is [`linalg::row_major_copy_pays`]'s (tall or square, at least +/// `64 * 64` elements), less the square matrices below `128`. There the dual runs +/// within 1.05x of a column-major input, so the copy -- two passes over `m * n` +/// and an `m * n` buffer -- is a pure loss: `64 x 64` f64 one thread costs 1.24x +/// a column-major input through the copy and 1.01x through the dual on Apple M5 +/// Max, 1.32x and 1.03x on Zen 5 (tensor4all/faer#314). From `128` up the dual's +/// row scaling `A[(k, j)] *= inv` walks a strided row every step, one cache line +/// per element, and the dual leaves the band: 1.09x at `512` and `1024` on four +/// threads (M5), 1.11x at `128` and 1.08-1.27x from `256` on Zen 5 (the scaling +/// alone goes from 0.22 to 3.26 ms at `1024`, one thread). On a tall matrix the dual also makes the long swap strided and the trailing block +/// wide (1.4-1.9x), and off the square the dual leaves the 1.05x band early +/// (`80 x 64`: 1.06x), so only the square shapes skip the copy. +#[inline] +pub(crate) fn row_major_takes_the_copy(nrows: usize, ncols: usize) -> bool { + nrows >= ncols + && ncols > 1 + && linalg::row_major_copy_pays(nrows, ncols) + && !(nrows == ncols && ncols < 128) +} #[derive(Copy, Clone, Debug)] pub struct FullPivLuInfo { pub transposition_count: usize, @@ -598,13 +621,11 @@ pub fn lu_in_place<'out, I: Index, T: ComplexField>( "LU with full pivoting prefers column-major or row-major matrix. Found matrix with generic strides." ); } - // Same row-major adaptation as the partial-pivoting factorisation above. + // Same row-major adaptation as the partial-pivoting factorisation above, except + // on the small square shapes the transpose dual factors in place at the + // column-major speed. let (M, N) = mat.shape(); - if mat.row_stride().unsigned_abs() != 1 - && M >= N - && N > 1 - && linalg::row_major_copy_pays(M, N) - { + if mat.row_stride().unsigned_abs() != 1 && row_major_takes_the_copy(M, N) { let mut work = linalg::column_major_copy(mat.rb(), par.reborrow()); let mut work = MatMut::from_column_major_slice_mut(&mut work, M, N); let ret = lu_in_place( diff --git a/faer/src/linalg/mod.rs b/faer/src/linalg/mod.rs index 39d0b7bb..4d183d4e 100644 --- a/faer/src/linalg/mod.rs +++ b/faer/src/linalg/mod.rs @@ -639,12 +639,25 @@ mod tests { cm: &[T], rm: &[T], ) { - let flipped = flipped && m > 0 && n > 0; // Above the copy floor both orders take the same arithmetic and must // agree bit for bit; below it the row-major side factors in place // through a strided panel, which is the same matrix but a different // summation order. let exact = crate::linalg::row_major_copy_pays(m, n); + compare_row_major_exact_if(m, n, flipped, what, cm, rm, exact); + } + /// [`compare_row_major`] with the exactness given by the caller, for an entry + /// point whose own guard decides when the copy is taken. + fn compare_row_major_exact_if( + m: usize, + n: usize, + flipped: bool, + what: &str, + cm: &[T], + rm: &[T], + exact: bool, + ) { + let flipped = flipped && m > 0 && n > 0; for j in 0..n { for i in 0..m { let k = if flipped { m - 1 - i } else { i }; @@ -825,7 +838,18 @@ mod tests { &mut ci_rm, rm_strides(m, n, flipped), ); - compare_row_major(m, n, flipped, "full pivoting lu", &cm, &rm); + // bit for bit where the copy is taken; the square shapes below 128 go + // through the transpose dual instead (the same pivots, the complex + // products in another order) + compare_row_major_exact_if( + m, + n, + flipped, + "full pivoting lu", + &cm, + &rm, + flu::row_major_takes_the_copy(m, n), + ); same_index("full pivoting lu row perm", &p_cm, &p_rm); same_index("full pivoting lu row perm inv", &pi_cm, &pi_rm); same_index("full pivoting lu col perm", &c_cm, &c_rm); @@ -862,6 +886,7 @@ mod tests { (16, 16), (17, 15), (64, 64), + (127, 127), (128, 128), (300, 17), ] { diff --git a/faer/tests/full_piv_lu_row_major_no_copy.rs b/faer/tests/full_piv_lu_row_major_no_copy.rs new file mode 100644 index 00000000..7bea51d9 --- /dev/null +++ b/faer/tests/full_piv_lu_row_major_no_copy.rs @@ -0,0 +1,112 @@ +//! tensor4all/faer#314: a small square row-major matrix is factored by the full-pivoting LU in +//! place, through the transpose dual, and not through a column-major copy of the whole matrix. +//! The allocator records the largest single allocation the calling thread makes while +//! `lu_in_place` runs; the copy is one `m * n` buffer, so the dual path stays below that size. +//! The shapes that still take the copy are the control: there the same counter must see it. +#![cfg(all(feature = "std", feature = "threadteam"))] +use faer::dyn_stack::{MemBuffer, MemStack}; +use faer::linalg::lu::full_pivoting::factor as flu; +use faer::{MatMut, Par}; +use std::alloc::{GlobalAlloc, Layout, System}; +use std::cell::Cell; + +struct LargestOnThisThread; +thread_local! { + static COUNTING: Cell = const { Cell::new(false) }; + static LARGEST: Cell = const { Cell::new(0) }; +} +fn record(size: usize) { + // `try_with`: the allocator also runs while a thread's locals are torn down + let _ = COUNTING.try_with(|on| { + if on.get() { + let _ = LARGEST.try_with(|l| l.set(Ord::max(l.get(), size))); + } + }); +} +unsafe impl GlobalAlloc for LargestOnThisThread { + unsafe fn alloc(&self, layout: Layout) -> *mut u8 { + record(layout.size()); + unsafe { System.alloc(layout) } + } + unsafe fn alloc_zeroed(&self, layout: Layout) -> *mut u8 { + record(layout.size()); + unsafe { System.alloc_zeroed(layout) } + } + unsafe fn realloc( + &self, + ptr: *mut u8, + layout: Layout, + new_size: usize, + ) -> *mut u8 { + record(new_size); + unsafe { System.realloc(ptr, layout, new_size) } + } + unsafe fn dealloc(&self, ptr: *mut u8, layout: Layout) { + unsafe { System.dealloc(ptr, layout) } + } +} +#[global_allocator] +static GLOBAL: LargestOnThisThread = LargestOnThisThread; + +/// The largest single allocation, in bytes, the calling thread makes while a row-major `m x n` +/// f64 matrix is factored (permutations and scratch are allocated before). +fn largest_allocation(m: usize, n: usize, mut par: Par<'_>) -> usize { + let mut a: Vec = (0..m * n) + .map(|k| { + ((k as u64).wrapping_mul(0x9E37_79B9_7F4A_7C15) >> 11) as f64 + / (1u64 << 53) as f64 + - 0.5 + }) + .collect(); + let (mut p, mut pi, mut q, mut qi) = ( + vec![0usize; m], + vec![0usize; m], + vec![0usize; n], + vec![0usize; n], + ); + let mut mem = MemBuffer::new(flu::lu_in_place_scratch::( + m, + n, + par.shape(), + Default::default(), + )); + LARGEST.with(|l| l.set(0)); + COUNTING.with(|c| c.set(true)); + flu::lu_in_place( + MatMut::from_row_major_slice_mut(&mut a, m, n), + &mut p, + &mut pi, + &mut q, + &mut qi, + par.reborrow(), + MemStack::new(&mut mem), + Default::default(), + ); + COUNTING.with(|c| c.set(false)); + LARGEST.with(|l| l.get()) +} + +#[test] +fn small_square_row_major_full_piv_lu_does_not_copy_the_matrix() { + let run = |mut par: Par<'_>| { + for n in [64, 96, 127] { + let whole = n * n * size_of::(); + let largest = largest_allocation(n, n, par.reborrow()); + assert!( + largest < whole, + "{n}x{n}: an allocation of {largest} bytes, the whole matrix is {whole}" + ); + } + // the control: these shapes keep the copy, and the counter sees it + for (m, n) in [(128, 128), (256, 32)] { + let whole = m * n * size_of::(); + let largest = largest_allocation(m, n, par.reborrow()); + assert!( + largest >= whole, + "{m}x{n}: largest allocation {largest} bytes, the copy is {whole}" + ); + } + }; + run(Par::Seq); + faer::team::with_team(4, |par| run(par)); +} From 32c615f4463a20d6d523d455b2c6a10cf1079aa5 Mon Sep 17 00:00:00 2001 From: Hiroshi Shinaoka Date: Mon, 9 Nov 2026 18:28:35 +0900 Subject: [PATCH 2/6] fix(lu): keep the copy for generic strides; document the layout-dependent tie order (#314) Review of #342. The square n < 128 exclusion from the row-major copy also applied to views with neither stride 1, e.g. 96x96 with strides (2, 192): those skipped the copy and fell to the separate matrix update and pivot search instead of the fused SIMD kernel. `row_major_takes_the_copy` now takes the strides and skips the copy only when the rows are contiguous (col_stride == 1); every other stride keeps the previous dispatch. The transpose dual breaks a magnitude tie by the smallest row first, the column-major search by the smallest column, so on ties the pivots and the transposition count differ by layout (zero 64x64 with A[(1,2)] = A[(0,4)], max_rank 1: (1,2) with two transpositions vs (0,4) with one). `lu_in_place` documents it under "Pivot ties" and a test pins the case; the adaptation test no longer says "the same pivots" unconditionally. Tests: the allocation test adds 96x96 with strides (2,192) and (192,2) (copy kept; fails at 96 bytes with the old guard) and a padded row-major view (100,1) (no copy). Co-Authored-By: Claude Opus 5.5 --- faer/src/linalg/lu/full_pivoting/factor.rs | 73 ++++++++--- faer/src/linalg/mod.rs | 13 +- faer/tests/full_piv_lu_row_major_no_copy.rs | 135 ++++++++++++++++++-- 3 files changed, 185 insertions(+), 36 deletions(-) diff --git a/faer/src/linalg/lu/full_pivoting/factor.rs b/faer/src/linalg/lu/full_pivoting/factor.rs index 40914947..b98fde67 100644 --- a/faer/src/linalg/lu/full_pivoting/factor.rs +++ b/faer/src/linalg/lu/full_pivoting/factor.rs @@ -571,28 +571,43 @@ pub fn lu_in_place_scratch( .array(2) .and(StackReq::new::<(usize, usize, T::Real)>(par.degree())) } -/// Whether `lu_in_place` factors a row-major `nrows x ncols` matrix through a -/// column-major copy rather than in place through the transpose dual -/// (`lu_in_place_unblocked(.., transpose = true)`). +/// Whether `lu_in_place` factors a matrix with these strides through a +/// column-major copy rather than in place. /// -/// The copy is [`linalg::row_major_copy_pays`]'s (tall or square, at least -/// `64 * 64` elements), less the square matrices below `128`. There the dual runs -/// within 1.05x of a column-major input, so the copy -- two passes over `m * n` -/// and an `m * n` buffer -- is a pure loss: `64 x 64` f64 one thread costs 1.24x -/// a column-major input through the copy and 1.01x through the dual on Apple M5 -/// Max, 1.32x and 1.03x on Zen 5 (tensor4all/faer#314). From `128` up the dual's -/// row scaling `A[(k, j)] *= inv` walks a strided row every step, one cache line -/// per element, and the dual leaves the band: 1.09x at `512` and `1024` on four +/// The copy is [`linalg::row_major_copy_pays`]'s (no unit row stride, tall or +/// square, at least `64 * 64` elements), less one layout: a square matrix below +/// `128` whose rows are contiguous (`col_stride == 1`, any row stride, so a +/// padded or reversed row-major view counts). That one goes through the +/// transpose dual (`lu_in_place_unblocked(.., transpose = true)`), whose fused +/// SIMD update needs exactly that unit column stride. A generic-stride view +/// (neither stride `1`, or a column stride of `-1`) keeps the copy: in place it +/// would fall to the separate matrix update and pivot search. +/// +/// The square exclusion is what was measured (tensor4all/faer#314): `64 x 64` +/// f64 one thread costs 1.24x a column-major input through the copy and 1.01x +/// through the dual on Apple M5 Max, 1.32x and 1.03x on Zen 5. The dual was +/// timed at `64`, `80`, `96`, `112`, `120` and `127` on M5 and at `64` and `96` +/// on Zen 5; the sizes in between are taken on the strength of their +/// neighbours. From `128` up the dual's row scaling `A[(k, j)] *= inv` walks a +/// strided row every step, one cache line per element, and the dual leaves the +/// 1.05x band on at least one machine: 1.09x at `512` and `1024` on four /// threads (M5), 1.11x at `128` and 1.08-1.27x from `256` on Zen 5 (the scaling -/// alone goes from 0.22 to 3.26 ms at `1024`, one thread). On a tall matrix the dual also makes the long swap strided and the trailing block -/// wide (1.4-1.9x), and off the square the dual leaves the 1.05x band early -/// (`80 x 64`: 1.06x), so only the square shapes skip the copy. +/// alone goes from 0.22 to 3.26 ms at `1024`, one thread). On a tall matrix the +/// dual also makes the long swap strided and the trailing block wide +/// (1.4-2.1x), and off the square it leaves the band early (`80 x 64`: 1.06x), +/// so only the square shapes skip the copy. #[inline] -pub(crate) fn row_major_takes_the_copy(nrows: usize, ncols: usize) -> bool { - nrows >= ncols +pub(crate) fn row_major_takes_the_copy( + nrows: usize, + ncols: usize, + row_stride: isize, + col_stride: isize, +) -> bool { + row_stride.unsigned_abs() != 1 + && nrows >= ncols && ncols > 1 && linalg::row_major_copy_pays(nrows, ncols) - && !(nrows == ncols && ncols < 128) + && !(col_stride == 1 && nrows == ncols && ncols < 128) } #[derive(Copy, Clone, Debug)] pub struct FullPivLuInfo { @@ -601,6 +616,24 @@ pub struct FullPivLuInfo { /// [rank truncation](FullPivLuParams#rank-truncation) option stopped the factorization pub rank: usize, } +/// # Pivot ties +/// +/// When several entries share the largest magnitude, the pivot depends on the +/// layout. The search scans the matrix it factors column by column and keeps +/// the first strict maximum: the smallest column, then the smallest row. A +/// row-major square matrix below `128` with contiguous rows is factored through +/// its transpose (the condition is `row_major_takes_the_copy` in this module's +/// source), so there the tie goes to the smallest *row*, then the smallest +/// column, of the original matrix. Every other input -- column-major, generic +/// strides, or a row-major matrix that is copied to column-major first -- uses +/// the column-then-row order. +/// +/// For example, a zero `64 x 64` f64 matrix with `A[(1, 2)] = A[(0, 4)] = 2` +/// and `max_rank = Some(1)` pivots on `(1, 2)` (two transpositions) when +/// column-major and on `(0, 4)` (one transposition) when row-major. Both are +/// valid full-pivoting factorizations of the same matrix; the permutations, the +/// transposition count and so the factors differ between the two layouts only +/// on such ties. pub fn lu_in_place<'out, I: Index, T: ComplexField>( mat: MatMut<'_, T>, row_perm: &'out mut [I], @@ -622,10 +655,10 @@ pub fn lu_in_place<'out, I: Index, T: ComplexField>( ); } // Same row-major adaptation as the partial-pivoting factorisation above, except - // on the small square shapes the transpose dual factors in place at the - // column-major speed. + // on the small square shapes with contiguous rows the transpose dual factors + // in place at the column-major speed. let (M, N) = mat.shape(); - if mat.row_stride().unsigned_abs() != 1 && row_major_takes_the_copy(M, N) { + if row_major_takes_the_copy(M, N, mat.row_stride(), mat.col_stride()) { let mut work = linalg::column_major_copy(mat.rb(), par.reborrow()); let mut work = MatMut::from_column_major_slice_mut(&mut work, M, N); let ret = lu_in_place( diff --git a/faer/src/linalg/mod.rs b/faer/src/linalg/mod.rs index 4d183d4e..2b153f2e 100644 --- a/faer/src/linalg/mod.rs +++ b/faer/src/linalg/mod.rs @@ -839,8 +839,12 @@ mod tests { rm_strides(m, n, flipped), ); // bit for bit where the copy is taken; the square shapes below 128 go - // through the transpose dual instead (the same pivots, the complex - // products in another order) + // through the transpose dual instead, the complex products in another + // order. The dual breaks a magnitude tie by the smallest row first, the + // column-major search by the smallest column (`lu_in_place`, "Pivot + // ties"), so the pivots agree here only because this fixture does not + // hit such a tie; `full_piv_lu_row_major_no_copy` pins one where the + // two orders part. compare_row_major_exact_if( m, n, @@ -848,7 +852,10 @@ mod tests { "full pivoting lu", &cm, &rm, - flu::row_major_takes_the_copy(m, n), + { + let (_, rs, cs) = rm_strides(m, n, flipped); + flu::row_major_takes_the_copy(m, n, rs, cs) + }, ); same_index("full pivoting lu row perm", &p_cm, &p_rm); same_index("full pivoting lu row perm inv", &pi_cm, &pi_rm); diff --git a/faer/tests/full_piv_lu_row_major_no_copy.rs b/faer/tests/full_piv_lu_row_major_no_copy.rs index 7bea51d9..57a9c0ae 100644 --- a/faer/tests/full_piv_lu_row_major_no_copy.rs +++ b/faer/tests/full_piv_lu_row_major_no_copy.rs @@ -1,12 +1,14 @@ -//! tensor4all/faer#314: a small square row-major matrix is factored by the full-pivoting LU in -//! place, through the transpose dual, and not through a column-major copy of the whole matrix. -//! The allocator records the largest single allocation the calling thread makes while -//! `lu_in_place` runs; the copy is one `m * n` buffer, so the dual path stays below that size. -//! The shapes that still take the copy are the control: there the same counter must see it. +//! tensor4all/faer#314: a small square row-major matrix is factored by the +//! full-pivoting LU in place, through the transpose dual, and not through a +//! column-major copy of the whole matrix. The allocator records the largest +//! single allocation the calling thread makes while `lu_in_place` runs; the +//! copy is one `m * n` buffer, so the dual path stays below that size. +//! The shapes that still take the copy are the control: there the same counter +//! must see it. #![cfg(all(feature = "std", feature = "threadteam"))] use faer::dyn_stack::{MemBuffer, MemStack}; use faer::linalg::lu::full_pivoting::factor as flu; -use faer::{MatMut, Par}; +use faer::{Auto, MatMut, Par}; use std::alloc::{GlobalAlloc, Layout, System}; use std::cell::Cell; @@ -28,10 +30,12 @@ unsafe impl GlobalAlloc for LargestOnThisThread { record(layout.size()); unsafe { System.alloc(layout) } } + unsafe fn alloc_zeroed(&self, layout: Layout) -> *mut u8 { record(layout.size()); unsafe { System.alloc_zeroed(layout) } } + unsafe fn realloc( &self, ptr: *mut u8, @@ -41,6 +45,7 @@ unsafe impl GlobalAlloc for LargestOnThisThread { record(new_size); unsafe { System.realloc(ptr, layout, new_size) } } + unsafe fn dealloc(&self, ptr: *mut u8, layout: Layout) { unsafe { System.dealloc(ptr, layout) } } @@ -48,10 +53,17 @@ unsafe impl GlobalAlloc for LargestOnThisThread { #[global_allocator] static GLOBAL: LargestOnThisThread = LargestOnThisThread; -/// The largest single allocation, in bytes, the calling thread makes while a row-major `m x n` -/// f64 matrix is factored (permutations and scratch are allocated before). -fn largest_allocation(m: usize, n: usize, mut par: Par<'_>) -> usize { - let mut a: Vec = (0..m * n) +/// The largest single allocation, in bytes, the calling thread makes while an +/// `m x n` f64 matrix with element strides `(rs, cs)` (both positive) is +/// factored (permutations and scratch are allocated before). +fn largest_allocation_strided( + m: usize, + n: usize, + (rs, cs): (usize, usize), + mut par: Par<'_>, +) -> usize { + let len = (m - 1) * rs + (n - 1) * cs + 1; + let mut a: Vec = (0..len) .map(|k| { ((k as u64).wrapping_mul(0x9E37_79B9_7F4A_7C15) >> 11) as f64 / (1u64 << 53) as f64 @@ -70,10 +82,22 @@ fn largest_allocation(m: usize, n: usize, mut par: Par<'_>) -> usize { par.shape(), Default::default(), )); + // SAFETY: `a` holds `(m - 1) * rs + (n - 1) * cs + 1` elements, and the + // strides address distinct elements (`rs >= n * cs` or `cs >= m * rs` in + // every caller) + let view = unsafe { + MatMut::from_raw_parts_mut( + a.as_mut_ptr(), + m, + n, + rs as isize, + cs as isize, + ) + }; LARGEST.with(|l| l.set(0)); COUNTING.with(|c| c.set(true)); flu::lu_in_place( - MatMut::from_row_major_slice_mut(&mut a, m, n), + view, &mut p, &mut pi, &mut q, @@ -85,6 +109,10 @@ fn largest_allocation(m: usize, n: usize, mut par: Par<'_>) -> usize { COUNTING.with(|c| c.set(false)); LARGEST.with(|l| l.get()) } +/// [`largest_allocation_strided`] for a contiguous row-major matrix. +fn largest_allocation(m: usize, n: usize, par: Par<'_>) -> usize { + largest_allocation_strided(m, n, (n, 1), par) +} #[test] fn small_square_row_major_full_piv_lu_does_not_copy_the_matrix() { @@ -94,19 +122,100 @@ fn small_square_row_major_full_piv_lu_does_not_copy_the_matrix() { let largest = largest_allocation(n, n, par.reborrow()); assert!( largest < whole, - "{n}x{n}: an allocation of {largest} bytes, the whole matrix is {whole}" + "{n}x{n}: an allocation of {largest} bytes, the whole matrix \ + is {whole}" ); } + // contiguous rows with a padded leading dimension are still row-major + let whole = 96 * 96 * size_of::(); + let largest = + largest_allocation_strided(96, 96, (100, 1), par.reborrow()); + assert!( + largest < whole, + "96x96, strides (100, 1): an allocation of {largest} bytes, the \ + whole matrix is {whole}" + ); // the control: these shapes keep the copy, and the counter sees it for (m, n) in [(128, 128), (256, 32)] { let whole = m * n * size_of::(); let largest = largest_allocation(m, n, par.reborrow()); assert!( largest >= whole, - "{m}x{n}: largest allocation {largest} bytes, the copy is {whole}" + "{m}x{n}: largest allocation {largest} bytes, the copy is \ + {whole}" + ); + } + // a generic-stride view (neither stride 1) is not row-major: it keeps + // the copy (tensor4all/faer#342 review), where in place it would miss + // the fused SIMD update + for (rs, cs) in [(2, 192), (96 * 2, 2)] { + let largest = + largest_allocation_strided(96, 96, (rs, cs), par.reborrow()); + assert!( + largest >= whole, + "96x96, strides ({rs}, {cs}): largest allocation {largest} \ + bytes, the copy is {whole}" ); } }; run(Par::Seq); faer::team::with_team(4, |par| run(par)); } + +/// The pivot search keeps the first strict maximum of the matrix it factors: +/// the smallest column, then the smallest row. The transpose dual factors the +/// transpose, so a row-major input breaks the same tie by the smallest row +/// first (`lu_in_place`, "Pivot ties"). +#[test] +fn full_piv_lu_tie_order_depends_on_the_layout() { + let n = 64; + let at = |i: usize, j: usize| { + if (i, j) == (1, 2) || (i, j) == (0, 4) { + 2.0 + } else { + 0.0 + } + }; + let first_pivot = |row_major: bool| { + let mut a = vec![0.0f64; n * n]; + for i in 0..n { + for j in 0..n { + a[if row_major { i * n + j } else { i + j * n }] = at(i, j); + } + } + let view = if row_major { + MatMut::from_row_major_slice_mut(&mut a, n, n) + } else { + MatMut::from_column_major_slice_mut(&mut a, n, n) + }; + let (mut p, mut pi, mut q, mut qi) = ( + vec![0usize; n], + vec![0usize; n], + vec![0usize; n], + vec![0usize; n], + ); + let mut mem = MemBuffer::new(flu::lu_in_place_scratch::( + n, + n, + Par::Seq.shape(), + Default::default(), + )); + let params = flu::FullPivLuParams { + max_rank: Some(1), + ..>::auto() + }; + let (info, _, _) = flu::lu_in_place( + view, + &mut p, + &mut pi, + &mut q, + &mut qi, + Par::Seq, + MemStack::new(&mut mem), + params.into(), + ); + (p[0], q[0], info.transposition_count) + }; + assert_eq!(first_pivot(false), (1, 2, 2)); + assert_eq!(first_pivot(true), (0, 4, 1)); +} From 17270db2308c944e4aa82fa3e2b223a7331d0dbe Mon Sep 17 00:00:00 2001 From: Hiroshi Shinaoka Date: Mon, 9 Nov 2026 19:03:09 +0900 Subject: [PATCH 3/6] docs(lu): limit the #314 boundary claims to the measured sizes; record the clean-head run The worklog and PATCHSET said square 64..=127 passes on both machines; the timed sizes are 64, 80, 96, 112, 120 and 127 on Apple M5 Max and 64, 96 on Zen 5, the rest is interpolated. The earlier run measured b16d6b8 plus an uncommitted copy-bypass switch; the record now adds a paired run of the committed 32c615f source against b16d6b8 (f64/c64, 1T/4T, three repetitions: dual 0.89-1.02x a column-major input on the 24 changed cells) and names the switch run as not reproducible from a revision. The lu_in_place tie note no longer says the factors differ only on ties. Co-Authored-By: Claude Opus 5.5 --- PATCHSET.md | 2 +- .../2026-11-09-full-piv-lu-row-major-dual.md | 89 +++++++++++++++---- faer/src/linalg/lu/full_pivoting/factor.rs | 6 +- 3 files changed, 74 insertions(+), 23 deletions(-) diff --git a/PATCHSET.md b/PATCHSET.md index 43b4a817..8c15e512 100644 --- a/PATCHSET.md +++ b/PATCHSET.md @@ -121,7 +121,7 @@ gate has changed and two kinds of adoption are now in the ledger. | Empty complex tridiagonal EVD with eigenvectors | `tridiagonal_self_adjoint_evd` had no empty-input return, unlike `self_adjoint_evd`. For a complex type with eigenvectors requested, its phase loop wrote `scale[0]` into an empty buffer and panicked; real types and the values-only path did not reach that write. It now returns `Ok(())` for `n == 0`. The function is upstream faer code (sarah-quiñones); nothing was reported upstream. | KEEP | tensor4all/faer#31 | Correctness only. Found by the #31 LAPACK validation harness (`lapack_tridiagonal_evd`, 36 panicking calls per complex type). Test `empty_tridiagonal_with_vectors` (f64/c64/c32) fails before the change and passes after. Upstream report: [`upstream-reports/faer-021-tridiag-evd-empty-complex`](upstream-reports/faer-021-tridiag-evd-empty-complex/). | | Row-major QR layout adaptation | `qr_in_place` factors a matrix whose columns are not the contiguous direction through a column-major copy: the panel is not layout-aware (`make_householder_in_place` reads the pivot column and scales it, and the rank-one update beside it walked the block with `z!`), so a row-major matrix cost one cache line per element -- at `4096x256` f64 the panel was 25x its column-major cost and the whole deficit, and the factorization did not scale (47.8 / 42.9 / 43.0 ms at 1T / 4T / 8T). `qr_in_place` now copies into a column-major buffer, factors that, and copies the factor back, which is what LAPACK's row-major interface pays internally; the arithmetic is the copy's, so the factor is the column-major factor of the same matrix bit for bit. Tiled copy; with threads it is split along the destination's major dimension in whole 16-element tiles (columns into the buffer, rows on the write-back, so a thread writes whole lines of its own: tensor4all/faer#92), and the work buffer is allocated uninitialised and initialised by the threads that copy into it, so its pages are first touched in parallel rather than by a serial `vec![zero; m * n]` on the caller (`linalg::column_major_copy`, shared by the four row-major adaptations). One `m * n` heap buffer per call, because `qr_in_place_scratch` cannot know the strides. The copy has no cost gate yet (tensor4all/faer#92 part 2, after #91's dispatch measurements). Paired, best of 21, CPUs 4 / 4-7 / 4-11, fork `bfe4041` -> `044f947`, C order ms (before -> after, F order after): `4096x256` f64 47.83 / 42.85 / 42.96 -> 19.44 / 10.22 / 8.39 (15.77 / 8.37 / 7.23), `1024x64` f64 1.090 / 2.687 / 2.919 -> 0.587 / 0.613 / 0.678, `512x512` f64 9.44 / 9.56 / 12.73 -> 5.24 / 3.11 / 3.24, `256x256` f64 1.455 / 1.596 / 2.315 -> 0.921 / 0.833 / 0.904. Against the campaign's MKL C-order column (27.73 / 8.66 / 8.20) `4096x256` f64 goes from 1.72x / 4.95x / 5.24x behind to 0.70x / 1.18x / 1.02x. Test: `test_row_major_factor_is_the_column_major_one` -- the factor, `tau` and rank bit for bit at f64 and c64, tall/wide/square, sequential and threaded, backwards strides, and a zero leading column. A one-row matrix is left on the direct path (it eliminates at most one column). The `perf-warn` that told row-major callers to transpose is gone. | KEEP | tlinalg-rs#46 | Heap allocation inside a public call, the same boundary the fork's two row splits already cross. Own code: no OpenBLAS or MKL source is used or adapted. Exactly the adapted shapes (tall or square) are the column-major factorisation bit for bit; the guards' exclusions are not. Why, and what lost: `docs/worklogs/2026-10-15-qr-row-major-layout.md`. | | Row-major col-pivoting QR | `qr::col_pivoting::factor::qr_in_place` computes the norm of every one of the `n` columns before it starts and re-reduces them per step, and its incremental-norm SIMD path is gated on `A.row_stride() == 1` and unwraps a column-major view, so a row-major matrix pays one cache line per element on every shape. It now factors a column-major copy and writes the factor back, wide included. Paired, best of 21, CPUs 4 / 4-7 / 4-11, `f64`, C order before -> after (F order after): `256x256` 58.96 / 59.63 / 59.61 -> 1.52 / 1.40 / 1.19 (1.45 / 1.33 / 1.54), `1024x64` 22.23 / 22.29 / 22.45 -> 0.726 / 0.520 / 0.509 (0.472 / 0.522 / 0.536); `c64` `256x256` one thread 64.40 -> 3.29, `1024x64` 24.69 -> 1.64. So the C/F ratio at `256x256` one thread is 38.8x -> 1.05 and the C order is *faster* than the F order at eight threads. Nothing in the campaign records a column-pivoted QR, so this is the record's own driver. Test: `linalg::tests::test_row_major_adaptation_matches_the_column_major_one`. | KEEP | tlinalg-rs#46 | Heap allocation inside a public call, as for the no-pivoting row. Why, and what lost: `docs/worklogs/2026-10-15-row-major-adaptation-lu-colpiv.md`. | -| Row-major LU, partial pivoting | `lu::partial_pivoting::factor::lu_in_place` walks each pivot column twice per step (search, then scaling) and swaps elements *inside* a column, which is a cache line per element when the columns are not contiguous. It now factors a column-major copy and writes the factor and its permutations back. Paired, best of 21, `f64`, C order before -> after (F order after): `512x512` 12.14 / 7.02 / 6.34 -> 3.37 / 2.01 / 1.89 (3.04 / 1.88 / 1.81) at 1T / 4T / 8T, `1024x1024` 59.80 / 27.62 / 22.93 -> 22.91 / 10.04 / 8.21 (19.89 / 8.19 / 6.65), `256x256` one thread 1.530 -> 0.486; `c64` `512x512` one thread 18.08 -> 9.50. Test: `linalg::tests::test_row_major_adaptation_matches_the_column_major_one` -- factor, permutations and transposition count bit for bit. The campaign's `lu` cell takes the composite, so `to_owned` below is what moves it. | KEEP | tlinalg-rs#46 | The full-pivoting flavour is deliberately **not** adapted: it already factors a row-major input through `mat.transpose_mut()`, and the copy measured a no-op there (512x512 f64 one thread 13.758 against 13.712 C order); `10ee1c9` adapted it after all, behind the same `row_major_copy_pays` floor. Since tensor4all/faer#314 the square row-major matrices below 128 skip that copy and go through the transpose dual, which runs within 1.05x of a column-major input there (64x64 f64 one thread: 1.24x through the copy, 1.01x through the dual on Apple M5 Max; 1.32x and 1.03x on Zen 5); from 128 up and on every tall shape the copy stays, because the dual's strided row scaling (and, tall, its strided long swap and wide trailing block) costs more than the copy. Guard `full_pivoting::factor::row_major_takes_the_copy`; tests `faer/tests/full_piv_lu_row_major_no_copy.rs` (no `m * n` allocation on the dual shapes, the copy still seen on 128x128 and 256x32) and the full-pivoting arm of `test_row_major_adaptation_matches_the_column_major_one` (127x127 added). Why: `docs/worklogs/2026-11-09-full-piv-lu-row-major-dual.md`. Why: `docs/worklogs/2026-10-15-row-major-adaptation-lu-colpiv.md`. | +| Row-major LU, partial pivoting | `lu::partial_pivoting::factor::lu_in_place` walks each pivot column twice per step (search, then scaling) and swaps elements *inside* a column, which is a cache line per element when the columns are not contiguous. It now factors a column-major copy and writes the factor and its permutations back. Paired, best of 21, `f64`, C order before -> after (F order after): `512x512` 12.14 / 7.02 / 6.34 -> 3.37 / 2.01 / 1.89 (3.04 / 1.88 / 1.81) at 1T / 4T / 8T, `1024x1024` 59.80 / 27.62 / 22.93 -> 22.91 / 10.04 / 8.21 (19.89 / 8.19 / 6.65), `256x256` one thread 1.530 -> 0.486; `c64` `512x512` one thread 18.08 -> 9.50. Test: `linalg::tests::test_row_major_adaptation_matches_the_column_major_one` -- factor, permutations and transposition count bit for bit. The campaign's `lu` cell takes the composite, so `to_owned` below is what moves it. | KEEP | tlinalg-rs#46 | The full-pivoting flavour is deliberately **not** adapted: it already factors a row-major input through `mat.transpose_mut()`, and the copy measured a no-op there (512x512 f64 one thread 13.758 against 13.712 C order); `10ee1c9` adapted it after all, behind the same `row_major_copy_pays` floor. Since tensor4all/faer#314 the square row-major matrices below 128 with contiguous rows (`col_stride == 1`) skip that copy and go through the transpose dual; generic strides keep the copy. Timed at 64, 80, 96, 112, 120 and 127 on Apple M5 Max (f64/c64, 1T/4T, `32c615f` against `b16d6b8`, per-repetition median of three: dual 0.89-1.02x a column-major input, the copy 0.98-1.34x) and at 64 and 96 on Zen 5 (0.94-1.04x); the sizes in between were not timed. From 128 up and on every tall shape the copy stays, because the dual's strided row scaling (and, tall, its strided long swap and wide trailing block) costs more than the copy. On a magnitude tie the dual picks the smallest row first, the column-major search the smallest column, so tied pivots differ by layout (documented on `lu_in_place`, "Pivot ties"). Guard `full_pivoting::factor::row_major_takes_the_copy`; tests `faer/tests/full_piv_lu_row_major_no_copy.rs` (no `m * n` allocation on the dual shapes, the copy still seen on 128x128, 256x32 and generic-stride 96x96; the tie case) and the full-pivoting arm of `test_row_major_adaptation_matches_the_column_major_one` (127x127 added). Why: `docs/worklogs/2026-11-09-full-piv-lu-row-major-dual.md`. Why: `docs/worklogs/2026-10-15-row-major-adaptation-lu-colpiv.md`. | | Row-major `to_owned` | `MatRef::to_owned` builds its result with `Mat::from_fn`, which fills column by column, so a row-major source is read with a stride of `ncols` and every element costs its own cache line: `512x512` f64 1.01 ms against 0.057, `1024x1024` 6.36 against 0.36, `4096x256` 5.08 against 0.37, `256x4096` 4.30 against 0.36. It now stages a tile when the source's columns are not contiguous and the source is already canonical. After: 0.275 / 2.007 / 2.045 / 1.774 ms. This is what the campaign's `lu` and `qr` cells pay -- a probe at `n = 512` f64 one thread C order reads `lu` 6.1417 -> 4.156 and 4.5329 -> 3.033 at eight threads, `qr` 12.5143 -> 11.185 -- and `full_piv_lu` 12.95 -> 12.63, `cholesky`, `eigh` and `svd_values` unchanged, because those three copy through their own routes rather than `to_owned`. | KEEP | tlinalg-rs#46 | The tile adds a `Mat::zeros` pass and is serial (`to_owned` has no `par`). Gated on `T::IS_CANONICAL`: `canonical()` reinterprets memory and does not apply `Conj::apply`, so a conjugated source (`A.adjoint().to_owned()`, which the sparse Cholesky tests use) keeps the elementwise path -- the ungated attempt broke five of them. The attempted generic `cloned` tile is **dropped**: initializing every destination entry by cloning the first source element amplifies heap-owned payloads by `m * n` before overwrite. `cloned` retains the original elementwise implementation; `test_cloned_does_not_replicate_the_first_heap_value` checks that each source is cloned once. The `Option` tile still moves each cloned numeric value into `to_owned`/factorization destinations without a second clone. Why: `docs/worklogs/2026-10-15-row-major-adaptation-lu-colpiv.md`. | | Row-major `copy_from_triangular_lower` | `MatMut::copy_from_triangular_lower` (what `Llt::new` and `Ldlt::new` copy their input with) ran the elementwise zip, which for a row-major source into a column-major destination reads with a stride of the source's row length: `1024x1024` f64 one thread 4.16 ms against 0.30 ms from a column-major source. It now copies each 16-column band as its diagonal tile plus the rectangle under it through `copy_across_layouts_block`, the tile `to_owned` uses: 0.65 ms. The campaign's `cholesky` cell, 1024 f64 8T, three processes per arm against MKL in the same process: row-major 5.445 -> 2.236 ms (MKL 3.33 / 3.39), column-major 1.674 -> 1.649 (MKL 2.65 / 2.67). Tests: `linalg::tests::tiled_triangular_copy_writes_what_the_elementwise_one_does` and `linalg::tests::llt_and_partial_piv_lu_factor_both_orders_bit_for_bit` (the factor of the same matrix through the two orders, bit for bit). | KEEP | faer-benchmark slide re-take (1024, 8T) | Serial, like `to_owned`: the method has no `par`. A conjugated source keeps the elementwise path. `copy_from` (the full-matrix copy) still runs elementwise: 8.6 ms row-major against 0.47 ms at `1024x1024`; it is not on the `cholesky` or `lu` path. Why: `docs/worklogs/2026-11-06-row-major-triangular-copy.md`. | | 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. | diff --git a/docs/worklogs/2026-11-09-full-piv-lu-row-major-dual.md b/docs/worklogs/2026-11-09-full-piv-lu-row-major-dual.md index 0e0f30da..21395a64 100644 --- a/docs/worklogs/2026-11-09-full-piv-lu-row-major-dual.md +++ b/docs/worklogs/2026-11-09-full-piv-lu-row-major-dual.md @@ -5,21 +5,41 @@ - `full_pivoting::factor::lu_in_place` used to take the shared row-major copy (`linalg::row_major_copy_pays`, at least `64 * 64` elements, tall or square) before its own transpose dual (`lu_in_place_unblocked(.., transpose = true)`) was even tried. - It now skips the copy for **square matrices below 128** and keeps it everywhere else - (`row_major_takes_the_copy`). + It now skips the copy for **square matrices below 128 whose rows are contiguous** + (`col_stride == 1`, any row stride, so a padded or reversed row-major view counts) and + keeps it everywhere else (`row_major_takes_the_copy`, which takes the strides). +- The stride condition came from the #342 review: the first version tested only the + shape, so a generic-stride view such as 96x96 with strides `(2, 192)` also skipped the + copy and, in place, fell to the separate matrix update and pivot search instead of + the fused SIMD kernel (that kernel needs the unit column stride the dual sees as a + unit row stride). Generic strides, and a column stride of `-1`, keep the copy as on + `main`. That path was not timed; the allocation test pins the dispatch. - The rule is the issue's acceptance criterion applied shape by shape: the dual has to - run within 1.05x of a column-major input. Square `64..=127` passes on both machines - measured (Apple M5 Max, this change: f64 C/A 1.009 / 1.039 at 64 and 1.007 / 1.005 - at 128, one / four threads; in an ungated scan, 64-120 square at 0.98-1.04 in f64 and c64; Zen 5, the - issue's earlier comment: n = 64 and 96 at 0.94-1.04 over f64/c64/f32 at 1T and 8T). + run within 1.05x of a column-major input. **Timed sizes**: on Apple M5 Max, square + 64, 80, 96, 112, 120 and 127 in f64 and c64 at one and four threads (the paired run + below: per-repetition median C/A 0.888-1.017 over those 24 cells); on Zen 5, the + issue's earlier comment, n = 64 and 96 at 0.94-1.04 over f64/c64/f32 at 1T and 8T. + The sizes in between (65-79, 81-95, ...) and every size other than 64 and 96 on + Zen 5 were **not timed**; they are taken on the strength of their neighbours. Square `128` is mixed on Zen 5 (f64 1T 1.114) and from 256 up the dual leaves the band on at least one machine (M5 four threads 1.087 at 512, 1.089 at 1024; Zen 5 1.08-1.27), so the boundary sits below 128. On M5 one thread the square dual stays at 1.03-1.04 up to 1024, so a single-thread-only rule could go further; not taken. - Off the square the band ends early (`80 x 64` f64 1.06x, `128 x 64` 1.11-1.14x on - M5, the same ungated scan), so near-square shapes keep the copy. + M5, an ungated scan run before measurement 1 at a 1-min load of about 7), so + near-square shapes keep the copy. - No change to the dual kernel in this slice, and no change to `row_major_copy_pays` or to the other families that share it. +- **Pivot ties depend on the layout.** The pivot search keeps the first strict maximum + of the matrix it factors (smallest column, then smallest row). The dual factors the + transpose, so a row-major input breaks a tie by the smallest original *row* first. + Example (`faer/tests/full_piv_lu_row_major_no_copy.rs`): zero 64x64 f64 with + `A[(1, 2)] = A[(0, 4)] = 2`, `max_rank = Some(1)`, `Par::Seq` pivots on `(1, 2)` with + two transpositions column-major and on `(0, 4)` with one row-major. #314 allows a + different pivot order on ties if it is documented; it is, under "Pivot ties" on + `lu_in_place`. On the adaptation test's fixtures, which hit no such tie, the + pivots, permutations and transposition counts agree; that is all the test shows. The + first version of this change said "the same pivots" without the tie exception. ## Where the dual loses (stage breakdown, M5 Max, f64, one thread) @@ -44,17 +64,47 @@ Per-stage timers inside `lu_in_place_unblocked` (local patch, not committed): ## Verification conclusions and constraints - Tests: `faer/tests/full_piv_lu_row_major_no_copy.rs` records the largest single - allocation of the calling thread during the call; on 64, 96 and 127 square it stays - below `m * n * 8`, and on 128x128 and 256x32 (copy kept) it reaches it, so the - counter is shown to see the copy. It is an integration test because the lib's test - binary already has a `#[global_allocator]` (`sparse::linalg::amd`). With the - square exclusion removed from the guard the test fails at `64x64` (an allocation of - 32768 bytes, the whole matrix). The full-pivoting arm of - `test_row_major_adaptation_matches_the_column_major_one` now asks the new guard - whether the two orders must agree bit for bit (copy) or up to rounding (dual), and - covers 127x127. -- Measurement: tensor4all/faer#314 comment (fbench `01a1cb8`, fork row only, fork - `b16d6b8` plus a local switch that bypasses the copy, f64, `--reps 5`, three + allocation of the calling thread during the call; on 64, 96 and 127 square, and on + a padded 96x96 row-major view (strides `(100, 1)`), it stays below `m * n * 8`; on + 128x128 and 256x32 (copy kept) and on 96x96 views with strides `(2, 192)` and + `(192, 2)` it reaches it, so the counter is shown to see the copy. It is an + integration test because the lib's test binary already has a `#[global_allocator]` + (`sparse::linalg::amd`). Red checks: with the square exclusion removed the test + fails at `64x64` (an allocation of 32768 bytes, the whole matrix); with only the + `col_stride == 1` condition removed it fails at `96x96, strides (2, 192)` (largest + allocation 96 bytes). A second test pins the tie example above. The full-pivoting + arm of `test_row_major_adaptation_matches_the_column_major_one` now asks the new + guard whether the two orders must agree bit for bit (copy) or up to rounding (dual), + and covers 127x127. +- Whole crate suite (`cargo test -p t4a-faer --no-fail-fast`, this host): 957 passed, + 1 failed, `lq_allocates_nothing_60x300` (allocation counts `[13, 12, 17, 17]`), which + fails the same way on `b16d6b8` and is one of the host failures listed in + `2026-11-04-gevd-repeated-complex-pairs.md`. +- **Measurement 2 (the one that covers the merged code)**: the source of `32c615f` + (clean committed tree; `faer/src` tree `38ce727`, unchanged by the later doc-only + commits) against `main` `b16d6b8`, both built into fbench `01a1cb8` with the same + one-line local `set_threads` patch (returns Ok for `NAME == "none"`), family + `full_piv_lu`, row `faer-fork`, the `in_place` line, `--reps 5 --prime-ms 250 + --statistic median`. Cells: square 64, 80, 96, 112, 120, 127 (changed) and 128, 256 + (unchanged code, the control), f64 and c64, 1T and 4T, three repetitions; per cell + and repetition the two binaries run back to back (head first on odd repetitions, + base first on even). Each process reports its own column-major A and row-major C (or + B), so the ratio is taken within a process and the median over the three + repetitions is reported. Gate as below (no `rustc`, busy `cargo`, or other + `target/` binary before and after); 11 cells were discarded and rerun, 96 kept, 1-min + load before a kept cell 7.0-26.4, 18:41-19:01 JST. `fbench verify` passed all 64 + cells per thread count (1T, 4T) on the head binary. The full table is in the #342 + description. Result: C/A 0.888-1.017 on all 24 changed cells, base B/A 0.978-1.343 + (median per cell); the unchanged 128/256 controls read C/A 0.976-1.181 at head and + B/A 1.005-1.150 at base. Single repetitions swing far more (one control cell 0.650, + another 1.616; the two binaries' A arms of one repetition differ by up to 1.8x), so + per-repetition medians are the only statistic here, and they are a scan, not a + campaign `candidate_compare.py` result (that was not run). +- Measurement 1 (the large-shape and boundary evidence above, 64-1024 square and the + tall shapes): tensor4all/faer#314 comment (fbench `01a1cb8`, fork row only, fork + `b16d6b8` plus a local, never committed `FAER_FLU_COPY=off` environment switch that + skipped the copy branch, i.e. the dual on every shape; that switch is in no published + revision, so arm C of that run cannot be rebuilt from a checkout. f64, `--reps 5`, three repetitions per cell, the arms run back to back within each repetition and the cells cycled per repetition rather than grouped by arm). - Gate: a 1-minute load below 4 is not reachable on this host (its idle baseline is @@ -65,5 +115,6 @@ Per-stage timers inside `lu_in_place_unblocked` (local patch, not committed): two other agents' builds just before; the per-cell values are in the issue comment. - macOS cannot pin threads, so the 4T cells mix core types: the two column-major controls of a 4T cell differ by up to 15% (1024x64), against at most 1% at 1T. The 4T - ratios on large shapes are the least reliable numbers here. + ratios on large shapes are the least reliable numbers here. Measurement 2 shows the + same at small sizes. - Upstream faer was not measured. diff --git a/faer/src/linalg/lu/full_pivoting/factor.rs b/faer/src/linalg/lu/full_pivoting/factor.rs index b98fde67..609c98d5 100644 --- a/faer/src/linalg/lu/full_pivoting/factor.rs +++ b/faer/src/linalg/lu/full_pivoting/factor.rs @@ -631,9 +631,9 @@ pub struct FullPivLuInfo { /// For example, a zero `64 x 64` f64 matrix with `A[(1, 2)] = A[(0, 4)] = 2` /// and `max_rank = Some(1)` pivots on `(1, 2)` (two transpositions) when /// column-major and on `(0, 4)` (one transposition) when row-major. Both are -/// valid full-pivoting factorizations of the same matrix; the permutations, the -/// transposition count and so the factors differ between the two layouts only -/// on such ties. +/// valid full-pivoting factorizations of the same matrix, but the +/// permutations, the transposition count and so the factors differ between the +/// two layouts. pub fn lu_in_place<'out, I: Index, T: ComplexField>( mat: MatMut<'_, T>, row_perm: &'out mut [I], From 63790b0f86d73b68113fd3118ea6ba7756dabbde Mon Sep 17 00:00:00 2001 From: Hiroshi Shinaoka Date: Mon, 9 Nov 2026 19:45:56 +0900 Subject: [PATCH 4/6] fix(lu): keep the smaller row when the SIMD pivot search merges tied accumulators (#314) Second review of #342. The SIMD full-pivot search walks each column with four interleaved accumulators (three in the fused update) and merged them with a strict >, so an equal maximum held by the accumulator merged first won even with the larger row: zero 64x64 row-major with A[(2,6)] = A[(2,8)] = 2 and max_rank 1 pivoted on (2,8) on a two-lane build. best_score now also takes the other accumulator on an equal score with a smaller row. Inherited from upstream; it changes the tied pivot for every layout, to the order the scalar search and reduce_2d's comment already describe. The lu_in_place "Pivot ties" section described the copy exclusion as the whole condition. It now states the dispatch: the copy when row_major_takes_the_copy holds, else the input when |row_stride| < |col_stride| and its transpose otherwise (an 8x16 row-major matrix is never copied and picks (0,4), column-major (1,2)). Tests: pivot_search_keeps_the_smaller_row_across_accumulators (every row pair of a 64-row column, plain and fused search, three alignments; failed at rows (0,2) before the fix) and the integration tie test (8x16, 256x32 copied, and the within-column cases; returned (2,8,2) before the fix). Co-Authored-By: Claude Opus 5.5 --- PATCHSET.md | 2 +- .../2026-11-09-full-piv-lu-row-major-dual.md | 46 +++++-- faer/src/linalg/lu/full_pivoting/factor.rs | 86 +++++++++++-- faer/tests/full_piv_lu_row_major_no_copy.rs | 120 ++++++++++-------- 4 files changed, 175 insertions(+), 79 deletions(-) diff --git a/PATCHSET.md b/PATCHSET.md index 8c15e512..0759bc19 100644 --- a/PATCHSET.md +++ b/PATCHSET.md @@ -121,7 +121,7 @@ gate has changed and two kinds of adoption are now in the ledger. | Empty complex tridiagonal EVD with eigenvectors | `tridiagonal_self_adjoint_evd` had no empty-input return, unlike `self_adjoint_evd`. For a complex type with eigenvectors requested, its phase loop wrote `scale[0]` into an empty buffer and panicked; real types and the values-only path did not reach that write. It now returns `Ok(())` for `n == 0`. The function is upstream faer code (sarah-quiñones); nothing was reported upstream. | KEEP | tensor4all/faer#31 | Correctness only. Found by the #31 LAPACK validation harness (`lapack_tridiagonal_evd`, 36 panicking calls per complex type). Test `empty_tridiagonal_with_vectors` (f64/c64/c32) fails before the change and passes after. Upstream report: [`upstream-reports/faer-021-tridiag-evd-empty-complex`](upstream-reports/faer-021-tridiag-evd-empty-complex/). | | Row-major QR layout adaptation | `qr_in_place` factors a matrix whose columns are not the contiguous direction through a column-major copy: the panel is not layout-aware (`make_householder_in_place` reads the pivot column and scales it, and the rank-one update beside it walked the block with `z!`), so a row-major matrix cost one cache line per element -- at `4096x256` f64 the panel was 25x its column-major cost and the whole deficit, and the factorization did not scale (47.8 / 42.9 / 43.0 ms at 1T / 4T / 8T). `qr_in_place` now copies into a column-major buffer, factors that, and copies the factor back, which is what LAPACK's row-major interface pays internally; the arithmetic is the copy's, so the factor is the column-major factor of the same matrix bit for bit. Tiled copy; with threads it is split along the destination's major dimension in whole 16-element tiles (columns into the buffer, rows on the write-back, so a thread writes whole lines of its own: tensor4all/faer#92), and the work buffer is allocated uninitialised and initialised by the threads that copy into it, so its pages are first touched in parallel rather than by a serial `vec![zero; m * n]` on the caller (`linalg::column_major_copy`, shared by the four row-major adaptations). One `m * n` heap buffer per call, because `qr_in_place_scratch` cannot know the strides. The copy has no cost gate yet (tensor4all/faer#92 part 2, after #91's dispatch measurements). Paired, best of 21, CPUs 4 / 4-7 / 4-11, fork `bfe4041` -> `044f947`, C order ms (before -> after, F order after): `4096x256` f64 47.83 / 42.85 / 42.96 -> 19.44 / 10.22 / 8.39 (15.77 / 8.37 / 7.23), `1024x64` f64 1.090 / 2.687 / 2.919 -> 0.587 / 0.613 / 0.678, `512x512` f64 9.44 / 9.56 / 12.73 -> 5.24 / 3.11 / 3.24, `256x256` f64 1.455 / 1.596 / 2.315 -> 0.921 / 0.833 / 0.904. Against the campaign's MKL C-order column (27.73 / 8.66 / 8.20) `4096x256` f64 goes from 1.72x / 4.95x / 5.24x behind to 0.70x / 1.18x / 1.02x. Test: `test_row_major_factor_is_the_column_major_one` -- the factor, `tau` and rank bit for bit at f64 and c64, tall/wide/square, sequential and threaded, backwards strides, and a zero leading column. A one-row matrix is left on the direct path (it eliminates at most one column). The `perf-warn` that told row-major callers to transpose is gone. | KEEP | tlinalg-rs#46 | Heap allocation inside a public call, the same boundary the fork's two row splits already cross. Own code: no OpenBLAS or MKL source is used or adapted. Exactly the adapted shapes (tall or square) are the column-major factorisation bit for bit; the guards' exclusions are not. Why, and what lost: `docs/worklogs/2026-10-15-qr-row-major-layout.md`. | | Row-major col-pivoting QR | `qr::col_pivoting::factor::qr_in_place` computes the norm of every one of the `n` columns before it starts and re-reduces them per step, and its incremental-norm SIMD path is gated on `A.row_stride() == 1` and unwraps a column-major view, so a row-major matrix pays one cache line per element on every shape. It now factors a column-major copy and writes the factor back, wide included. Paired, best of 21, CPUs 4 / 4-7 / 4-11, `f64`, C order before -> after (F order after): `256x256` 58.96 / 59.63 / 59.61 -> 1.52 / 1.40 / 1.19 (1.45 / 1.33 / 1.54), `1024x64` 22.23 / 22.29 / 22.45 -> 0.726 / 0.520 / 0.509 (0.472 / 0.522 / 0.536); `c64` `256x256` one thread 64.40 -> 3.29, `1024x64` 24.69 -> 1.64. So the C/F ratio at `256x256` one thread is 38.8x -> 1.05 and the C order is *faster* than the F order at eight threads. Nothing in the campaign records a column-pivoted QR, so this is the record's own driver. Test: `linalg::tests::test_row_major_adaptation_matches_the_column_major_one`. | KEEP | tlinalg-rs#46 | Heap allocation inside a public call, as for the no-pivoting row. Why, and what lost: `docs/worklogs/2026-10-15-row-major-adaptation-lu-colpiv.md`. | -| Row-major LU, partial pivoting | `lu::partial_pivoting::factor::lu_in_place` walks each pivot column twice per step (search, then scaling) and swaps elements *inside* a column, which is a cache line per element when the columns are not contiguous. It now factors a column-major copy and writes the factor and its permutations back. Paired, best of 21, `f64`, C order before -> after (F order after): `512x512` 12.14 / 7.02 / 6.34 -> 3.37 / 2.01 / 1.89 (3.04 / 1.88 / 1.81) at 1T / 4T / 8T, `1024x1024` 59.80 / 27.62 / 22.93 -> 22.91 / 10.04 / 8.21 (19.89 / 8.19 / 6.65), `256x256` one thread 1.530 -> 0.486; `c64` `512x512` one thread 18.08 -> 9.50. Test: `linalg::tests::test_row_major_adaptation_matches_the_column_major_one` -- factor, permutations and transposition count bit for bit. The campaign's `lu` cell takes the composite, so `to_owned` below is what moves it. | KEEP | tlinalg-rs#46 | The full-pivoting flavour is deliberately **not** adapted: it already factors a row-major input through `mat.transpose_mut()`, and the copy measured a no-op there (512x512 f64 one thread 13.758 against 13.712 C order); `10ee1c9` adapted it after all, behind the same `row_major_copy_pays` floor. Since tensor4all/faer#314 the square row-major matrices below 128 with contiguous rows (`col_stride == 1`) skip that copy and go through the transpose dual; generic strides keep the copy. Timed at 64, 80, 96, 112, 120 and 127 on Apple M5 Max (f64/c64, 1T/4T, `32c615f` against `b16d6b8`, per-repetition median of three: dual 0.89-1.02x a column-major input, the copy 0.98-1.34x) and at 64 and 96 on Zen 5 (0.94-1.04x); the sizes in between were not timed. From 128 up and on every tall shape the copy stays, because the dual's strided row scaling (and, tall, its strided long swap and wide trailing block) costs more than the copy. On a magnitude tie the dual picks the smallest row first, the column-major search the smallest column, so tied pivots differ by layout (documented on `lu_in_place`, "Pivot ties"). Guard `full_pivoting::factor::row_major_takes_the_copy`; tests `faer/tests/full_piv_lu_row_major_no_copy.rs` (no `m * n` allocation on the dual shapes, the copy still seen on 128x128, 256x32 and generic-stride 96x96; the tie case) and the full-pivoting arm of `test_row_major_adaptation_matches_the_column_major_one` (127x127 added). Why: `docs/worklogs/2026-11-09-full-piv-lu-row-major-dual.md`. Why: `docs/worklogs/2026-10-15-row-major-adaptation-lu-colpiv.md`. | +| Row-major LU, partial pivoting | `lu::partial_pivoting::factor::lu_in_place` walks each pivot column twice per step (search, then scaling) and swaps elements *inside* a column, which is a cache line per element when the columns are not contiguous. It now factors a column-major copy and writes the factor and its permutations back. Paired, best of 21, `f64`, C order before -> after (F order after): `512x512` 12.14 / 7.02 / 6.34 -> 3.37 / 2.01 / 1.89 (3.04 / 1.88 / 1.81) at 1T / 4T / 8T, `1024x1024` 59.80 / 27.62 / 22.93 -> 22.91 / 10.04 / 8.21 (19.89 / 8.19 / 6.65), `256x256` one thread 1.530 -> 0.486; `c64` `512x512` one thread 18.08 -> 9.50. Test: `linalg::tests::test_row_major_adaptation_matches_the_column_major_one` -- factor, permutations and transposition count bit for bit. The campaign's `lu` cell takes the composite, so `to_owned` below is what moves it. | KEEP | tlinalg-rs#46 | The full-pivoting flavour is deliberately **not** adapted: it already factors a row-major input through `mat.transpose_mut()`, and the copy measured a no-op there (512x512 f64 one thread 13.758 against 13.712 C order); `10ee1c9` adapted it after all, behind the same `row_major_copy_pays` floor. Since tensor4all/faer#314 the square row-major matrices below 128 with contiguous rows (`col_stride == 1`) skip that copy and go through the transpose dual; generic strides keep the copy. Timed at 64, 80, 96, 112, 120 and 127 on Apple M5 Max (f64/c64, 1T/4T, `32c615f` against `b16d6b8`, per-repetition median of three: dual 0.89-1.02x a column-major input, the copy 0.98-1.34x) and at 64 and 96 on Zen 5 (0.94-1.04x); the sizes in between were not timed. From 128 up and on every tall shape the copy stays, because the dual's strided row scaling (and, tall, its strided long swap and wide trailing block) costs more than the copy. On a magnitude tie a row-major input that is not copied (this exclusion, but also wide and small ones) is searched transposed and takes the smallest row first, a column-major or copied one the smallest column, so tied pivots differ by layout (documented on `lu_in_place`, "Pivot ties"). The SIMD search's accumulator merge now keeps the smaller row on a tie within a column, for every layout (it kept whichever accumulator merged first; inherited from upstream). Guard `full_pivoting::factor::row_major_takes_the_copy`; tests `faer/tests/full_piv_lu_row_major_no_copy.rs` (no `m * n` allocation on the dual shapes, the copy still seen on 128x128, 256x32 and generic-stride 96x96; the tie cases) and `pivot_search_keeps_the_smaller_row_across_accumulators` and the full-pivoting arm of `test_row_major_adaptation_matches_the_column_major_one` (127x127 added). Why: `docs/worklogs/2026-11-09-full-piv-lu-row-major-dual.md`. Why: `docs/worklogs/2026-10-15-row-major-adaptation-lu-colpiv.md`. | | Row-major `to_owned` | `MatRef::to_owned` builds its result with `Mat::from_fn`, which fills column by column, so a row-major source is read with a stride of `ncols` and every element costs its own cache line: `512x512` f64 1.01 ms against 0.057, `1024x1024` 6.36 against 0.36, `4096x256` 5.08 against 0.37, `256x4096` 4.30 against 0.36. It now stages a tile when the source's columns are not contiguous and the source is already canonical. After: 0.275 / 2.007 / 2.045 / 1.774 ms. This is what the campaign's `lu` and `qr` cells pay -- a probe at `n = 512` f64 one thread C order reads `lu` 6.1417 -> 4.156 and 4.5329 -> 3.033 at eight threads, `qr` 12.5143 -> 11.185 -- and `full_piv_lu` 12.95 -> 12.63, `cholesky`, `eigh` and `svd_values` unchanged, because those three copy through their own routes rather than `to_owned`. | KEEP | tlinalg-rs#46 | The tile adds a `Mat::zeros` pass and is serial (`to_owned` has no `par`). Gated on `T::IS_CANONICAL`: `canonical()` reinterprets memory and does not apply `Conj::apply`, so a conjugated source (`A.adjoint().to_owned()`, which the sparse Cholesky tests use) keeps the elementwise path -- the ungated attempt broke five of them. The attempted generic `cloned` tile is **dropped**: initializing every destination entry by cloning the first source element amplifies heap-owned payloads by `m * n` before overwrite. `cloned` retains the original elementwise implementation; `test_cloned_does_not_replicate_the_first_heap_value` checks that each source is cloned once. The `Option` tile still moves each cloned numeric value into `to_owned`/factorization destinations without a second clone. Why: `docs/worklogs/2026-10-15-row-major-adaptation-lu-colpiv.md`. | | Row-major `copy_from_triangular_lower` | `MatMut::copy_from_triangular_lower` (what `Llt::new` and `Ldlt::new` copy their input with) ran the elementwise zip, which for a row-major source into a column-major destination reads with a stride of the source's row length: `1024x1024` f64 one thread 4.16 ms against 0.30 ms from a column-major source. It now copies each 16-column band as its diagonal tile plus the rectangle under it through `copy_across_layouts_block`, the tile `to_owned` uses: 0.65 ms. The campaign's `cholesky` cell, 1024 f64 8T, three processes per arm against MKL in the same process: row-major 5.445 -> 2.236 ms (MKL 3.33 / 3.39), column-major 1.674 -> 1.649 (MKL 2.65 / 2.67). Tests: `linalg::tests::tiled_triangular_copy_writes_what_the_elementwise_one_does` and `linalg::tests::llt_and_partial_piv_lu_factor_both_orders_bit_for_bit` (the factor of the same matrix through the two orders, bit for bit). | KEEP | faer-benchmark slide re-take (1024, 8T) | Serial, like `to_owned`: the method has no `par`. A conjugated source keeps the elementwise path. `copy_from` (the full-matrix copy) still runs elementwise: 8.6 ms row-major against 0.47 ms at `1024x1024`; it is not on the `cholesky` or `lu` path. Why: `docs/worklogs/2026-11-06-row-major-triangular-copy.md`. | | 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. | diff --git a/docs/worklogs/2026-11-09-full-piv-lu-row-major-dual.md b/docs/worklogs/2026-11-09-full-piv-lu-row-major-dual.md index 21395a64..647ade0a 100644 --- a/docs/worklogs/2026-11-09-full-piv-lu-row-major-dual.md +++ b/docs/worklogs/2026-11-09-full-piv-lu-row-major-dual.md @@ -30,16 +30,36 @@ near-square shapes keep the copy. - No change to the dual kernel in this slice, and no change to `row_major_copy_pays` or to the other families that share it. -- **Pivot ties depend on the layout.** The pivot search keeps the first strict maximum - of the matrix it factors (smallest column, then smallest row). The dual factors the - transpose, so a row-major input breaks a tie by the smallest original *row* first. - Example (`faer/tests/full_piv_lu_row_major_no_copy.rs`): zero 64x64 f64 with +- **Pivot ties depend on the layout.** The search runs on the input as it is when + `|row_stride| < |col_stride|` (or on its column-major copy when the copy is taken), + and on its transpose otherwise, and keeps the smallest column, then the smallest row, + of that matrix. So a row-major input that is not copied breaks a tie by row first. The + square exclusion is not the only such case: wide row-major matrices and those below + `64 * 64` elements are never copied and were searched transposed before this change + too (#342 review: 8x16 row-major picks `(0, 4)`, column-major `(1, 2)`). Example + (`faer/tests/full_piv_lu_row_major_no_copy.rs`): zero 64x64 f64 with `A[(1, 2)] = A[(0, 4)] = 2`, `max_rank = Some(1)`, `Par::Seq` pivots on `(1, 2)` with two transpositions column-major and on `(0, 4)` with one row-major. #314 allows a different pivot order on ties if it is documented; it is, under "Pivot ties" on - `lu_in_place`. On the adaptation test's fixtures, which hit no such tie, the - pivots, permutations and transposition counts agree; that is all the test shows. The - first version of this change said "the same pivots" without the tie exception. + `lu_in_place`, which now describes that dispatch rather than the copy exclusion alone. + On the adaptation test's fixtures, which hit no such tie, the pivots, permutations and + transposition counts agree; that is all the test shows. +- **The SIMD search did not keep the smaller row on a tie within one column** + (inherited from upstream, found by the second #342 review). Each column is searched + with four interleaved accumulators (three in the fused update), and they were merged + with a strict `>`, so an equal maximum in the accumulator merged first won even with + the larger row: zero 64x64 row-major with `A[(2, 6)] = A[(2, 8)] = 2` pivoted on + `(2, 8)` on Apple M5 (two f64 lanes), and column-major `A[(6, 2)] = A[(8, 2)]` on + `(8, 2)`. `best_score` now takes the other accumulator on a strictly larger score or + on an equal one with a smaller row. Red/green: `pivot_search_keeps_the_smaller_row_across_accumulators` + (every row pair of a 64-row column, both searches, three alignments) failed at rows + `(0, 2)` with row 2 returned, and the integration test returned `(2, 8, 2)`, before the + fix; both pass after it. This changes which pivot is taken on such ties for every + layout, including column-major; it is the order the scalar search and the existing + doc comment in `reduce_2d` already described. Not timed: the merge runs once per + column per step, against an O(m) column walk. The parallel step's merge across + column blocks keeps the earlier block with a strict `>`, which is the smaller column; + no tie test runs on a team. ## Where the dual loses (stage breakdown, M5 Max, f64, one thread) @@ -76,10 +96,16 @@ Per-stage timers inside `lu_in_place_unblocked` (local patch, not committed): arm of `test_row_major_adaptation_matches_the_column_major_one` now asks the new guard whether the two orders must agree bit for bit (copy) or up to rounding (dual), and covers 127x127. -- Whole crate suite (`cargo test -p t4a-faer --no-fail-fast`, this host): 957 passed, - 1 failed, `lq_allocates_nothing_60x300` (allocation counts `[13, 12, 17, 17]`), which - fails the same way on `b16d6b8` and is one of the host failures listed in +- Whole crate suite (`cargo test -p t4a-faer --no-fail-fast`, this host): the one + failure at this change's heads is `lq_allocates_nothing_60x300` (allocation counts + `[13, 12, 17, 17]`); the same command on `b16d6b8` (run once, 2026-11-09, separate + target directory) failed on that test alone, with 955 passing. The counts are in the + #342 description. It is one of the host failures listed in `2026-11-04-gevd-repeated-complex-pairs.md`. +- Padded and reversed row-major views (`(100, 1)`, `(-100, 1)`) go through the dual by + the guard; the allocation test shows that for `(100, 1)` and the adaptation test + covers reversed rows for correctness. Their speed was **not measured**; the timed + cells are contiguous row-major only. - **Measurement 2 (the one that covers the merged code)**: the source of `32c615f` (clean committed tree; `faer/src` tree `38ce727`, unchanged by the later doc-only commits) against `main` `b16d6b8`, both built into fbench `01a1cb8` with the same diff --git a/faer/src/linalg/lu/full_pivoting/factor.rs b/faer/src/linalg/lu/full_pivoting/factor.rs index 609c98d5..83fa19d3 100644 --- a/faer/src/linalg/lu/full_pivoting/factor.rs +++ b/faer/src/linalg/lu/full_pivoting/factor.rs @@ -30,7 +30,16 @@ fn best_score( score: RealReg>, indices: T::SimdIndex, ) -> (RealReg>, T::SimdIndex) { - let is_better = (**simd).gt(score, best_score); + // merges two accumulators of one column, which hold interleaved rows: on equal + // scores the smaller row wins, whichever accumulator holds it, so a tie goes to + // the smallest row as in the scalar search (tensor4all/faer#342) + let is_better = (**simd).or_mask( + (**simd).gt(score, best_score), + (**simd).and_mask( + (**simd).ge(score, best_score), + T::simd_index_less_than(&simd.0, indices, best_indices), + ), + ); ( RealReg(simd.select(is_better, score.0, best_score.0)), simd.iselect(is_better, indices, best_indices), @@ -618,22 +627,32 @@ pub struct FullPivLuInfo { } /// # Pivot ties /// -/// When several entries share the largest magnitude, the pivot depends on the -/// layout. The search scans the matrix it factors column by column and keeps -/// the first strict maximum: the smallest column, then the smallest row. A -/// row-major square matrix below `128` with contiguous rows is factored through -/// its transpose (the condition is `row_major_takes_the_copy` in this module's -/// source), so there the tie goes to the smallest *row*, then the smallest -/// column, of the original matrix. Every other input -- column-major, generic -/// strides, or a row-major matrix that is copied to column-major first -- uses -/// the column-then-row order. +/// When several entries of the remaining block share the largest magnitude, +/// the pivot depends on the layout, because the search runs on one of two +/// matrices: +/// +/// 1. If `row_major_takes_the_copy` holds (in this module's source: `|row +/// stride| != 1`, tall or square, at least `64 * 64` elements, not a square +/// matrix below `128` with `col_stride == 1`), the input is copied to +/// column-major and the copy is searched as it is. +/// 2. Otherwise the input itself is searched as it is when `|row_stride| < +/// |col_stride|`, and its transpose is searched when `|row_stride| >= +/// |col_stride|`. +/// +/// Among the tied entries of the matrix that is searched, the pivot is the one +/// in the smallest column, then the smallest row, of the remaining block at +/// that step (positions after the earlier swaps). So a column-major input, or a +/// row-major one that is copied, breaks a tie by column first, and a row-major +/// input that is not copied (a square one below `128`, a wide one, one below +/// `64 * 64` elements) breaks it by row first, on the SIMD search as on the +/// scalar one. /// /// For example, a zero `64 x 64` f64 matrix with `A[(1, 2)] = A[(0, 4)] = 2` /// and `max_rank = Some(1)` pivots on `(1, 2)` (two transpositions) when -/// column-major and on `(0, 4)` (one transposition) when row-major. Both are -/// valid full-pivoting factorizations of the same matrix, but the -/// permutations, the transposition count and so the factors differ between the -/// two layouts. +/// column-major and on `(0, 4)` (one transposition) when row-major; so does an +/// `8 x 16` matrix with the same two entries. Both are valid full-pivoting +/// factorizations of the same matrix, but the permutations, the transposition +/// count and so the factors differ between the two layouts. pub fn lu_in_place<'out, I: Index, T: ComplexField>( mat: MatMut<'_, T>, row_perm: &'out mut [I], @@ -886,6 +905,45 @@ mod tests { } } + /// Within one column the SIMD search splits the rows over four accumulators (three in the + /// fused update), lane by lane. Two equal maxima in one column can land in different + /// accumulators with the larger row in the one merged first; the merge must still keep the + /// smaller row (tensor4all/faer#342 review: `A[(6, 2)] = A[(8, 2)]` returned row 8 on a + /// two-lane build). Every pair of rows is tried, so the test holds whatever the lane count. + #[test] + fn pivot_search_keeps_the_smaller_row_across_accumulators() { + let m = 64; + for r0 in 0..m { + for r1 in r0 + 1..m { + let mut a = crate::Mat::::zeros(m, 3); + a[(r0, 1)] = 2.0; + a[(r1, 1)] = 2.0; + assert_eq!( + best_in_matrix(a.as_ref()), + (r0, 1, 2.0), + "rows {r0}, {r1}" + ); + + // the fused update with a zero update leaves the matrix as it is + let lhs = crate::Col::::zeros(m); + let rhs = crate::Row::::zeros(3); + for align in [0, 1, 3] { + let mut b = a.clone(); + assert_eq!( + rank_one_update_and_best_in_matrix( + b.as_mut(), + lhs.as_ref(), + rhs.as_ref(), + align + ), + (r0, 1, 2.0), + "update, rows {r0}, {r1}, align {align}" + ); + } + } + } + } + #[test] fn test_gh238() { #[rustfmt::skip] diff --git a/faer/tests/full_piv_lu_row_major_no_copy.rs b/faer/tests/full_piv_lu_row_major_no_copy.rs index 57a9c0ae..80a5cbc3 100644 --- a/faer/tests/full_piv_lu_row_major_no_copy.rs +++ b/faer/tests/full_piv_lu_row_major_no_copy.rs @@ -162,60 +162,72 @@ fn small_square_row_major_full_piv_lu_does_not_copy_the_matrix() { faer::team::with_team(4, |par| run(par)); } -/// The pivot search keeps the first strict maximum of the matrix it factors: -/// the smallest column, then the smallest row. The transpose dual factors the -/// transpose, so a row-major input breaks the same tie by the smallest row -/// first (`lu_in_place`, "Pivot ties"). -#[test] -fn full_piv_lu_tie_order_depends_on_the_layout() { - let n = 64; - let at = |i: usize, j: usize| { - if (i, j) == (1, 2) || (i, j) == (0, 4) { - 2.0 - } else { - 0.0 - } +/// The first pivot `lu_in_place` picks on a zero matrix with two entries of +/// magnitude 2, row-major or column-major, at `max_rank = Some(1)`: `(row, +/// column, transposition count)`. +fn first_pivot( + m: usize, + n: usize, + ties: [(usize, usize); 2], + row_major: bool, +) -> (usize, usize, usize) { + let mut a = vec![0.0f64; m * n]; + for (i, j) in ties { + a[if row_major { i * n + j } else { i + j * m }] = 2.0; + } + let view = if row_major { + MatMut::from_row_major_slice_mut(&mut a, m, n) + } else { + MatMut::from_column_major_slice_mut(&mut a, m, n) }; - let first_pivot = |row_major: bool| { - let mut a = vec![0.0f64; n * n]; - for i in 0..n { - for j in 0..n { - a[if row_major { i * n + j } else { i + j * n }] = at(i, j); - } - } - let view = if row_major { - MatMut::from_row_major_slice_mut(&mut a, n, n) - } else { - MatMut::from_column_major_slice_mut(&mut a, n, n) - }; - let (mut p, mut pi, mut q, mut qi) = ( - vec![0usize; n], - vec![0usize; n], - vec![0usize; n], - vec![0usize; n], - ); - let mut mem = MemBuffer::new(flu::lu_in_place_scratch::( - n, - n, - Par::Seq.shape(), - Default::default(), - )); - let params = flu::FullPivLuParams { - max_rank: Some(1), - ..>::auto() - }; - let (info, _, _) = flu::lu_in_place( - view, - &mut p, - &mut pi, - &mut q, - &mut qi, - Par::Seq, - MemStack::new(&mut mem), - params.into(), - ); - (p[0], q[0], info.transposition_count) + let (mut p, mut pi, mut q, mut qi) = ( + vec![0usize; m], + vec![0usize; m], + vec![0usize; n], + vec![0usize; n], + ); + let mut mem = MemBuffer::new(flu::lu_in_place_scratch::( + m, + n, + Par::Seq.shape(), + Default::default(), + )); + let params = flu::FullPivLuParams { + max_rank: Some(1), + ..>::auto() }; - assert_eq!(first_pivot(false), (1, 2, 2)); - assert_eq!(first_pivot(true), (0, 4, 1)); + let (info, _, _) = flu::lu_in_place( + view, + &mut p, + &mut pi, + &mut q, + &mut qi, + Par::Seq, + MemStack::new(&mut mem), + params.into(), + ); + (p[0], q[0], info.transposition_count) +} + +/// `lu_in_place`, "Pivot ties": the matrix that is searched is the input when +/// `|row_stride| < |col_stride|` (or its column-major copy) and its transpose +/// otherwise; a tie goes to the smallest column, then the smallest row, of the +/// matrix searched. +#[test] +fn full_piv_lu_tie_order_depends_on_the_layout() { + // two columns: the searched column wins + let ties = [(1, 2), (0, 4)]; + assert_eq!(first_pivot(64, 64, ties, false), (1, 2, 2)); + assert_eq!(first_pivot(64, 64, ties, true), (0, 4, 1)); + // a wide row-major matrix is never copied, so it is searched transposed + // too (tensor4all/faer#342 review) + assert_eq!(first_pivot(8, 16, ties, false), (1, 2, 2)); + assert_eq!(first_pivot(8, 16, ties, true), (0, 4, 1)); + // a tall row-major matrix that takes the copy is searched column-major + assert_eq!(first_pivot(256, 32, ties, true), (1, 2, 2)); + + // one searched column, two rows: the smaller row wins, across the SIMD + // accumulators (tensor4all/faer#342 review: (2, 8) was returned) + assert_eq!(first_pivot(64, 64, [(2, 6), (2, 8)], true), (2, 6, 2)); + assert_eq!(first_pivot(64, 64, [(6, 2), (8, 2)], false), (6, 2, 2)); } From 809eaa2991a0d11d936ce5ad11b4243e8eaa42c1 Mon Sep 17 00:00:00 2001 From: Hiroshi Shinaoka Date: Mon, 9 Nov 2026 19:49:27 +0900 Subject: [PATCH 5/6] docs(lu): say the #314 timings predate the tie fix; record the base suite run Co-Authored-By: Claude Opus 5.5 --- docs/worklogs/2026-11-09-full-piv-lu-row-major-dual.md | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/docs/worklogs/2026-11-09-full-piv-lu-row-major-dual.md b/docs/worklogs/2026-11-09-full-piv-lu-row-major-dual.md index 647ade0a..55f18e50 100644 --- a/docs/worklogs/2026-11-09-full-piv-lu-row-major-dual.md +++ b/docs/worklogs/2026-11-09-full-piv-lu-row-major-dual.md @@ -106,9 +106,11 @@ Per-stage timers inside `lu_in_place_unblocked` (local patch, not committed): the guard; the allocation test shows that for `(100, 1)` and the adaptation test covers reversed rows for correctness. Their speed was **not measured**; the timed cells are contiguous row-major only. -- **Measurement 2 (the one that covers the merged code)**: the source of `32c615f` - (clean committed tree; `faer/src` tree `38ce727`, unchanged by the later doc-only - commits) against `main` `b16d6b8`, both built into fbench `01a1cb8` with the same +- **Measurement 2 (the guard as merged, before the tie fix)**: the source of `32c615f` + (clean committed tree; `faer/src` tree `38ce727`). `63790b0` then changed `faer/src` + again (the `best_score` tie fix, one more mask per accumulator merge, and doc + comments); that head was **not re-timed**, so the numbers below do not cover it. + It was measured against `main` `b16d6b8`, both built into fbench `01a1cb8` with the same one-line local `set_threads` patch (returns Ok for `NAME == "none"`), family `full_piv_lu`, row `faer-fork`, the `in_place` line, `--reps 5 --prime-ms 250 --statistic median`. Cells: square 64, 80, 96, 112, 120, 127 (changed) and 128, 256 From 00bbde97dafcb0181bf137c982a5011ec2350b26 Mon Sep 17 00:00:00 2001 From: Hiroshi Shinaoka Date: Mon, 9 Nov 2026 20:25:22 +0900 Subject: [PATCH 6/6] docs(lu): the complex pivot tie follows the search's score, and test it (#342) --- faer/src/linalg/lu/full_pivoting/factor.rs | 26 ++++++-- faer/tests/full_piv_lu_row_major_no_copy.rs | 74 +++++++++++++++++++++ 2 files changed, 93 insertions(+), 7 deletions(-) diff --git a/faer/src/linalg/lu/full_pivoting/factor.rs b/faer/src/linalg/lu/full_pivoting/factor.rs index 83fa19d3..a7fe026a 100644 --- a/faer/src/linalg/lu/full_pivoting/factor.rs +++ b/faer/src/linalg/lu/full_pivoting/factor.rs @@ -639,13 +639,25 @@ pub struct FullPivLuInfo { /// |col_stride|`, and its transpose is searched when `|row_stride| >= /// |col_stride|`. /// -/// Among the tied entries of the matrix that is searched, the pivot is the one -/// in the smallest column, then the smallest row, of the remaining block at -/// that step (positions after the earlier swaps). So a column-major input, or a -/// row-major one that is copied, breaks a tie by column first, and a row-major -/// input that is not copied (a square one below `128`, a wide one, one below -/// `64 * 64` elements) breaks it by row first, on the SIMD search as on the -/// scalar one. +/// Among the tied entries of the matrix that is searched, the pivot is, for a +/// **real** `T`, the one in the smallest column, then the smallest row, of the +/// remaining block at that step (positions after the earlier swaps), on the SIMD +/// search as on the scalar one. So a column-major input, or a row-major one that +/// is copied, breaks a tie by column first, and a row-major input that is not +/// copied (a square one below `128`, a wide one, one below `64 * 64` elements) +/// breaks it by row first. +/// +/// For a **complex** `T` the index order among tied entries is not part of the +/// contract, because the two searches do not score an entry the same way. The +/// scalar one compares `|re| + |im|`; the SIMD one, used when `T` has SIMD +/// support and the matrix searched is column-major, compares the components of +/// the absolute value separately, so an entry with the larger real or imaginary +/// part wins over one with a larger `|re| + |im|`. On an Apple M5 a zero `16 x +/// 16` c64 matrix with `A[(0, 1)] = 3 + 4i` and `A[(2, 1)] = 5` and `max_rank = +/// Some(1)` pivots on `(2, 1)` under SIMD and on `(0, 1)` without it, although +/// `|re| + |im|` is larger at `(0, 1)`. Ties are still resolved inside the +/// matrix that is searched, so the layout argument above holds; only the index +/// order within one searched matrix is unspecified for complex `T`. /// /// For example, a zero `64 x 64` f64 matrix with `A[(1, 2)] = A[(0, 4)] = 2` /// and `max_rank = Some(1)` pivots on `(1, 2)` (two transpositions) when diff --git a/faer/tests/full_piv_lu_row_major_no_copy.rs b/faer/tests/full_piv_lu_row_major_no_copy.rs index 80a5cbc3..b029cb8a 100644 --- a/faer/tests/full_piv_lu_row_major_no_copy.rs +++ b/faer/tests/full_piv_lu_row_major_no_copy.rs @@ -6,6 +6,7 @@ //! The shapes that still take the copy are the control: there the same counter //! must see it. #![cfg(all(feature = "std", feature = "threadteam"))] +use faer::c64; use faer::dyn_stack::{MemBuffer, MemStack}; use faer::linalg::lu::full_pivoting::factor as flu; use faer::{Auto, MatMut, Par}; @@ -231,3 +232,76 @@ fn full_piv_lu_tie_order_depends_on_the_layout() { assert_eq!(first_pivot(64, 64, [(2, 6), (2, 8)], true), (2, 6, 2)); assert_eq!(first_pivot(64, 64, [(6, 2), (8, 2)], false), (6, 2, 2)); } + +/// First pivot `(row, column, transpositions)` on a zero `m x n` c64 matrix with +/// two entries, at `max_rank = Some(1)`. +fn first_pivot_c64( + m: usize, + n: usize, + ties: [(usize, usize, c64); 2], + row_major: bool, +) -> (usize, usize, usize) { + let mut a = vec![c64::new(0.0, 0.0); m * n]; + for (i, j, v) in ties { + a[if row_major { i * n + j } else { i + j * m }] = v; + } + let view = if row_major { + MatMut::from_row_major_slice_mut(&mut a, m, n) + } else { + MatMut::from_column_major_slice_mut(&mut a, m, n) + }; + let (mut p, mut pi, mut q, mut qi) = ( + vec![0usize; m], + vec![0usize; m], + vec![0usize; n], + vec![0usize; n], + ); + let mut mem = MemBuffer::new(flu::lu_in_place_scratch::( + m, + n, + Par::Seq.shape(), + Default::default(), + )); + let params = flu::FullPivLuParams { + max_rank: Some(1), + ..>::auto() + }; + let (info, _, _) = flu::lu_in_place( + view, + &mut p, + &mut pi, + &mut q, + &mut qi, + Par::Seq, + MemStack::new(&mut mem), + params.into(), + ); + (p[0], q[0], info.transposition_count) +} + +/// `lu_in_place`, "Pivot ties" for a complex `T`: the index order is not part of +/// the contract, and the score the search uses is visible. `(0, 1)` holds +/// `3 + 4i` and `(2, 1)` holds `5`: `|re| + |im|` is larger at `(0, 1)`, the +/// largest single component at `(2, 1)`. The scalar search compares the sum, the +/// SIMD one (column-major input) the components, so the pivot differs by build +/// (tensor4all/faer#342 review). +#[test] +fn full_piv_lu_complex_tie_follows_the_search_score() { + use faer::traits::ComplexField; + let ties = [(0, 1, c64::new(3.0, 4.0)), (2, 1, c64::new(5.0, 0.0))]; + let simd = const { ::SIMD_CAPABILITIES.is_simd() }; + assert_eq!( + first_pivot_c64(16, 16, ties, false), + if simd { (2, 1, 2) } else { (0, 1, 1) } + ); + // the same two entries swapped: both searches keep the larger component + assert_eq!( + first_pivot_c64( + 16, + 16, + [(0, 1, c64::new(5.0, 0.0)), (2, 1, c64::new(3.0, 4.0))], + false + ), + (0, 1, 1) + ); +}