BLAS beta==0 write-only semantics (fixes 0*NaN=NaN poisoning from cold scratch)

bench/REDUCED_SWEEP_RESULTS.md

@@ -30,16 +30,14 @@ top of the shuffle cost.
 
 ## Win-condition (what the picker encodes)
 
-Reduced is competitive only when **both** hold:
-
-1. `n_out ≤ blockDim / 32` — every output gets its own warp, so there is no
-   serial output loop layered on top of the shuffle; and
-2. `K_contract ≥ 32` — the contraction is long enough to (a) use a full warp's
-   lanes and (b) amortize the shuffle tail.
-
-`glass::suggested_use_reduced<n_out, K_contract, blockDim>()` returns `true` only
-in that corner and `false` otherwise — i.e. it recommends **serial almost
-always**. Treat the `*_reduced` family as a niche tool, not a default.
+**Retired to constant `false` (2026-07-08).** The quiet-GPU resweep measured
+**0 of 48** configurations where reduced clears the ±5% tie margin — even the
+former long-contraction corner (`n_out ≤ blockDim/32 && K_contract ≥ 32`, the
+2×64×2 cells that measured 1.11–1.12× on the earlier shared-GPU sweep)
+collapsed into the noise band. `glass::suggested_use_reduced<>()` therefore
+returns `false` unconditionally on sm_120; its signature is kept as the seam
+for a future data-derived corner on different hardware (e.g. Jetson Orin).
+Treat the `*_reduced` family as an expressiveness tool, never a speed default.
 
 ## Caveat for the tensor / congruence families
 

bench/bench_paper_fusion.cu

