diff --git a/CHANGELOG.md b/CHANGELOG.md index 662ce5ad..3cfafc0a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -17,8 +17,10 @@ Changes since `t4a-faer` 0.32.0 (tag `faer-v0.32.0`). Numbers are tensor4all/fae `StackReq::or(colamd::order_scratch(..), ..)` either way, so the ordering was paid for twice. On the issue's reproducer (n = 2000, nnz = 7990, identity ordering, `--release`) the largest single allocation of a simplicial analysis drops from 332712 B, which was exactly the `colamd::order_scratch` request, to 239928 B, which is the analysis workspace (92784 B less); the supernodal arm still - allocates 468080 B, which did not change; that figure is #337's reported run and was not re-measured here, and its call site - was not traced, so which helper it belongs to is not asserted. The `None` path requests what it + allocates 468080 B, which did not change. That figure was attributed afterwards (measured at `e756f13` and at `b989876`, two runs + each): it is not a four-array request but one `Vec` of 58510 elements, allocated three times per symbolic analysis inside + `cholesky::supernodal::ghost_factorize_supernodal_symbolic`; the 64000 B helper only sizes a `MemBuffer` and never appears as its own + allocation. See tensor4all/faer#340. The `None` path requests what it requested before, and no caller `MemStack` variant is added: that is an API decision, not part of this fix (#322). ## Correctness diff --git a/PATCHSET.md b/PATCHSET.md index 562ddf0a..fd9a765c 100644 --- a/PATCHSET.md +++ b/PATCHSET.md @@ -192,7 +192,7 @@ gate has changed and two kinds of adoption are now in the ledger. | Zero-sized elements in the owned reservation and cleanup paths (#286) | `RawMatUnit::try_with_capacity` pads `row_capacity` to a multiple of `align / size` so the row stride in bytes is a multiple of the allocation's alignment. `Mat` has no numeric bound, so for a zero-sized `T` `size_of::()` is 0 while `align_of::()` is at least 1: `align > size` held and `align / size` was a division by zero, so `safe Mat::<()>::with_capacity(1, 1)` panicked with `attempt to divide by zero` at `faer/src/mat/matown.rs:76`. A zero-sized element's byte stride is 0 for every capacity, so the padding is skipped when `size == 0`; nothing else changes, including the empty `layout` and the no-allocation branch the size arithmetic below already takes. Making that path reachable exposed a second defect on the same invariant: `DropMat` advances its column pointer in bytes but `init_with` and `truncate`'s column branch passed it the element count (`faer/src/mat/matown.rs:359`, `:545`), so an aligned or zero-sized element was dropped at a misaligned address (Miri: unaligned reference at `:274`; the aligned non-ZST case is already UB on `a055727`, i.e. not ZST-specific). Both now pass `row_capacity * size_of::()`, as the other two `DropMat` sites already did. Both lines are upstream faer (`c01bf0b0`, `e95ba701`, still in 0.31.0): upstream has the same defects, was not reported, and feedback is deferred. | EXPERIMENT | tensor4all/faer#286 | Correctness only, no timing claim (one branch outside any inner loop, one multiply at cleanup). Tests `zero_sized_elements_reserve_without_panicking` (`Mat::<()>::with_capacity(1, 1)`, `Mat::<()>::with_capacity(0, 0)`, ZST growth, and an `f64` reservation) and `cleanup_of_aligned_and_zero_sized_elements_uses_a_byte_stride` (a `#[repr(align(8))]` drop-counting element, a `#[repr(align(64))]` ZST) fail before the change and pass after; `cargo +nightly miri test -p t4a-faer --test owned_capacity` is clean. Record: [worklog](docs/worklogs/2026-11-07-zst-reservation-division.md). | | `norm_max` propagates NaN (#316) | `linalg::reductions::norm_max` (and `MatRef`/`ColRef`/`RowRef::norm_max`) skipped NaN entries: the SIMD kernel accumulates with `simd.max`, a `select(acc > x, acc, x)` that keeps the non-NaN operand, and the pairwise halves and the scalar loop combine with `fmax` (IEEE `maxNum`). It returned the largest finite magnitude, and 0 for an all-NaN input, while `norm_l1`/`norm_l2` returned NaN. It now returns NaN when any entry, or either part of a complex entry, is NaN, as LAPACK's `xLANGE('M')` and Julia's `norm(x, Inf)` do. The SIMD kernel keeps its max accumulators and adds a second set that sums the magnitudes: every term is `>= 0`, `+inf` or NaN, so the sum is NaN exactly when an entry is NaN, and overflow gives `inf`, not NaN. The pairwise combine propagates NaN, and the scalar loop returns NaN at the first NaN entry. An empty matrix still returns 0. This changes behaviour: internal callers that take `norm_max` of a matrix holding NaN (QZ, tridiagonal EVD, bidiagonal SVD, Bunch-Kaufman scaling) now see NaN instead of a finite norm. Both files are byte-identical to the upstream sync base `0ecd846`, so upstream faer carries the same behaviour there; whether later upstream changed it is not checked, and feedback is deferred. | EXPERIMENT | tensor4all/faer#316 | Correctness only, no timing claim. The SIMD inner loop does one extra vector add per batch. New tests in `linalg::reductions::norm_max::tests`: f32/f64/c32/c64 with NaN in a real entry and in the real or imaginary part of a complex entry, at the first, middle and last position, and all-NaN, across shapes on both sides of the SIMD tail and the pairwise split. Each case goes through the contiguous (SIMD) path and a row-stride-2 view (scalar loop), and through the `MatRef`/`ColRef`/`RowRef` methods. Finite and `inf` inputs give the same result on both paths, and `norm_l1`/`norm_l2` agree on NaN. Removing the SIMD sum line or the scalar NaN return fails `test_norm_max_propagates_nan`. | | `stats` column/row reductions with a reversed output axis (#317) | `col_mean`/`col_varm`, and through them `row_mean`/`row_varm`, normalize a negative row stride with `mat.reverse_rows()` before the kernel runs. Rows are the output axis of a column reduction, and `out` (and, for `varm`, `col_mean`) were not reversed with them. Entry `i` was written to `nrows - 1 - i`, and `varm` subtracted the mean of the mirrored row. All four internal functions (`col_{mean,varm}_{ignore,propagate}`) now reverse `out`, and `col_mean` for `varm`, together with `mat`. Reversing the reduced axis (`reverse_cols` in `col_*`) was already harmless. Public API unchanged. `meanvar.rs` is byte-identical to the upstream sync base `0ecd846`, so upstream faer carries the same defect there; later upstream is not checked, and feedback is deferred. | EXPERIMENT | tensor4all/faer#317 | Correctness only, no timing claim. New tests in `stats::meanvar::tests`: `col_*` on `reverse_rows()` and `row_*` on `reverse_cols()`, mean and varm, `Propagate` and `Ignore`, f32/f64/c32/c64, on column-major (fallback kernel) and row-major (row-major SIMD kernel) storage with one NaN. Each is compared with the same matrix copied to an owned one. The issue's reproducer is pinned as well. Removing the output reversal from any one of the four functions fails both tests. | -| Explicit sparse column ordering stops paying for COLAMD (#322) | `factorize_symbolic_lu_with_column_ordering` and `factorize_symbolic_qr_with_column_ordering` build their internal `MemBuffer` from `StackReq::or(colamd::order_scratch::(m, n, A_nnz), ..)`, but COLAMD runs only in the `None` branch, so a caller that brings the ordering still allocates a COLAMD-sized workspace on every symbolic analysis (dyn-stack's `MemBuffer` takes it uninitialized, so the cost is the allocation, not a memset). Both functions include `order_scratch` only when `column_ordering.is_none()` now; the `None` request and the analysis request are unchanged. Both entry points are this fork's own (#239, `3889bd1`), so upstream faer does not carry the defect and no `upstream-reports/` entry is added. | KEEP | tensor4all/faer#322 | Allocation only, no numerical change: the two functions' permutations, symbolic structures and every downstream numeric path are untouched, and no benchmark was run. The issue's reproducer (n = 2000, nnz = 7990, identity ordering, `--release`, counting global allocator) goes from 332712 B largest allocation to 239928 B for simplicial LU and QR: the 332712 B was exactly `colamd::order_scratch`, and 239928 B is the analysis request (4 + 5 n-element arrays, the transpose arrays, and the largest simplicial arm), while the supernodal arm stays at 468080 B. Its call site was not traced and it is not attributed to any helper here. The figure itself is #337's reported run of issue #322's counting program and was not re-measured for this change. `faer/tests/sparse_symbolic_colamd_workspace.rs` pins both directions on a 100 x 100 sample: with an explicit ordering the largest allocation stays below `order_scratch`, without one it reaches it; the first fails before the change and both pass after. `cargo test -p t4a-faer -j6` passes except `lq_no_alloc::lq_allocates_nothing_60x300`, which fails identically on `origin/main` (dense LQ, unrelated), plus `no_step_alloc`, which counts heap allocations, and `triangular_solve::aliased_rhs::staged_solve_is_bitwise_the_in_place_one`, which counts staged solve nodes through `stage::STAGED` (`faer/src/linalg/triangular_solve.rs:509,579`) and asserts that a 1 KiB stride was **not** staged (the message names the failure), not that anything allocated. Both pass standalone at this revision and both were seen to fail inside the whole-suite run; what makes them fail there was not identified. `no_step_alloc` is its own test binary, so another test binary cannot perturb its counter directly, which rules out the obvious explanation without supplying a better one. | +| Explicit sparse column ordering stops paying for COLAMD (#322) | `factorize_symbolic_lu_with_column_ordering` and `factorize_symbolic_qr_with_column_ordering` build their internal `MemBuffer` from `StackReq::or(colamd::order_scratch::(m, n, A_nnz), ..)`, but COLAMD runs only in the `None` branch, so a caller that brings the ordering still allocates a COLAMD-sized workspace on every symbolic analysis (dyn-stack's `MemBuffer` takes it uninitialized, so the cost is the allocation, not a memset). Both functions include `order_scratch` only when `column_ordering.is_none()` now; the `None` request and the analysis request are unchanged. Both entry points are this fork's own (#239, `3889bd1`), so upstream faer does not carry the defect and no `upstream-reports/` entry is added. | KEEP | tensor4all/faer#322 | Allocation only, no numerical change: the two functions' permutations, symbolic structures and every downstream numeric path are untouched, and no benchmark was run. The issue's reproducer (n = 2000, nnz = 7990, identity ordering, `--release`, counting global allocator) goes from 332712 B largest allocation to 239928 B for simplicial LU and QR: the 332712 B was exactly `colamd::order_scratch`, and 239928 B is the analysis request (4 + 5 n-element arrays, the transpose arrays, and the largest simplicial arm), while the supernodal arm stays at 468080 B. It was attributed afterwards (tensor4all/faer#340, measured at `e756f13` and `b989876`, two runs each): it is one `Vec` of 58510 elements, allocated three times per symbolic analysis inside `cholesky::supernodal::ghost_factorize_supernodal_symbolic`; the 64000 B four-array helper only sizes a `MemBuffer` and never appears as an allocation of its own. The figure itself is #337's reported run of issue #322's counting program and was not re-measured for this change. `faer/tests/sparse_symbolic_colamd_workspace.rs` pins both directions on a 100 x 100 sample: with an explicit ordering the largest allocation stays below `order_scratch`, without one it reaches it; the first fails before the change and both pass after. `cargo test -p t4a-faer -j6` passes except `lq_no_alloc::lq_allocates_nothing_60x300`, which fails identically on `origin/main` (dense LQ, unrelated), plus `no_step_alloc`, which counts heap allocations, and `triangular_solve::aliased_rhs::staged_solve_is_bitwise_the_in_place_one`, which counts staged solve nodes through `stage::STAGED` (`faer/src/linalg/triangular_solve.rs:509,579`) and asserts that a 1 KiB stride was **not** staged (the message names the failure), not that anything allocated. Both pass standalone at this revision and both were seen to fail inside the whole-suite run; what makes them fail there was not identified. `no_step_alloc` is its own test binary, so another test binary cannot perturb its counter directly, which rules out the obvious explanation without supplying a better one. | | Complex `abs_max` drops a NaN in the real part (#320) | `ComplexField::simd_abs_max` takes the larger of the two part magnitudes with an unordered `re > im` select, and that select returns `im` when `re` is NaN: a NaN real part read as the finite imaginary magnitude, while a NaN imaginary part did return NaN. The imaginary side already lands on `im`, so the real side takes an explicit `!(re == re)` mask in the generic `Complex` SIMD path, and the scalar fallbacks of `ComplexImpl` / `ComplexImpl` test both parts before the comparison. No in-tree type reaches either path (native `c32` / `c64` go through the real coercion in `norm_max`), so no dispatch, allocation or numerical result on an existing input changes. | KEEP | tensor4all/faer#320 | Correctness; the `faer-traits` unit test `complex_abs_max_propagates_nan` drives `SimdCtx::abs_max` on `Complex` and fails before / passes after (real-part NaN, imaginary-part NaN, both, and three finite controls). No performance claim: two mask ops per SIMD chunk on a path no in-tree type reaches. | | Complex QZ: a zero chased down `B`'s diagonal updates the columns of `Q` (#324) | When `hessenberg_to_qz_unblocked` (the unblocked complex QZ path, which the default parameters run for `n < max(15, blocking_threshold)`, i.e. `n < 75`, and which the blocked path also calls at the end) finds `T[(j, j)]` negligible with `j == ilo` or a negligible `H[(j, j - 1)]`, it chases the zero to the bottom with rotations on rows `jch, jch + 1` of `H` and `T`. `Q` must take the adjoint rotation on its columns, as the QZ sweep, the other zero-chase branch and the blocked path do (LAPACK `zhgeqz`: `ZROT(N, Q(1,JCH), 1, Q(1,JCH+1), 1, C, DCONJG(S))`); it rotated rows `jch, jch + 1` of `Q` with `s` instead. `Q` stayed unitary but `A = Q S Z^H` failed (`O(1)`), so every left eigenvector of `gevd_cplx` was wrong while the eigenvalues and the right eigenvectors, which use only `S`, `T` and `Z`, were right. It now rotates columns `jch, jch + 1` with `conj(s)`. | KEEP | tensor4all/faer#324 | Correctness only, no timing and no performance claim: the same number of rotations on a column pair instead of a row pair. Whether the branch runs depends on rounding: the issue's 12 x 12 `c32` pencil takes it on x86-64 with AVX2/FMA (left residual `1.962e-1` before, `2.079e-7` after; right `1.889e-7` both) and not on aarch64 or on x86-64 without AVX (Rosetta), where it passes before and after; `c64` does not take it. Tests: `linalg::gevd::qz_cplx::tests::unblocked_zero_chase_updates_q_columns` gives the unblocked iteration an exact zero at `j == ilo` and below a zero subdiagonal at `c32`/`c64`, `n = 2, 3, 6, 12`, and checks `Q`, `Z` unitary and `A = Q S Z^H`, `B = Q T Z^H` to `10 n eps` (fails before on aarch64: `2.51` against `7.0e-15`); `faer/tests/gevd_infinite_eigenvalue.rs` embeds the issue's bit patterns and checks the left and right residuals at `c32` and `c64` (fails before on x86-64 AVX2 only); before folding a column into the maxima it normalizes the eigenvector by its norm and `(alpha, beta)` by `max(|alpha|, |beta|)`, which leaves the residual unchanged, and requires the two eigenvector norms, that scale and the residual denominator to be finite and positive; a numerator formed from the values as they come loses the residual instead, and no later division brings it back: `u = v = 1e-200 * e0` with `alpha = 1e-150` underflows it to `0`, and `alpha = f64::MAX / 2` with `beta = 0` sends the denominator to `inf`, whose quotient is `0` as well. It then asserts both residuals finite, with the column number and values in the message, because `f64::max`/`f64::min` drop a `NaN` and the `0/0` of an all-zero eigenvector column would otherwise read as a perfect residual. Measured: with `ul = Mat::::zeros(n, n)` after the `gevd_cplx` call the test passed before these asserts and fails after with `column 0: left eigenvector has norm 0e0, alpha = Complex { re: -2.6592135429382324, im: 2.4517483711242676 }, beta = Complex { re: 0.0, im: 0.0 }`, and the same mutation on `ur` fails with `column 0: right eigenvector has norm 0e0`; replacing each column's eigenvectors with `k * e0`, `k = 1.125 * (f64::MAX / den)`, passed with the single-division form (`nu * den = inf`, residual recorded as `0e0` against a staged `2.978e-1`) and fails with `c32 left residual 2.977833357222273e-1 >= 1.430511474609375e-5` after it; and `residual_oracle_rejects_unnormalized_eigenpairs` hands the oracle the review's two columns directly (`alpha = f64::MAX / 2, beta = 0, u = v = 1e-300 * e0` and `alpha = 1e-150, beta = 0, u = v = 1e-200 * e0`): on the oracle without normalization both read `left residual 0e0 < 1.430511474609375e-5` and `right residual 0e0 < 1.430511474609375e-5`, and with it `2.977833357222273e-1` and `2.941381743992797e-1`. The dispatch claim was measured too, with a temporary `eprintln!` at both `hessenberg_to_qz_unblocked` call sites and default `params`: a 20 x 20 pencil prints `n=20 nmin=75 blocking_threshold=75 -> unblocked`, and a 100 x 100 one prints `n=100 nmin=75 -> final unblocked call after the blocked loop`. The same code is in upstream faer `main` `5833f4e` at `faer/src/linalg/gevd/qz_cplx/mod.rs:123-127`, from `57e25ac` (generalized evd initial impl); verified by reading that source, not by running it, and no `upstream-reports/` entry is added (the deferred-feedback decision of `AGENTS.md`). |