matmul: overwrite a MaybeUninit destination (#136)

PATCHSET.md

@@ -133,6 +133,7 @@ gate has changed and two kinds of adoption are now in the ledger.
 | Team dispatch cost read from threadteam | `utils::cost::team_section_cost` and `MEASURED_MAX_WIDTH` (faer's own linear extrapolation above width 8, #97) are deleted: the extrapolation lives in threadteam (`WidthCosts` 0.1.1, linear above a row's largest anchor; 0.1.2, a measured Zen2 team row for `x86_64-avx2` at widths 2 to 64). `ParShape::dispatch_cost` reads `DispatchCosts::cost(Team, w)`; `par_pays_time`, `team::kernel::kernel_width` and the tiled team product's floor ask `TeamExec::dispatch_cost(w)`. `kernel_width` searches from the full width down instead of stopping at the first narrow width that does not pay (measured rows fall with the width). Requirement `threadteam = "0.1.2"`. Decisions change at widths 7 and 8 for rows whose largest anchor is 6, at width 8 on `x86_64-avx512` (its 0.98 us anchor) and at every width on `x86_64-avx2`. | EXPERIMENT (pending the maintainer's acceptance decision) | tensor4all/threadteam-rs#6, tensor4all/faer#91 | No timing (gate inputs only; acceptance run left to the maintainer). Tests: `utils::cost::tests::the_team_cost_is_threadteams_table_at_every_width`, `utils::team::tests::kernel_width_finds_the_widest_paying_width_on_a_falling_row`. Record: [worklog](docs/worklogs/2026-10-29-team-cost-from-threadteam.md). |
 | Every region at the widest width that pays | The full-width gates (`par_pays`, `ParThreshold::region_pays`, `ParShape::pays_full` at the call sites, the block apply's `n_tasks`) chose between the token's whole degree and one thread, so on a wide team a region that pays at 8 or 16 ran serially. Every region now asks `utils::cost::par_width` / `ParThreshold::region_width` for the widest width that pays (every width searched, falling rows included; the GEMM kernel's `kernel_width` is the same rule) and runs on a sub-team of that width (`utils::thread::with_width`, a split into one child: no worker dispatch). Splits whose task count is capped (triangular solve jobs of 64 columns, block apply columns, `upgrade_householder_factor`'s rows per participant) search inside their cap. The factorizations that latched `par = Par::Seq` for the rest of the call (Hessenberg, tridiagonal, unblocked bidiagonal, full pivoting LU, Bunch-Kaufman, column pivoting QR) decide per step; `Never` still makes the whole call sequential; the post-loop compact-factor upgrades of the Hessenberg and tridiagonal reductions now ask their own gate (they ran under the last step's latch, i.e. sequentially). The blocked bidiagonal panel gate is a floor (sequential or the whole `par`), so the panel's products choose their own widths. The Bunch-Kaufman blocked swap loop, split over the whole degree without a gate, now has one (row swaps at the LU swap rate). `Fixed` keeps its rules (the whole degree above the floor). `par_pays`, `region_pays` and the Seq-only GEMM predicate are removed; `par-trace` records the chosen width. | EXPERIMENT (pending the maintainer's acceptance run on broccoli / evo) | tensor4all/faer#91, tensor4all/faer#95 | No timing here (performance evidence is the monitor's: broccoli 1/8/16/32/64T, evo). Tests: `utils::team::tests::{par_width_takes_the_widest_width_that_pays, with_width_runs_on_a_sub_team_of_that_width, the_families_narrow_their_regions_and_match_sequential}`, the rewritten `region_rules_follow_the_serial_time`, `fixed_and_never_keep_their_meaning`, `par_width_reads_the_degree_of_the_token`. Record: [worklog](docs/worklogs/2026-10-29-width-per-region.md). |
 | No heap allocation inside a parallel region (#111) | Every region on the team route allocated on the heap. `for_each_child`, behind every `for_each` and `for_each_with_par`, collected all the items into a `Vec`, one entry per column for `par_col_iter_mut`, and then a slot `Vec`. The full-pivoting LU step and both Bunch-Kaufman full-pivoting updates built an `nthreads`-long vector of pivot candidates per step, and `join_sized` boxed its two closures. The allocations were suspected of causing #109's floor-level loss (`full_piv_lu f64 128` at 8T on evo, +6.6 to +10.4%). The evo re-measurement of this change shows they were not: the cell is at 1.0025 of `main` with them gone. The change stands on its own as a cleanup of the mechanism. Now faer's dispatch and these workspaces allocate nothing. An iterator adapter that allocates in its own `split_at`, such as rayon's `map_with`, still does. `for_each_child` splits the iterator's producer on the leader into the blocks (`split_at`, remainder-relative) and each block's leader iterates its piece. Up to 64 blocks, with a producer of at most 256 bytes, are dealt from an inline array in one fork; wider or heavier splits halve recursively, one section per level. The pivot candidates come from the call's scratch: `lu_in_place_scratch` and `cholesky_in_place_scratch` grow by `par.degree()` entries. `join_sized` takes two typed slots. | KEEP (evo: no regression; `full_piv_lu f64 128` 8T at 1.0025 of `main`) | tensor4all/faer#111, tensor4all/faer#95 | Test `faer/tests/no_step_alloc.rs`: its own binary with a process-wide counting allocator. On teams of 4, 8 and 70 (the last takes the recursive path), with the regions forced parallel and a warm-up first, it counts zero allocations for full-pivoting LU 96, Bunch-Kaufman with full pivoting 96 and `join_sized`, and checks that each one dispatched. On `origin/main` the same test counts 190 allocations in one full-pivoting LU 96 call at width 4. Record: [worklog](docs/worklogs/2026-10-29-no-region-allocation.md). |
+| Matrix product into uninitialized storage (#136) | `matmul_with_conj_uninit(dst: MatUninitMut, lhs, conj_lhs, rhs, conj_rhs, alpha, par) -> MatMut` computes `alpha * op(lhs) * op(rhs)` into a destination of `MaybeUninit<T>` elements and returns it as an initialized `MatMut`, for callers (strided-rs#198) that would otherwise fill the output before an overwriting product. Requested upstream as sarah-quinones/faer-rs#269; the view type `MatUninitMut` (`UninitMut`) and the first tests come from shinaoka's upstream PR #270 (same author). The destination's prior contents are never read and no `&T`/`&mut T` to an unwritten element is formed. The native types on x86_64 with AVX2 or AVX-512 go to `t4a-private-gemm-x86` with `Accum::Replace` (on a team, through `gemm_team`); other scalars write each dot product through `ptr::write`; an empty depth writes zeros. The `gemm` and `nano-gemm` crates never see uninitialized storage (`gemm`'s GEMV forms a `&mut [T]` of the destination, its Apple AMX kernel copies partial tiles of it before checking `alpha`, and `nano-gemm` scales it by `beta`), so on other targets the native types are zero-filled first and multiplied as by `matmul_with_conj`. Any strides, negative included; a batch is one view per slot. A panic leaves the storage `MaybeUninit` (possibly partly written) with nothing exposed as initialized. Existing `MatMut` APIs are unchanged. | EXPERIMENT | tensor4all/faer#136 | Tests in `linalg::matmul::tests`: NaN-prefilled destinations with padded leading dimensions for f32/f64/c32/c64, row- and column-major, conjugations, at `1x1x1` up to `130x70x90` and across the kernel's 512 depth block (`9x600x7`); the padding must stay NaN and the product must match `matmul_with_conj`. Also a team, a batched buffer, negative strides, an empty depth, a matrix-vector shape and `fx128`. Miri runs the seven small tests; six of them use `MaybeUninit::uninit()` storage (the batched one is NaN-filled), and under Miri the native types take the zero-fill route. No uninitialized read was reported. The NaN and team tests are too large for Miri. Indicative only (primerose, shared, one shape, best of 30, three runs): `2048x64x2048` f64 `Par::Seq` into existing storage is 0.89 to 0.93 of a zero fill followed by `Accum::Replace`. The saving is the output-sized fill, so it shrinks as `k` grows. Not timed on a quiet host. Record: [worklog](docs/worklogs/2026-10-29-matmul-uninit.md). |
 
 ## Ground rules
 

docs/worklogs/2026-10-29-matmul-uninit.md

@@ -0,0 +1,37 @@
+# Matrix product into uninitialized storage (#136)
+
+## Goal
+
+strided-rs#198 holds output buffers as `MaybeUninit<T>` and wants a product that overwrites
+them without first filling them, which the upstream request sarah-quinones/faer-rs#269 asks
+for. This fork adds `matmul_with_conj_uninit` taking a `MatUninitMut` and returning the
+initialized `MatMut`. The view type and the first tests come from shinaoka's upstream PR #270.
+Nothing was posted upstream.
+
+## Decisions
+
+- **Kernel routes.** Overwriting is not enough; the destination must never be read and never
+  be borrowed as `&T`/`&mut T` before it is written.
+  - `t4a-private-gemm-x86` with `Accum::Replace` meets both. It writes every element in its
+    first depth block (later blocks add to what that block wrote) and loads the destination
+    only when accumulating. It is the route for the native types on x86_64 with AVX2 or
+    AVX-512, directly or through `gemm_team` on a team.
+  - The `gemm` crate with `read_dst = false` was the first plan for the other targets. It was
+    dropped after the independent review traced two exceptions in `gemm-common` 0.19: its GEMV
+    builds a `&mut [T]` over the destination before writing, and its experimental Apple AMX
+    kernel copies the destination into a staging tile before it checks `alpha`.
+  - `nano-gemm` scales the destination by `beta`, so it reads it.
+  - On those targets the native types are therefore zero-filled through raw writes and then
+    multiplied by the ordinary path. This is correct everywhere, and gains nothing there.
+- **Other scalars.** One dot product per element, written with `ptr::write`. A team splits
+  the columns.
+- **Panics.** The `MatMut` is formed only after the whole write. A panic leaves the caller's
+  storage `MaybeUninit`, possibly partly written, and never exposed as initialized.
+
+## Limitations
+
+- Miri cannot run the SIMD kernel. Under Miri the native types take the zero-fill route, so
+  Miri checks the generic and fill paths, and the NaN-prefill tests check that the x86 kernel
+  does not carry old contents into the product. Neither proves the absence of a discarded read
+  inside the kernel; that rests on reading its epilogues.
+- The timing is one shape on a shared host (in `PATCHSET.md`), not a quiet-host measurement.

faer/src/lib.rs

@@ -1153,7 +1153,7 @@ pub type cx128 = traits::cx128;
 /// `Complex<f64>`
 pub type fx128 = traits::fx128;
 pub use col::{Col, ColMut, ColRef};
-pub use mat::{Mat, MatMut, MatRef};
+pub use mat::{Mat, MatMut, MatRef, MatUninitMut};
 pub use row::{Row, RowMut, RowRef};
 #[allow(unused_imports, dead_code)]
 mod internal_prelude {
@@ -1253,7 +1253,7 @@ pub mod prelude {
 	pub use super::prelude_sp::*;
 	pub use super::{Par, ParShape, Scale, c32, c64, col, mat, row, unzip, zip};
 	pub use col::{Col, ColMut, ColRef};
-	pub use mat::{Mat, MatMut, MatRef};
+	pub use mat::{Mat, MatMut, MatRef, MatUninitMut};
 	pub use reborrow::{IntoConst, Reborrow, ReborrowMut};
 	pub use row::{Row, RowMut, RowRef};
 	/// see [`Default`]