@@ -315,6 +315,8 @@ static void run_shape(Handles H, const char* dt) {
         std::vector<Row> rows = {
             {"fused_tb32",  [&]{ fused(32); }},
             {"fused_tb128", [&]{ fused(128); }},
+            {"fused_tb256", [&]{ fused(256); }},   // wide-block probe: the (36,12)/(48,16)
+            {"fused_tb512", [&]{ fused(512); }},   // chain-wins cells may be TB-starved
             {"chain",       chain},
         };
         for (auto& r : rows) {

bench/tune.py

@@ -247,8 +247,13 @@ def build_reduced(sms):
 
 
 def predicate_use_reduced(n_out, K, blockDim):
-    """Mirror of glass::suggested_use_reduced<n_out,K_contract,blockDim>()."""
-    return (n_out <= blockDim // 32) and (K >= 32)
+    """Mirror of glass::suggested_use_reduced<n_out,K_contract,blockDim>().
+
+    Retired to constant False 2026-07-08: the quiet-GPU resweep measured 0/48
+    reduced wins under the ±5% margin, so the predicate's true-corner was
+    emptied. If a sweep on new hardware finds wins, derive a new corner from
+    that data and update BOTH this mirror and the .cuh predicate."""
+    return False
 
 
 def analyze_reduced(text, margin):

docs/agent_debugging_guide.md

@@ -127,26 +127,37 @@ self-adapts to the launch.
 - **How to localize:** run at `THREADS=1` (forces full serialization — if it now passes, you have
   a race or a coverage gap) versus many threads, and bisect.
 
-### 1c. `beta = 0` GEMM still READS C → `0 * NaN = NaN` poisoning
-
-GLASS's beta-form GEMM computes `C[cidx] = alpha*res + beta*C[cidx]` — it **reads `C` even when
-`beta == 0`** (see `gemm_impl` in `src/base/L3/gemm.cuh`, and the docstrings explicitly say
-"C is read; caller must initialize it"). If `C` is an **uninitialized scratch slot**, its leftover
-bit pattern can be `NaN`/`Inf`, and `0 * NaN == NaN` poisons the result even though beta is 0.
-
-- **Where this applies:** every beta-taking path — `glass::gemm` / `gemm_impl`,
-  `gemm_tiled`, `gemm_strided`, and the `_1d` batched variants — when called with
-  `beta = 0` into a destination the caller did not initialize.
-- **The fix:** the **caller must initialize `C`** before a `beta = 0` write into cold scratch
-  (a tiny `for (i=rank; i<size; i+=size) C[i] = 0;` + `__syncthreads()`), OR use a GLASS
-  overload that has no `beta*C` term. `gemm.cuh` provides implicit-`beta=0` overloads (the
-  `gemm` forms documented "with implicit `beta = 0`") that write `C` directly and never
-  read it — prefer those for write-into-fresh-scratch.
-- **The tell:** a discrepancy that is **thread-count-dependent and vanishes when isolated / run at
-  1 thread** is the signature of reading cold scratch (the leftover bytes are nondeterministic
-  across launches; at 1 thread the slot is more reliably overwritten before re-read). Reproduce
-  flakes by running the check many times, not once. `test_gemm_strided` already parametrizes
-  `(alpha,beta) = (1.0, 0.0)` for exactly this reason — keep beta=0 cases in any new GEMM test.
+### 1c. `beta = 0` semantics — FIXED library-wide (2026-07-08): destinations are write-only
+
+**Historical bug class, now closed at the primitive level.** GLASS's beta-form ops used to compute
+`C[cidx] = alpha*res + beta*C[cidx]` unconditionally — reading `C` even when `beta == 0`, so an
+**uninitialized scratch slot** holding leftover `NaN`/`Inf` poisoned the result (`0 * NaN == NaN`).
+This is exactly how GRiD's generated RNEA NaN'd (2026-07-08): the v/a forward-pass gemvs write cold
+`s_vaf` shared memory with `beta = 0`, and a diverged rollout kernel had left NaN on the SM.
+
+- **The guarantee (BLAS/cuBLAS convention, now enforced):** when `beta == 0` the destination is
+  **write-only** — it need not hold a valid value on input. Every beta-taking blend routes through
+  `beta_blend(acc, beta, dst)` in `src/base/barrier.cuh`, which reads `dst` only under
+  `beta != 0`. Covered: `gemv` / `gemv_strided` / `gemv_segmented` / `gemv_reduced`, `gemm`
+  (untiled + tile4 + tiled) / `gemm_strided` / `gemm_reduced`, `syrk` / `syrk_reduced`, `symm`,
+  `congruence`, `axpby` — block, `warp::`, and `cgrps::` surfaces (shared impls).
+- **When you write a NEW beta-taking op:** use `beta_blend`, never a raw `alpha*res + beta*dst`.
+  The no-beta overloads remain preferred for write-into-fresh-scratch (they also skip the
+  runtime compare).
+- **The regression net:** `test_l2.py::test_gemv*_beta0_poisoned_y_no_read` and
+  `test_l3.py::test_gemm_rt_betaform_beta0_no_read` call the BETA overloads with `beta = 0` into a
+  NaN-poisoned destination and require a clean result; `test_gemm_rt_beta0_no_read` covers the
+  no-beta overloads. Keep poisoned-destination beta=0 cases in any new beta-form test.
+- **The `nvidia::` surface honors the same guarantee** — audited 2026-07-08: `l3_simt.cuh`
+  forwards to the fixed base cores, and the cuBLASDx route was probed empirically (poison ∈
+  {0, 1e30, NaN} → bit-identical output, f32+f64): `cublasdx::execute` does not blend a
+  `beta == 0` destination on MathDx 26.03. That vendor behavior is PINNED by
+  `test_nvidia_dispatch.py::test_dispatch_op[beta0_poison]` — if a future MathDx bump breaks it,
+  that test fails and the fix is to zero-fill `c_smem` instead of staging `C` when `beta == 0`.
+- **The historical tell (if you ever see it again):** a discrepancy that is thread-count-dependent,
+  nondeterministic across launches, and vanishes at 1 thread is the signature of reading cold
+  scratch. Sanitizers do NOT catch read-of-uninitialized **shared** memory — qNaN-poison the
+  arena to prove it (see GRiD/PDDP heavy-robustness notes).
 
 ### 1d. Storage-order / layout-flag mistakes
 

docs/source/index.rst

@@ -77,10 +77,11 @@ The **contraction-parallel + fused family** (``gemm_reduced`` / ``gemv_reduced``
 ``syrk_reduced``, the ``tensor_vec_contract`` / ``vec_tensor_vec`` /
 ``congruence_sym`` / ``bilinear`` ops, ``riccati_gain``, and compile-out
 robustness flags on ``potrf`` / ``ldlt`` / ``posv``) is there for
-**expressiveness and fusion only**: the ``*_reduced`` decomposition measured
-*slower* than the plain serial in-thread contraction in **47 of 48 swept shapes**
-(``bench/REDUCED_SWEEP_RESULTS.md``), and the ``suggested_use_reduced<>`` picker
-declines it almost everywhere — **prefer the plain ops for throughput**. See
+**expressiveness and fusion only**: on a quiet GPU the ``*_reduced``
+decomposition measured slower than the plain serial in-thread contraction in
+**48 of 48 swept shapes** (±5% margin, ``bench/REDUCED_SWEEP_RESULTS.md``), and
+the ``suggested_use_reduced<>`` picker now declines it everywhere — **prefer
+the plain ops for throughput**. See
 :doc:`user_guide/concepts/contraction_parallel` for the measurement and
 :doc:`user_guide/concepts/namespaces` for the naming convention.
 

docs/source/user_guide/concepts/contraction_parallel.rst

@@ -27,7 +27,9 @@ All ship in the three SIMT surfaces (``glass::`` block, ``glass::warp::``,
    RTX 5090 (sm_120), the contraction-parallel decomposition is **slower than the
    plain serial in-thread contraction in 47 of 48 swept shapes**, often by
    10–100× (full table: ``bench/REDUCED_SWEEP_RESULTS.md``), and the
-   :cpp:func:`glass::suggested_use_reduced` picker declines it almost everywhere.
+   :cpp:func:`glass::suggested_use_reduced` picker now declines it everywhere
+   (retired to constant ``false`` after the 2026-07-08 quiet resweep measured
+   0/48 wins under the ±5% tie margin).
    **Prefer the plain ops** (``gemm`` / ``gemv`` / ``syrk``) for throughput;
    reach for this family only for the fused forms (``tensor_vec_contract``,
    ``vec_tensor_vec``, ``congruence_sym``, ``bilinear``) that the serial surface
@@ -62,14 +64,16 @@ every configuration** — 47 of 48 swept shapes lose, often by 10–100×. The s
 ``gemm`` over shared-resident data is a tight per-thread loop that is very hard to
 beat at these sizes, while ``*_reduced`` pays a ~5-step shuffle latency per output
 and, at the typical short contraction (K = 14–21), leaves most of a warp's lanes
-idle. The *only* win in the swept space was ``n_out = 4, K = 64`` at ``blockDim
-≥ 128`` — and only by ~1.1×. So there is no "~14×"; the realistic story is "use
-serial unless you are in a tiny, well-characterized corner."
-
-:cpp:func:`glass::suggested_use_reduced` encodes that corner —
-``n_out <= blockDim/32`` (every output owns a warp) **and** ``K_contract >= 32``
-(long enough to fill a warp and amortize the shuffle) — and returns ``false``
-otherwise, i.e. recommends serial almost always:
+idle. The earlier shared-GPU sweep showed one marginal win (``n_out = 4,
+K = 64`` at ``blockDim ≥ 128``, ~1.1×); the quiet-GPU resweep of 2026-07-08
+collapsed even that cell into the ±5% noise band — **0 of 48 configurations
+pick reduced**. So there is no "~14×", and no corner either; the realistic
+story is "use serial."
+
+:cpp:func:`glass::suggested_use_reduced` encodes that measurement — it returns
+``false`` unconditionally on sm_120, and keeps its ``<n_out, K_contract,
+blockDim>`` signature as the seam where a retune on different hardware (e.g.
+Jetson Orin) can reinstate a data-derived corner without touching call sites:
 
 .. code-block:: cuda
 

docs/source/user_guide/tutorials/sweep_results.rst

@@ -139,5 +139,39 @@ kernel, composed GLASS calls never pay the API floor again.
 Both harnesses live in ``bench/`` and rerun via
 ``python3 bench/paper_sweeps.py`` (see ``bench/PAPER_SWEEPS.md``).
 
+Choosing among the dense solve paths (measured guidance)
+---------------------------------------------------------
+
+GLASS deliberately ships **no** auto-dispatch between its linear-system
+solvers — the right choice depends on structure and conditioning, which a
+compile-time table cannot see. Instead, ``bench/tune.py --legs solvers``
+measures the trade-offs on your GPU and records them in
+``bench/SOLVERS_SWEEP_RESULTS.md``. The RTX 5090 numbers (NPROB=8192,
+ns/problem):
+
+**SPD single solve — use ``posv``; the alternatives price as follows.**
+``gesv`` (pivoted LU, the robustness fallback) costs 1.0–2.1× ``posv`` in
+fp32 at ``N`` ≥ 16 (1.8–2.1× at 32–64); the invert-then-multiply
+anti-pattern (``inv`` + ``gemv``) costs 1.9–4.7× at ``N`` ≥ 16. **Below**
+``N`` = 16 all three paths are within a few nanoseconds of each other (and
+``inv``/``gesv`` can even edge out ``posv``, especially in fp64) — there,
+choose by *numerics*, not speed: Cholesky is backward-stable on SPD input and
+fails loudly (with ``CHECK``) on indefinite input, while an explicit inverse
+amplifies conditioning error silently. Speed only ever argues *for* ``posv``,
+never against it.
+
+**Block-tridiagonal chains — ``bdsv`` (direct) vs ``pcg`` (iterative) is
+problem-dependent; do not hard-code either.** On our diagonally-dominant test
+system PCG converges in ~3 iterations and wins 11 of 12 cells (up to 9×); at
+(BlockSize=12, Knots=16) fp32 the direct sweep wins 1.8×. PCG's cost scales
+linearly with its iteration count, so an ill-conditioned Riccati/KKT chain
+(10–100× more iterations) moves the crossover proportionally toward
+``bdsv`` — read the ``pcg iters`` column of your own sweep before
+generalizing, or measure with your actual matrices.
+
+**``syev`` / ``eig_clamp``** — the decompose–clamp–reconstruct op costs the
+same as the bare eigensolve (the clamp epilogue is free); budget ~0.9 µs
+fp32 / ~8.8 µs fp64 per 32×32 problem at saturation.
+
 See :doc:`../concepts/tuning` for how to emit a per-host override table from a
 sweep, and :doc:`../../api_reference/defaults` for the picker API.

src/base/L1/axpy.cuh

@@ -62,7 +62,7 @@ template <typename Bar, typename T, bool TRAILING_SYNC = true>
 __device__ void axpby_impl(Bar bar, uint32_t n, T alpha, T *x, T beta, T *y, T *z)
 {
     uint32_t rank = bar.rank(), size = bar.size();
-    for (uint32_t i = rank; i < n; i += size) z[i] = alpha*x[i] + beta*y[i];
+    for (uint32_t i = rank; i < n; i += size) z[i] = beta_blend(alpha*x[i], beta, y[i]);
     if constexpr (TRAILING_SYNC) bar.sync();
 }
 

src/base/L2/gemv.cuh

@@ -14,7 +14,7 @@ __device__ void gemv_impl(uint32_t rank, uint32_t size,
                 T a = ROW_MAJOR_A ? A[col*n + row] : A[col + row*m];
                 res += a * x[col];
             }
-            y[row] = alpha*res + beta*y[row];
+            y[row] = beta_blend(alpha*res, beta, y[row]);
         }
     } else {
         for (uint32_t row = rank; row < m; row += size) {
@@ -23,7 +23,7 @@ __device__ void gemv_impl(uint32_t rank, uint32_t size,
                 T a = ROW_MAJOR_A ? A[row*n + col] : A[row + col*m];
                 res += a * x[col];
             }
-            y[row] = alpha*res + beta*y[row];
+            y[row] = beta_blend(alpha*res, beta, y[row]);
         }
     }
 }
@@ -129,7 +129,7 @@ __device__ void gemv_impl_ct(uint32_t rank, uint32_t size,
                 T a = ROW_MAJOR_A ? A[col*N + row] : A[col + row*M];
                 res += a * x[col];
             }
-            y[row] = alpha*res + beta*y[row];
+            y[row] = beta_blend(alpha*res, beta, y[row]);
         }
     } else {
         for (uint32_t row = rank; row < M; row += size) {
@@ -138,7 +138,7 @@ __device__ void gemv_impl_ct(uint32_t rank, uint32_t size,
                 T a = ROW_MAJOR_A ? A[row*N + col] : A[row + col*M];
                 res += a * x[col];
             }
-            y[row] = alpha*res + beta*y[row];
+            y[row] = beta_blend(alpha*res, beta, y[row]);
         }
     }
 }
@@ -238,8 +238,8 @@ namespace warp {
      * of the `M×N` matrix `A` (each row an independent inner product). Set
      * `TRANSPOSE=true` for `Aᵀ * x` and `ROW_MAJOR=true` for row-major `A`. No shared
      * scratch, no `__syncthreads`; independent warps may run distinct problems
-     * concurrently. Full 32 lanes required. `C`/`y` is read (the `beta * y` term);
-     * use the no-beta overload to write into uninitialized destinations. NumPy
+     * concurrently. Full 32 lanes required. `y` is read only when `beta != 0`
+     * (BLAS semantics: `beta == 0` treats `y` as write-only). NumPy
      * equivalent: `y = alpha*A@x + beta*y` (or `alpha*A.T@x + beta*y` when transposed).
      *
      * @tparam T          Scalar type (e.g. `float`, `double`).

src/base/L2/gemv_reduced.cuh

@@ -33,7 +33,7 @@ namespace detail {
                     p[vlane] = acc;
                 }
                 const T res = reduced_tree32<T>(p);
-                y[o] = HAS_BETA ? (alpha*res + beta*y[o]) : (alpha*res);
+                y[o] = HAS_BETA ? beta_blend(alpha*res, beta, y[o]) : (alpha*res);
             }
             return;
         }
@@ -45,7 +45,7 @@ namespace detail {
                 for (uint32_t c = lane; c < CDIM; c += 32u)
                     partial += (TRANSPOSE ? A[c + o*M] : A[o + c*M]) * x[c];
                 const T res = warp::reduce<T>(partial);
-                if (lane == 0) y[o] = HAS_BETA ? (alpha*res + beta*y[o]) : (alpha*res);
+                if (lane == 0) y[o] = HAS_BETA ? beta_blend(alpha*res, beta, y[o]) : (alpha*res);
             }
         }
     }
@@ -66,7 +66,7 @@ namespace detail {
  * @param alpha  Scalar on the product.
  * @param A      Input matrix (M x N, column-major).
  * @param x      Input vector (length N if !TRANSPOSE else M).
- * @param beta   Scalar on the existing y (read; caller must initialize it).
+ * @param beta   Scalar on the existing y (read only when `beta != 0`).
  * @param y      In/out result (length M if !TRANSPOSE else N).
  */
 template <typename T, uint32_t M, uint32_t N, bool TRANSPOSE = false, bool TRAILING_SYNC = true>

src/base/L2/gemv_segmented.cuh

@@ -130,7 +130,7 @@ __device__ void gemv_segmented(
         if (ATOMIC_Y) {
             atomicAdd(&y[y_base + out], alpha * res);
         } else {
-            T val = alpha * res + beta * y[y_base + out];
+            T val = beta_blend(alpha * res, beta, y[y_base + out]);
             if (FUSE_SCALED_ADD)
                 val += S[static_cast<uint32_t>(seg_s_off[seg]) + out] * scalar[seg];
             y[y_base + out] = val;

src/base/L2/gemv_strided.cuh

@@ -37,7 +37,7 @@ __device__ void gemv_strided(T alpha, const T* A, const T* x, T beta, T* y)
         T res = static_cast<T>(0);
         for (uint32_t col = 0; col < N; col++)
             res += A[row + col * ROW_STRIDE] * x[col];
-        y[row] = alpha * res + beta * y[row];
+        y[row] = beta_blend(alpha * res, beta, y[row]);
     }
 }
 

src/base/L3/congruence.cuh

@@ -48,10 +48,10 @@ namespace detail {
                 }
                 const T res = reduced_tree32<T>(p);
                 const uint32_t idx = i + j*ROWS;
-                R[idx] = HAS_BETA ? (alpha*res + beta*R[idx]) : (alpha*res);
+                R[idx] = HAS_BETA ? beta_blend(alpha*res, beta, R[idx]) : (alpha*res);
                 if (SYMM && i != j) {
                     const uint32_t m = j + i*ROWS;
-                    R[m] = HAS_BETA ? (alpha*res + beta*R[m]) : (alpha*res);
+                    R[m] = HAS_BETA ? beta_blend(alpha*res, beta, R[m]) : (alpha*res);
                 }
             }
             return;
@@ -69,10 +69,10 @@ namespace detail {
                 const T res = warp::reduce<T>(partial);
                 if (lane == 0) {
                     const uint32_t idx = i + j*ROWS;
-                    R[idx] = HAS_BETA ? (alpha*res + beta*R[idx]) : (alpha*res);
+                    R[idx] = HAS_BETA ? beta_blend(alpha*res, beta, R[idx]) : (alpha*res);
                     if (SYMM && i != j) {
                         const uint32_t m = j + i*ROWS;
-                        R[m] = HAS_BETA ? (alpha*res + beta*R[m]) : (alpha*res);
+                        R[m] = HAS_BETA ? beta_blend(alpha*res, beta, R[m]) : (alpha*res);
                     }
                 }
             }