*2026-07-09. Compiled from a systematic review of the tree covering the scalar cores, per-ISA packetmath, GPU paths, math functions, core integration, sparse, Tensor execution/RNG, tests/CI, the external toolchain landscape, and the open-issue/MR backlog. File:line references are against master @ c4bed222b (2026-07-09) and may drift. Filed as a tracking work item; individual MRs per the sequence in §9 can be split out as child items. (Edit 1: rescoped Appendix A — hardware-local wide-accumulation kernels are in scope, only the generic GEBP framework redesign is deferred.)* **Scope:** element-wise operations — the scalar types, conversions, `numext` math, packetmath across all backends, GPU device paths, and the test/CI story — **plus hardware-local wide-accumulation product kernels** where an ISA provides widening bf16/fp16 operations and the existing per-backend override seams suffice (Appendix A.1). Only the **generic GEBP/GEMM framework redesign** for wider-than-result accumulators is deferred to a follow-up project (Appendix A.2). Rationale: the GEBP kernel framework requires significant redesign to support an accumulator type wider than the input/output type, but wherever hardware-accelerated bf16/fp16 operations with wider accumulation can be used *locally* — the pattern the AltiVec MMA kernel already proves in-tree — we should use them now. --- ## 1. Baselines: what Eigen can assume in mid-2026 | Area | Floor today | External state | Consequence | |---|---|---|---| | C++ | C++14 (`cxx_std_14`) | nvcc 13.3 still accepts `--std=c++14`; `_Float16`/`__fp16`/`__bf16` are compiler extensions usable under C++14 | Baseline holds. `std::float16_t`/`std::bfloat16_t` (C++23, GCC≥13 libstdc++) are optional-detection candidates only. C++17 bump is separately in flight (!2465 / #3050) | | CUDA | 11.4, sm_70 (`Macros.h:593,601`) | CUDA 13.3 current; 13.x dropped Maxwell/Pascal/Volta → sm_75 min there; `cuda_bf16.h` ships since 11.0; **CUDA ≥ 12.2 enables host-side emulated `__half`/`__nv_bfloat16` operators + implicit integer conversions by default** | `cuda_bf16.h` can be included unconditionally. sm_70-only paths are legacy once built with 13.x. Host operator coexistence needs an explicit audit (§5.4) | | HIP/ROCm | 5.6 (`Macros.h:617`) | ROCm 7.2.4 production, deployed base 6.3/6.4; `hip_bf16.h` (`__hip_bfloat16`, CUDA-parity API) is the documented path; the legacy `hip_bfloat16.h` class Eigen builds on is legacy-but-not-deprecated; `hip_bf16.h` had multiple-definition link bugs through ≈6.0.2 | Migrate HIP bf16 to `hip_bf16.h` gated **ROCm ≥ 6.1**, keep the legacy path below that. Raising the HIP floor to 6.x is realistic if needed | | x86 | F16C keyed on `__F16C__` | F16C universal on all AVX2 CPUs; AVX512FP16 = Sapphire/Granite Rapids only (AMD first in Zen6); AVX512BF16 = Cooper Lake, SPR+, Zen4/5 — **convert + dot-product only, no bf16 arithmetic exists on any x86 CPU** | Treat F16C as the fast-conversion baseline, AVX512FP16 as niche/server. x86 bf16 strategy = convert + `vdpbf16ps` (f32 accumulate; GEMM follow-up) | | ARM | `__ARM_FEATURE_FP16_*` | v8.2 fp16 arithmetic ubiquitous (all 2018+ silicon); v8.6 bf16 (`BFDOT`/`BFMMLA`) is **widening-only into f32** (Graviton3, Neoverse V1/N2, Apple M2+); SME2 in Apple M4 | fp16 packets assumable on all relevant AArch64. Any ARM bf16 compute must be bf16-in/f32-accumulate | | RISC-V | `zvfh`+`zfh` → RVV10FP16 | Zvfh ratified but optional (RVA23 mandates only Zvfhmin=conversions); bf16 exts (Zvfbfmin/Zvfbfwma) ratified, conversion + widening-FMA-to-f32 only | Same story as ARM/x86: no native bf16 arithmetic anywhere; f32 is the universal compute type for bf16 | | Host compilers | — | `_Float16`: GCC≥12 / Clang≥15 on x86-64 (SSE2), long-standing on AArch64, **absent on MSVC** (keep uint16 emulation forever); `__bf16` arithmetic: GCC≥13 / Clang≥17 | Feature-macro detection, never version-based; MSVC path immortal | A cross-cutting external fact that shapes the follow-up project: **every** contemporary ISA (VDPBF16PS, BFDOT/BFMMLA, Zvfbfwma) and every contemporary BLAS (oneDNN "strict" default, OpenBLAS sbgemm, ACL mixed-precision) treats f32 as the accumulator for 16-bit inputs. Eigen's per-op 16-bit-rounding semantics are the outlier. --- ## 2. State of the tree: what is already right (do-not-touch list) The audit found the *scalar* 16-bit story fundamentally sound, and several backends reference-quality. These are strengths to preserve; any refactor below must not regress them: 1. **Conversions are RTNE-correct everywhere it matters.** float→half (Giesen soft path `Half.h:670-715`, F16C imm=0, `__float2half`, hardware casts) and float→bf16 (`BFloat16.h:431-622`, exhaustively documented) both handle denormals, overflow→inf, and NaN sign preservation; NEON (`NEON/PacketMath.h:4683-4707`) and AltiVec (`AltiVec/PacketMath.h:2007-2061`) reimplement the same RTNE + NaN-canonicalization semantics in vector form. Tested bit-exactly. 2. **The via-float emulation is *provably* correctly rounded.** float's p=24 ≥ 2p+2 for both half (p=11) and bf16 (p=8), so compute-in-float + one rounding gives exact IEEE results for +,−,×,÷,√,fma (Figueroa 1995). This justifies `is_iec559=true` for half and means the x86/NEON-fallback packet paths are bit-identical to native IEEE fp16 — *more* accurate than chained native-fp16 arithmetic. Do not replace via-float single ops with native-precision loops "for speed" without checking accuracy semantics. 3. **The fused cast evaluator** (`CoreEvaluators.h:569-760`) — vectorized, bounds-safe, multi-ratio (1,2,4,8); `a.cast<float>()*b.cast<float>()` is single-pass with no temporary. It is the pattern for the missing NEON/RVV cast entries and is fp8-ready as-is. 4. **AltiVec/VSX bf16 stack** — complete packet set, correct masks via bit-interleave, dedicated GEMM/GEMV with f32 MMA accumulators and runtime P10 dispatch. The in-tree model implementation (relevant to the follow-up project). 5. **RVV10 fp16 packets and NEON native-fp16 branch** — fully native, idiomatic, both recently landed (!2087-!2093, !2273) and CI-covered on real hardware (zvfh board, arm64 runners). 6. **GPU fp16 floors modernization** (5b173498e, !2547) — all pre-sm_70/`__HIP_ARCH_HAS_*` guards are gone, clean `#error`s at the floors; device ops are unconditional native intrinsics. Don't re-add version guards. 7. **Tensor tree-reduction for 16-bit on CPU** (LeafSize<half>=200/bf16=128, `TensorReduction.h:226-236`) with a probabilistic error-bound test — the quality bar the dense `Redux.h` path currently fails (§4.6). 8. **Uniform 16-bit RNG** (`TensorRandom.h:60-76`) — bit-exact k·2⁻¹⁰ / k·2⁻⁷ grids, Sterbenz-exact on every storage backend; and core `random_impl<half/bf16>` correctly delegates to float with 10/7 mantissa bits. 9. **Test infrastructure bones**: bit-exact `VERIFY_*_BITS_EQUAL` macros, NaN-tolerant harness for 16-bit, packetmath's exhaustive cast matrix and exact rounding-op checks, per-MR smoke coverage of half/bf16 scalar+packet suites on x86-64, aarch64 (+fp16), riscv64 (zvfh). 10. **TF ABI constraints are documented and load-bearing**: the host-`__half_raw`-across-phases design (`Half.h:106-136`), the `__half_raw::x` / `__bfloat16_raw::value` member names, the `half_impl` ADL-via-base-class pattern, and the `construct_from_rep_tag` idiom (copied by ml_dtypes) are frozen surface. See §5.2 for the exact list. ### Backend coverage matrix (element-wise, current state) | Backend | half packets | bf16 packets | arithmetic impl | vectorized math fns | 16↔32 casts | segment/masked | |---|---|---|---|---|---|---| | SSE-only | — (scalar) | — (scalar) | — | — | — | — | | AVX/AVX2 | Packet8h | Packet8bf | via-f32 (F16C) | full via-f32 set | yes | no | | AVX512 | Packet16h | Packet16bf | via-f32 | full (bf16 log family needs DQ) | yes | no | | AVX512FP16 | Packet32h/16h/8h native | (Packet16bf) | **native fp16** | partial — erf/ndtri/bessel/sinh/cosh/log10/cbrt lost vs non-FP16 | yes (FP16 casts) | no | | NEON (AArch64) | Packet8hf/4hf native | Packet4bf via-f32 | native / via-f32 | **thin: only ptanh+logistic for hf**; full set for bf16 | half yes, **bf16 no** | n/a | | SVE | **none** | **none** | — | — | — | — | | SME | (defers to NEON, but drops NEON GEBP — §4 bug) | none | — | — | — | — | | AltiVec/VSX | **none** | Packet8bf | via-f32 (P10 hw cvt) | full | yes | partial-load yes | | RVV10 | Packet1Xh/2Xh native | PacketXbf | native / widening | defined but **trait-disabled for fp16 exp/log** | **none** | n/a | | SYCL | cl_half8 etc. | **none** | native sycl::half | broad | yes | **size-8 device gap** | | GPU CUDA/HIP | half2/Packet4h2 native | **none — pure scalar float-emulation even on sm_80/CDNA** | native fp16 | h2exp/h2log/h2sqrt/h2rsqrt; rounding intrinsics unused | yes (fp16) | n/a | | ZVector/LSX/MSA/HVX | none | none | — | — | — | — | --- ## 3. The demand signal (what users actually hit) Small but real: ~8 relevant open issues. The strongest external ask is bf16 GEMV/GEMM via OpenBLAS (#2936, ARM + Google/TF involvement — follow-up project, but note the already-merged sbgemm path has two correctness bugs, Appendix A.6). AVX512FP16 demand shows up as breakage reports from Sapphire-Rapids users (`-march=native`): #2910, #2951, #2829/#2835, #3013 — and that path currently has **zero runtime CI** (!2558 dropped the test jobs for lack of a runner). GPU fp16 issues are toolchain-robustness (#3077 H100 `cuda_bf16.h` ambiguity, #2929 CUDA 12.9). #2937 asks for `numext::nextafter` for half. Nobody has explicitly asked for `__nv_bfloat16` packets — that motivation is downstream (TF/JAX-style users and the dead bf16 GPU tests), not the issue tracker. In-flight work to coordinate with: **!2700** (fast-math constant poisoning — touches half/bf16 packet kernels; land/rebase first), **!2465/#3050** (C++17 bump — don't add C++14 contortions it will delete), **!2159** (gtest migration), **#3024** (GEBP DoublePacket/predux_half refactor — GEMM-adjacent), **!2681** (GPU module refactor). Power bf16 code is unowned since 2023 (ChipKerchner moved to RVV) — flag him for review on anything touching it. --- ## 4. Workstream A — correctness fixes (P0/P1, element-wise) Ordered roughly by severity × blast radius. Items 1–4 are the ones I'd treat as release-blocking quality bugs. ### A1. `NumTraits<Eigen::half>::epsilon()` is 8× too small (P0) `Half.h:896` returns `0x0800` = 2⁻¹³; true fp16 epsilon is 2⁻¹⁰ = `0x1400`, which `std::numeric_limits<half>::epsilon()` (`Half.h:290`) already returns — 96dc37a03 (2022) fixed numeric_limits and missed NumTraits; bf16 is consistent (both 2⁻⁷). Every eps-driven threshold for half is 8× too tight: JacobiSVD sweep/deflation (`JacobiSVD.h:381,486,825`), ColPivQR rank (`ColPivHouseholderQR.h:403-404`), RealSchur deflation, iterative-solver default tolerance, SparseQR pivot threshold. **Fix:** change to `0x1400` (or delete the override and inherit from numeric_limits). **Must land together with A12's SparseQR guard** — the larger epsilon makes SparseQR's already-broken in-type threshold 8× worse. Re-test `EigenSolver<half>` (currently NoConvergence) and half SVD after; note SVD has an independent root cause (A13). `dummy_precision` (`0x211f` ≈ 1e-2) stays byte-identical — it feeds isApprox and all test tolerances. ### A2. aarch64/GPU raw-bit interop performs *numeric* conversion (P0 — regression from !2273, 2026-06-01) Since !2273 made `__fp16` the default `Eigen::half` storage on **all** aarch64 builds (`__ARM_FP16_FORMAT_IEEE` is a default predefine; `Macros.h:393-401`), three raw-bit paths silently value-convert in the host phase of CUDA/HIP builds on Grace/Jetson (and analogously for `_Float16` under AVX512FP16): - `operator __half()` — `hr.x = x` assigns `__fp16` into `unsigned short` (`Half.h:240-246`): `half(1.5)` → `__half` bits `0x0001`. - `raw_uint16_to_half` GPUCC branch — `h.x = x` converts uint16→`__fp16` numerically (`Half.h:625-628`): `NumTraits<half>::infinity()` becomes `half(31744.0)` = inf-by-accident-not-bits… actually a finite garbage value; **every raw-bit constant in NumTraits/numeric_limits/packet code is corrupted**. - HIP `half_base(const __half&)` — `x = __half_as_ushort(h)` (`Half.h:178`), same class of bug. - Additionally `half(__fp16)` references `construct_from_rep_tag` on the *vendor* `__half_raw` in the device pass (`Half.h:215-221`, no `!EIGEN_GPU_COMPILE_PHASE` guard) — probable compile break with clang-as-CUDA-compiler on aarch64. **Fix (recommended):** make every raw-bit site storage-aware — route through the existing bit-preserving tagged constructors / `raw_half_as_uint16` (`Half.h:634-647` already handles all storages) and guard `Half.h:215-221` with `!defined(EIGEN_GPU_COMPILE_PHASE)`. The blunt alternative (zero `EIGEN_HAS_ARM64_FP16`/`EIGEN_HAS_BUILTIN_FLOAT16` under `EIGEN_GPUCC`) restores pre-!2273 behavior but forfeits native host storage on GPU builds. Either way: add a **host-side raw-bit round-trip test** compiled under `EIGEN_CUDACC` (`bit_cast<uint16_t>(h) == static_cast<__half>(h)` raw bits — 10 lines, would have caught this) and a compile-only aarch64+CUDA CI job. The Tensor RNG path is immune (uses `numext::bit_cast`, whole-object copy) — no blocker there. ### A3. RVV10 bf16 comparison masks are corrupted (P0 within its backend) `PacketMathBF16.h:243-260, 579-596` implement `pcmp_*` as `F32ToBf16(pcmp_*<f32>)`; the all-ones f32 true mask is a NaN pattern, and `vfncvtbf16` is an *arithmetic* conversion producing canonical NaN `0x7FC0` instead of `0xFFFF`. Anything consuming masks bitwise (generic `pselect`/`pand`/`psign`) computes garbage. NEON (`F32MaskToBf16Mask`) and AltiVec (`F32ToBf16Bool`) do this correctly. The packetmath test can't see it because `CHECK_CWISE2_MASK` only tests zero/nonzero (`packetmath_test_shared.h:134-144`). **Fix:** integer-truncate the f32 mask (`vnsrl` by 16), mirroring the fp16 path; **and strengthen `CHECK_CWISE2_MASK` to compare full bit patterns** (catches this class everywhere). Note CI never builds this backend today (march string lacks `_zvfbfmin_zvfbfwma` — §7). ### A4. Tensor normal RNG emits NaN/±inf for 16-bit types (P0 for TF-style users) `RandomToTypeNormal<T>` (`TensorRandom.h:196-210`) runs Leva's ratio-of-uniforms in T's precision; the 16-bit uniform is exactly 0 with p = 2⁻¹⁰/2⁻⁷, `log(0) = -inf` poisons the accept test, and measured output rates are 1.5e-5 (half) / 2.4e-4 (bf16) NaN+inf; tails are also truncated at 5.2σ/4.4σ. **Fix:** specialize for half/bf16 as `T(RandomToTypeNormal<float>(...))` — one rounding, float-quality tails; optionally guard u==0 in the float path (residual ~1e-9 rate). Add an all-finite large-N test (deterministic catch at N~1M) and half/bf16 subtests to `tensor_random_gpu.cu`. ### A5. Scalar math accuracy trio (P1) - **`scalar_log2_op`** computes `Scalar(LOG2E)*log(a)` in 16-bit — 3 roundings, 37% of half inputs misrounded, wrong on exact powers of two; packet lanes (via-float) disagree with scalar tail (`UnaryFunctors.h:504-505`). Route real scalars through `numext::log2` (bit-exact on all 31743 inputs); keep the LOG2E form for complex only. Also use `hlog2` on device / `::log2f` on host inside `half_impl::log2` (`Half.h:796-798`). - **`numext::rsqrt`** = `T(1)/sqrt(x)` — two 16-bit roundings; 26%/28% misrounds vs single-rounded float computation; device misses `hrsqrt` (`MathFunctions.h:2146-2157`). Add dedicated half/bf16 overloads (host: `half(1.0f/::sqrtf(float))`; device: `hrsqrt`); re-enable the disabled array_cwise rsqrt test (`test/array_cwise.cpp:213-224` FIXME). - **`numext::hypot`** runs the whole `positive_real_hypot` recurrence in 16-bit (~14% mismatch; also degrades `abs(complex<half>)`) (`MathFunctionsImpl.h:151-171`). Specialize: half → naive float (cannot overflow: 65504²·2 ≪ FLT_MAX); bf16 → run the recurrence in float. ### A6. bfloat16 `isnan`/`isinf`/`isfinite` are fast-math-fragile (P1) `BFloat16.h:634-651` classify via `float(a)` and std functions — constant-folded to false under `-ffinite-math-only`. Half already uses raw bit tests (`Half.h:753-763`). Port the bit-pattern classification to bf16 (`(bits&0x7fff)>0x7f80` etc.). This aligns with the #3045 FTZ/fast-math audit and !2700's hardening direction. ### A7. `blueNorm()` returns 0 for half (P1) Blue's 1978 constants assume a float/double-range exponent; with fp16 (iemin=-13, iemax=16) the "big" rescale `(ax·2⁻¹³)²` underflows to exactly 0 for unit-scale data (`StableNorm.h:109-141`); bf16 is structurally OK but drifts. Dispatch 16-bit Real to float accumulation (or to the `stableNorm` path). Related small item: `hypotNorm` drift for bf16. ### A8. Dense reductions: silent 16-bit accumulation (P1, policy decision) `VectorX<half>::Ones(70000).sum()` = 32768 (packet lanes stall at 2048); `.mean()` = 0 (divisor `Scalar(70000)` = inf) (`Redux.h:298-309,552`); Tensor `MeanReducer` has the same divisor overflow + a count-inflation quirk under tree combine (`TensorFunctors.h:130-143`). The Tensor sum path already does tree reduction and is fine. **Fix (minimum):** compute the mean divisor (or reciprocal) in float; document dense 16-bit `sum()` semantics. **Fix (better):** port the Tensor LeafSize tree-reduction idea to `Redux.h` for sub-32-bit scalars. Changing accumulator *type* is deliberately out of scope (GEMM-adjacent policy); tree reduction in-type is accuracy-improving without changing the contract. Sparse `sum()` mirrors whatever dense does (`SparseRedux.h:22-34`). ### A9. Compile-blockers from int literals & operator ambiguity (P0-trivial) All one-liners, zero risk to other types: - `ConjugateGradient.h:51`, `LeastSquareConjugateGradient.h:52`: `tol_error = 0;` → `RealScalar(0)` (CG for 16-bit currently **cannot be instantiated at all**). - `IncompleteLUT.h:455,473,475`, `IncompleteCholesky.h:266`: literal fills/assignments. - `AmbiVector.h:270,279` (+ default arg at :264): blocks `(A*B).pruned()`. - `Dot.h:201`: `RealScalar(1)/p` ambiguous for half/bf16 (`operator/(half,Index)` vs implicit-float) — breaks `lpNorm<p>` for general p → `RealScalar(1)/RealScalar(p)`. Add compile/smoke tests so these can't regress (§7). A targeted grep for `= 0;` / `= 1;` int-literal assignment to Scalar lvalues in SparseCore/IterativeLinearSolvers would catch stragglers. ### A10. `std::hash<Eigen::half>` breaks the ±0 contract (P2) Raw-bit hash gives hash(+0)≠hash(−0) while +0==−0 (`Half.h:1024-1033`). Canonicalize −0 before hashing; same check for bf16; add tests. ### A11. GPU packet-semantics polish (P2) `pselect(half2)` uses a numeric mask==0 test (−0.0 mask reads as false; partial-bit masks act all-or-nothing) (`GPU/PacketMath.h:718-724`); `pmin`/`pmax` NaN-lane semantics are the mirror image of the CPU scalar convention (`:818-826`); generic `preverse` is identity for multi-lane GPU packets — currently unreachable but a booby trap after 8290a5529 (`TensorReverse.h:204`). Fix alongside the §5.4 GPU cleanup MR. ### A12. Sparse 16-bit numerics (P1, land with A1) - **SparseQR** computes its default pivot threshold `20(m+n)·max2Norm·eps` entirely in RealScalar (`SparseQR.h:452-454`): half n=400 → threshold 4.79 → **rank 0/400 with info=Success**; n=2000 → `RealScalar(80000)`=inf; and `max2Norm` accumulates `squaredNorm` in half → inf for realistic columns. Compute both in float (round only the final value), or `static_assert` SparseQR for 16-bit Real. **The A1 epsilon fix makes this 8× worse if landed alone — same MR.** - **Sparse `isApprox`** squares norms in-type → any half matrix with ‖·‖²_F > 65504 compares approximately-equal to anything (`SparseFuzzy.h:27`); measured `isApprox(C,-C)==true`. Promote the accumulation for 16-bit Real or use scaled norms. - **BiCGSTAB** (the only 16-bit iterative solver that compiles) goes non-finite at data scale ~64 in half via in-type `squaredNorm`/`dot` (`BiCGSTAB.h:69,75,106`). Document the scale bound or promote the three quadratic quantities. - Note: SimplicialLLT/LDLT and SparseLU are epsilon-free and honest in 16-bit up to conditioning limits; `prune` uses dummy_precision, not the broken epsilon — leave alone, point users there. ### A13. `real_2x2_jacobi_svd` overflow makes half SVD garbage (P1, decomposition-adjacent) Near convergence `u=t/d` ~ 1/(2·eps) and `abs2(u)` overflows fp16 (2¹⁸ > 65504) → `c=s=0` "rotation" annihilates the work matrix; JacobiSVD/BDCSVD return all-zero/garbage singular values for half on generic inputs (`Jacobi.h:489-491`; also the benign identity-rotation variant in `makeJacobi`, `Jacobi.h:104-113`). bf16 immune (range). Guard the large-|u| branch (compute via `1/u`) or compute c,s in float for 16-bit Real; add a half/bf16 jacobisvd regression test. The A1 epsilon fix narrows but does not close the window. Strictly this is beyond element-wise scope, but it's a ~5-line numerics guard with an outsized quality payoff; recommend including. ### A14. Latent scalar-core guard bugs (P2, fold into §5 refactor) - HIP `__bfloat16_raw(unsigned short)` value-converts via `hip_bfloat16(float)` instead of bit-storing (`BFloat16.h:104`) — currently unreachable, silently wrong for any future caller. - `EIGEN_DEVICE_FUNC` push (`EIGEN_GPUCC`) / pop (`EIGEN_CUDACC`) mismatch in BFloat16.h — leaks the redefinition for the rest of every hipcc TU (`BFloat16.h:298-300` vs `:387-389`). - Dead `EIGEN_HAS_NATIVE_BF16` machinery (referenced `BFloat16.h:296,302`, defined nowhere). - `EIGEN_HAS_NATIVE_GPU_FP16` never `#undef`'d (`Half.h:375`); `EIGEN_HAS_BUILTIN_FLOAT16` defined `1` in one place, valueless in the other (`ConfigureVectorization.h:316` vs `:468`). --- ## 5. Workstream B — scalar cores: parity, interop, and the storage refactor ### 5.1 API parity between Half.h and BFloat16.h `bfloat16_impl` has `cbrt, sinh, cosh, asinh, acosh, fmin, fmax, nextafter` (+ mixed int operators); `half_impl` has none of them (#2937 asks for `nextafter` specifically). The generic promote-to-float templates make most of these *work* for half, but with different overload-resolution semantics and no device intrinsics. Add the missing half overloads following the house style (compute in float, round once — §2 item 2 is the invariant); consider `numext::nextafter<half>` bit-twiddling mirroring bf16's. Conversely bf16's non-IEEE `nextafter(-0,+)` → +0 quirk is test-pinned; leave it, document it. ### 5.2 Storage-selection refactor (enables fp8 later; do it once) The 3-way native-storage condition (`EIGEN_HAS_GPU_FP16 || EIGEN_HAS_ARM64_FP16 || EIGEN_HAS_BUILTIN_FLOAT16`) recurs ≥6 times in Half.h; BFloat16.h solves the same problems with a *different* idiom (public `EIGEN_CONSTEXPR` push/undef vs private `_EIGEN_MAYBE_CONSTEXPR`). Consolidate to: one storage-selection traits point per type, a generalized `construct_from_rep_tag` (bf16 lacks it — root enabler of the A2/A14 bug class), and one shared MAYBE_CONSTEXPR macro. While there: condition MAYBE_CONSTEXPR on `EIGEN_HAS_GPU_FP16` **only** — `__builtin_bit_cast` (GCC≥11/Clang≥9/MSVC≥19.27) makes the `__fp16`/`_Float16` branches constexpr-capable, restoring the constexpr that every Apple-silicon/Graviton build silently lost with !2273. **Frozen surface (from the extension-contract audit — refactor must keep these bit-for-bit and name-for-name):** `half_impl::__half_raw` with member `.x` (and its 3 storage variants) + nested `half::__half_raw`; `raw_uint16_to_half` / `raw_half_as_uint16` / `float_to_half_rtne` / `half_to_float`; the `numext::bit_cast` uint16↔half/bf16 specializations; the inheritance chain `half : half_base : __half_raw` (ADL through the base namespace is what makes `numext` math resolve — moving math out of `half_impl` is a silent downstream break); `__bfloat16_raw` with member `.value` (do **not** harmonize the x/value asymmetry — both are load-bearing); `float_to_bfloat16_rtne<false/true>` two-specialization design; implicit `operator float()`; sizeof==2 + trivially-copyable (static-asserted by bit_cast, relied on by TensorRandom/Serializer). Also: Half.h is MPL-2.0+Giesen, BFloat16.h is Apache-2.0/TF — keep the files separate; share mechanism via a third header if needed. Also from the contract audit: compound operators (`-=`, `/=`) are a hard requirement (TriangularSolverMatrix instantiates them); document the minimal custom-scalar contract in `doc/CustomizingEigen_CustomScalar.dox` (currently understates: no mention of numeric_limits, compound ops, dummy_precision). ### 5.3 Native-type interop (new capability, cheap) Today user `_Float16` data cannot round-trip to `Eigen::half` without manual bit_cast unless AVX512FP16/zfh chose that storage. Add always-available **explicit** conversions (ctor + conversion operator) to/from `_Float16` and `__fp16` gated on feature detection *independent of storage choice*; same for `__bf16` (GCC≥13/Clang≥17) on bfloat16; optional `std::float16_t`/`std::bfloat16_t` blocks under `__STDCPP_FLOAT16_T__`/`__STDCPP_BFLOAT16_T__`. Keep uint16 as the ABI; never expose builtin types in public signatures (MSVC). ### 5.4 GPU-facing scalar work - **`std::numeric_limits<Eigen::half>` is uncallable from device code** (members lose constexpr under GPU builds and carry no `EIGEN_DEVICE_FUNC`, `Half.h:287-295`); bf16's works on CUDA but is degraded on HIP by the EIGEN_CONSTEXPR blanking. Annotate with `EIGEN_DEVICE_FUNC`; the hip_bf16.h migration (§6.2) removes the blanking. - **CUDA ≥ 12.2 host-operator coexistence audit**: NVIDIA now ships host-side `__half`/`__nv_bfloat16` arithmetic + implicit integer conversions by default; Eigen's operators/conversions must coexist (the #3077 ambiguity was this class). Add a test TU that includes `cuda_fp16.h` + `cuda_bf16.h` + Eigen together across the CI CUDA matrix; prefer relying on the `__CUDA_NO_HALF_OPERATORS__` macro family (still supported in 13.3) if conflicts appear. - Scalar device intrinsics not yet used: `hrsqrt` (A5), `hrint`/`htrunc`, `hlog2` — small wins, fold into the device-math MR. --- ## 6. Workstream C — GPU: give bfloat16 a real device story; finish fp16 ### 6.1 CUDA bf16 (staged) Today `cuda_bf16.h` appears nowhere; `EIGEN_HAS_CUDA_BF16` is tested at `Eigen/Core:56` but **never defined**; every bf16 device op is software float-emulation even on Hopper; there are no bf16 packets, no `__ldg`, no `__shfl_sync`, no atomics; and the bf16 GPU test is triple-gated dead (§7). Staged adoption, each stage independently shippable: 1. **Wire the header**: include `cuda_bf16.h` next to `cuda_fp16.h` (`ConfigureVectorization.h:576-583`), define `EIGEN_HAS_CUDA_BF16` under `EIGEN_CUDACC` (floor 11.4 > 11.0, unconditional), resolve the `Eigen/Core:53-58` TODO. Add native scalar device arithmetic mirroring `Half.h:383-412` (the header self-emulates below sm_80 — no new arch guards needed beyond what it does internally), device bit-test isnan/isinf, `__ldg`, `__shfl_sync`, `loadConstant`. 2. **Packets**: `__nv_bfloat162` / a `Packet4bf2` mirroring half2/Packet4h2 (arithmetic, math via `h2exp`-family bf16 counterparts where available, else float), `GPU/TypeCasting.h` bf16↔float, `PacketType<bfloat16,GpuDevice>`. 3. **Reductions/atomics** — only after the half-reduction dispatch question (6.3) is resolved, so new kernels aren't dead on arrival. ### 6.2 HIP bf16: migrate off the legacy class The `__bfloat16_raw : public hip_bfloat16` inheritance (`BFloat16.h:101-105`) is why `EIGEN_CONSTEXPR` is blanked for **every hipcc TU** (`BFloat16.h:23-33`). Migrate to `hip_bf16.h`/`__hip_bfloat16` (CUDA-parity API) gated **ROCm ≥ 6.1** (multiple-definition link bug through ≈6.0.2); keep the legacy path for 5.6–6.0. This restores constexpr under hipcc, kills the A14 latent ctor bug, and makes the CUDA/HIP bf16 code paths symmetric — do it before or with 6.1 stage 2 so the packet work lands on both vendors at once. ### 6.3 fp16 device cleanups - **Resolve the dead optimized half-reduction machinery**: host-side dispatch reads `PacketType<half,GpuDevice>` in the host phase where it's the scalar template (`HasAdd=0`) since cc3573ab4 (2021), so `FullReductionKernelHalfFloat` + `atomicReduce(half2/Packet4h2)` + the "needed for GPU half reduction" packet specializations (~500 lines, `TensorReductionGpu.h:123-783`, `GPU/PacketMath.h:1461-1498`) are unreachable; half reductions run the generic path. Verify with a launch trace on a real GPU, then either restore host-visible dispatch (a GpuDevice-specific constant keyed on `CoeffReturnType==half`, independent of packet_traits) or delete. Don't leave it to rot. (Note: the fallback is *correct*; the fp16-accumulation semantics of the optimized kernels are themselves an accuracy question — if restored, consider f32 partials.) - **Use the rounding/transcendental intrinsics**: `h2floor/h2ceil/h2rint/h2trunc` for pfloor/pceil/print/ptrunc (currently ~5-10-op generic emulation), `h2rcp` for preciprocal, `hlog2/h2log2/h2exp2/h2log10` where accuracy allows; keep sin/cos/tanh on float deliberately (comment why). - **Traits hygiene**: set `HasCmp=1` (pcmp_*/pselect exist but the flag is off, forcing scalar fallback for comparison functors); derive `PacketType<half,GpuDevice>` flags from `packet_traits<half>` instead of the hand-copied enum that has already drifted (HasExpm1/HasLog1p/HasPow); A11 semantics fixes. - **Dead-code sweep**: nested `EIGEN_GPU_COMPILE_PHASE` guards inside the device-only block, the 70500 SDK threshold (`Macros.h:653` — equivalent to `>0` at the 11.4 floor), HIP shfl FIXMEs. --- ## 7. Workstream D — CPU backend completeness (element-wise) ### x86 1. **Segment/masked load-store for 16-bit packets** — none exists (`masked_*_available=false` everywhere; `has_packet_segment` only float/double/int); AVX512BW k-mask `vmovdqu16` does it in one instruction. Biggest generic win for cwise tails on AVX512. 2. **AVX512FP16 fidelity gaps** (native path, real-user hardware, zero runtime CI): - `psqrt/prsqrt/preciprocal` use Newton-from-`rsqrt_ph` *unconditionally* (`PacketMathFP16.h:899-948`) while exact `_mm512_sqrt_ph`/`div_ph` exist — breaks the float convention (Newton only under `EIGEN_FAST_MATH`). Regression vs the non-FP16 build, which is provably correctly rounded. - `predux/predux_mul` reduce in fp16 (`_mm512_reduce_add_ph`) — overflow at 65504 and O(lanes) error where non-FP16 builds reduce in f32. Convert-and-reduce-in-f32 for add/mul; keep native min/max. - Enabling AVX512FP16 silently **drops** vectorized erf/ndtri/bessel/sinh/cosh/log10/cbrt for half (traits zeroed "should be implemented in future", `PacketMathFP16.h:56-62`) — instantiate the via-float wrappers at 8h/16h and 2×16h-split for 32h. - `pmadd` gated on `EIGEN_VECTORIZE_FMA` instead of FP16 (FMA-PH is mandatory in AVX512FP16); memory-round-trip ptranspose for 16h/8h sub-blocks. 3. **AVX512BF16**: the single use (`_mm512_cvtneps_pbh`) is GCC-only via a stale `EIGEN_GNUC_STRICT_AT_LEAST(10,1)` guard (`AVX512/PacketMath.h:2658-2684`) — relax to any compiler defining `__AVX512BF16__` (clang has supported the intrinsics and casts for years); consider `_mm512_cvtne2ps_pbh` to halve conversion cost. (`vdpbf16ps` is GEMM-territory → Appendix A.) Note the hw-convert path flushes subnormals — divergence from the soft path worth a comment/test decision. 4. **Housekeeping**: `EIGEN_VECTORIZE_AVS512VL` typo makes AVX512DQ+VL fast casts dead (`AVX/TypeCasting.h:235,275`); `_mm256_srai_epi64` mis-gated on AVX512FP16 instead of AVX512VL (`AVX/PacketMath.h:522-526`); redundant nested FP16 guards; bf16 log-family behind AVX512DQ with a stale "fails test" comment — re-test and either enable or document. 5. **SSE-only**: half/bf16 are fully scalar. Decide explicitly: either a minimal F16C-based Packet4h (F16C implies ≥AVX in practice, so arguably moot) or — recommended — document SSE-only 16-bit as scalar-by-design and rely on the generic path. No silent middle ground. ### NEON / SVE / SME 1. **hf transcendentals**: `pexp/plog/psin/pcos` trait-disabled (`NEON/PacketMath.h:5417-5420`) so half `.exp()` runs one-scalar-at-a-time on the most widely deployed native-fp16 platform; the widen-to-f32 pattern is already in-tree (`ptanh`, `UnaryFunctors.h:24-52`). Add the full `F16_PACKET_FUNCTION` set at 4hf/8hf widths and flip the traits. 2. **Packet4bf flags**: `HasSqrt/HasRsqrt` (+`HasCbrt`) are 0 although implementations are instantiated — one-line trait fix. 3. **Casts**: no vectorized bf16↔float on NEON (verified LinearTraversal vs the half version's vectorized cast) — `type_casting_traits` + `pcast` reusing `F32ToBf16/Bf16ToF32` are nearly free. Same for RVV10 (all four 16↔32 pairs missing; its declared `type_casting_traits<_Float16,int16>` are keyed on the wrong type and dead). 4. **Guard hygiene**: make `EIGEN_HAS_ARM64_FP16_SCALAR_ARITHMETIC` a 0/1 value-tested macro like its VECTOR sibling; document/rename `EIGEN_HAS_ARM64_FP16` as a storage/format macro (it lacks the GPU-phase exclusion its siblings have — part of the A2 fix); `template<>` missing on Packet8hf pcmp in the fallback macro; duplicate `eigen_packet_wrapper` ID 19 (Packet4bf vs MSVC Packet1d); fp16 native-branch predux accumulates in fp16 while the *fallback* branch reduces in f32 — an accuracy inversion between build modes; consider f32 predux-add in the native branch for consistency. 5. **SVE**: zero fp16/bf16 despite fp16 being architecturally mandatory in SVE — add `PacketXh` (svfloat16_t) mirroring PacketXf + fp16↔f32 casts; bf16 packet as NEON-style widen/narrow later. 6. **SME**: the SME build includes NEON PacketMath but **replaces** the NEON GEBP header (`Eigen/Core:443-444`), silently dropping the NEON fp16 (and float/double) lane-FMA GEBP traits — include both (SME traits pin only float/float). Widening FMOPA/BFMOPA kernels are follow-up-project material. ### RVV10 Flip `HasExp/HasLog` (+Expm1/Log1p/Log2) for fp16 — the via-float kernels are already compiled, just trait-disabled (`PacketMathFP16.h:72-73` vs `:839-859`); add the casts (above); A3 mask fix. (The bf16 GEBP per-madd renarrow is follow-up material, Appendix A.) ### SYCL `PacketWrapper<cl_half8, 8>` has no device specialization → vectorized half reductions can't compile if instantiated (`InteropHeaders.h:129-218`, `TensorReductionSycl.h:168`) — add the size-8 specialization or clamp half's advertised packet size; no bf16 at all (note as gap; SYCL has `ext::oneapi::bfloat16` if demand appears). No SYCL CI exists — flag any investment here as untestable today. ### Non-goals (document, don't build) HVX/LSX/MSA/ZVector 16-bit packets (no hardware story or no demand), AArch32 fp16 vectorization, `std::complex<half/bf16>` (formally unspecified for non-builtin T; zero usage in-tree — state "unsupported" in docs, or add one smoke test if we want to keep it de-facto working; the A5 hypot fix helps it either way). --- ## 8. Workstream E — tests and CI ### Tests (ordered by leverage) 1. **Exhaustive 2¹⁶ ULP tests** (the single highest-leverage addition): for each unary math function and both 16-bit types, sweep all 65536 bit patterns against a double-computed reference rounded to the type, with an explicit per-function max-ULP budget; plus an exhaustive half→float exactness check and a structured float→half RNE check (all rounding-boundary neighborhoods). Cheap (65k × ~40 fns ≈ milliseconds), definitive, and it converts the A5-class bugs from "measured by an audit" to "caught by CI". Wire into the smoke list. 2. **Re-enable the skipped 16-bit conformance blocks**: psin/pcos special-values + large-arg (scoped to the type's range) for half/bf16 (`packetmath.cpp:1190-1235`, existing TODO with my name on it); bf16 plog edge cases (`:1131-1160`); rsqrt scalar test after A5. 3. **Mask-exactness in packetmath**: `CHECK_CWISE2_MASK` full-bit comparison (catches A3 class on every backend). 4. **GPU**: host-side raw-bit round-trip test under `EIGEN_CUDACC` (A2); half/bf16 subtests in `gpu_basic.cu` (currently zero non-Tensor device 16-bit coverage); 16-bit uniform+normal in `tensor_random_gpu.cu`; register `tensor_of_bfloat16_gpu` in the **CUDA** test block (today HIP-only, and 8/9 subtests additionally dead behind the never-defined `EIGEN_HAS_CUDA_BF16` — fix both with §6.1 stage 1; drop the stale reduction-exclusion comment); device-callable numeric_limits check; the CUDA≥12.2 header-coexistence TU (§5.4). 5. **Instantiation coverage**: one 16-bit CALL_SUBTEST per decomposition suite (residual-norm at dummy_precision; documents which solvers are *supposed* to work at 16-bit); redux rows; serializer row; sparse compile+smoke set (A9/A12); tensor_contraction/thread_pool 16-bit instantiations (currently zero — the TF-hot path). 6. **Scalar-file symmetry**: bf16 file lacks denormal/fmod/increment tests that half has; half file lacks the round-trip/nextafter bit-tests bf16 has — cross-port using the `VERIFY_*_BITS_EQUAL` style. ### CI 1. **AVX512FP16 runtime coverage via Intel SDE** (`sde64 -spr`) — the path is compile-only today, demand is proven by breakage reports, and the build job already documents a GCC -O3 miscompile forced to -O2. Even a packetmath_13-only SDE job covers most of the risk. Prerequisite for any further AVX512FP16 work. 2. **F16C on mainstream jobs**: `-DEIGEN_TEST_F16C=on` on the existing gcc/clang AVX2 jobs (probing confirmed gcc `-mavx512dq` doesn't even define `__F16C__`; today the hardware conversion path runs only in two nightly clang jobs). 3. **POWER10 MMA**: one ppc64le job with `-mcpu=power10` under qemu (the 745-line bf16 MMA stack currently never even compiles in CI). 4. **RVV bf16**: add `_zvfbfmin_zvfbfwma` to the riscv64 march string (backend never built in CI; A3 sat invisible for exactly this reason). Plan to drop `allow_failure` on the riscv jobs when the runner is trusted. 5. **aarch64 + CUDA compile-only job** (A2 class; Grace/Jetson is a real deployment target). 6. **GPU job policy**: all GPU test jobs are `allow_failure:true` and run only on schedule/label — master merges run zero GPU tests. Once stable, move to the exit-42 retry-whitelist the CPU jobs use; consider a periodic sm_90 runner (Hopper is compile-only today). --- ## 9. Suggested MR sequence Small, independently-landable MRs; per repo convention each lands with its tests (and benchmarks where perf-relevant). Order minimizes rebase pain against !2700 (packet kernels) and keeps the epsilon/SparseQR coupling. | # | MR | Contents | Depends on | |---|---|---|---| | 1 | `Core: fix NumTraits<half>::epsilon` | A1 + A12 SparseQR float-threshold + EigenSolver/SVD re-check + changelog note | — (label `all-tests`) | | 2 | `Core: make half raw-bit interop storage-aware` | A2 (4 sites) + host round-trip test + ARM guard hygiene (D-NEON item 4) | — | | 3 | `Core: 16-bit scalar math accuracy` | A5 (log2/rsqrt/hypot) + A6 bf16 bit-classification + exhaustive-ULP test harness (E1) | after !2700 lands/rebases | | 4 | `Core: half/bf16 scalar API parity` | §5.1 (cbrt…nextafter, #2937) + A10 hash + scalar test symmetry (E6) | — | | 5 | `Sparse: 16-bit literal fixes + smoke tests` | A9 + A12 isApprox/BiCGSTAB docs-or-guards + sparse smoke set | 1 | | 6 | `Tensor: fix 16-bit normal RNG` | A4 + finiteness/range tests + GPU RNG tests | — | | 7 | `Core: dense 16-bit reduction policy` | A8 (mean divisor; optional LeafSize port — benchmark) + blueNorm A7 | — | | 8 | `RVV10: fix bf16 masks; enable fp16 math traits` | A3 + D-RVV + CI march string (E-CI 4) + mask-exactness test (E3) | — | | 9 | `NEON: complete fp16/bf16 element-wise` | D-NEON 1-3 (+SME GEBP include fix) | 3 | | 10 | `AVX512FP16: exactness + completeness` + SDE CI job | D-x86 2, E-CI 1 | E-CI 1 first | | 11 | `GPU: CUDA bf16 stage 1` | §6.1(1) + test gate fixes (E4) + coexistence TU | — | | 12 | `GPU: HIP bf16 migration (ROCm≥6.1)` | §6.2 + A14 HIP items | 11 (symmetry) | | 13 | `GPU: bf16 packets (stage 2)` | §6.1(2) | 11, 12 | | 14 | `GPU: fp16 cleanup` | §6.3 (dead reductions decision, intrinsics, traits, A11) | !2681 coordination | | 15 | `Core: scalar storage refactor + native interop` | §5.2 + §5.3 + A14 guard cleanups + constexpr restoration | 1-4 landed; coordinate with !2465 (C++17) | | 16 | `Jacobi: guard 2x2 SVD rotation overflow` | A13 + jacobisvd 16-bit test | 1 | | 17 | Docs | 16-bit scalar user page: backend matrix, accumulation semantics, which solvers work, mixed-precision guidance, custom-scalar contract update | most of the above | | — | CI-only MRs | F16C flags, POWER10 qemu, aarch64-CUDA compile, GPU allow_failure policy | anytime | **Hardware-local GEMM wave (Appendix A.1, in scope):** runs after the correctness MRs above, gated on its testing prerequisite. Sequence: (i) k-scaling GEMM accuracy test + 16-bit GEMM benchmarks; (ii) RVV10 bf16 accumulator fix (smallest); (iii) OpenBLAS sbgemm correctness fix, then sbgemv (#2936); (iv) NEON BFDOT/BFMMLA bf16 kernel; (v) AVX512BF16 `dpbf16_ps` kernel; (vi) SME FMOPA/BFMOPA; (vii) GPU / ThreadPool f32-accumulating contraction kernels. Each is a per-backend kernel override (AltiVec pattern) — no generic-GEBP changes, so they don't collide with the deferred A.2 redesign or #3024. Coordination rules: rebase packet-kernel MRs (3, 8, 9, 10, 14) onto !2700's `EIGEN_FAST_MATH_CONSTANT_BARRIER` idiom; anything touching AltiVec/RVV shared code → ChipKerchner review; Tensor-visible numerics changes (6, 7, 13) → call out TF golden-value impact in the MR description; avoid C++14 contortions that !2465 will delete. --- ## Appendix A — 16-bit matrix products & wider accumulators *Scoping decision (2026-07-09, refined): Eigen's generic GEBP framework pins the accumulator to the result scalar (`gebp_traits` sets `AccPacket=ResPacket` = packet-of-ResScalar, `GeneralBlockPanelKernel.h:486,524-528`), which makes bf16/fp16 GEMM of very little value beyond tiny k — O(k·eps₁₆) error growth, plus fp16 mid-accumulation overflow at |x|>65504, plus (on x86) a ~5-15-uop convert tax per FMA. Supporting a wider accumulator **generically** requires a significant GEBP redesign and is deferred (A.2). But wherever hardware-accelerated bf16/fp16 with wider accumulation can be exploited **locally** — by overriding a whole `gebp_kernel`/dedicated kernel for one scalar type, the pattern AltiVec MMA already proves in-tree — we should, now (A.1).* **Shared context.** Hardware direction is unanimous: `VDPBF16PS` (x86), `BFDOT/BFMMLA` (ARM v8.6), `Zvfbfwma` (RISC-V), SME `FMOPA/BFMOPA`, GPU tensor ops — all widen into f32; no ISA offers non-widening bf16 arithmetic. Industry default (oneDNN "strict", OpenBLAS sbgemm, cuBLAS) is f32 accumulation with 16-bit I/O; reduced-precision accumulation is opt-in. The in-tree model is the AltiVec bf16 stack: `gebp_kernel<bfloat16>` override → MMA kernel with f32 `__vector_quad` accumulators, result staged in f32 and rounded once per kc-panel (`MatrixProduct.h:3662`, `MatrixProductMMAbfloat16.h:326-365`); its pre-P10 VSX fallback proves the convert-panels-to-f32 variant. ### A.1 In scope now — hardware-local wide accumulation (the AltiVec pattern) Each of these overrides a kernel for one scalar type on one backend; none threads a wider `AccPacket` through the generic traits. All store to the 16-bit result once per panel, so accuracy becomes O(k·eps_f32) + one final rounding. **Prerequisite for all of them** (land first): a k-scaling accuracy test against an f64 reference with a (c₁+c₂√k)·eps tolerance model — current product tests compare against a same-type reference at dummy_precision and cannot see accumulator-precision changes — plus 16-bit GEMM benchmarks per repo convention. 1. **RVV10 bf16 GEBP accumulator** — the cheapest fix in the whole area: the dedicated micro-kernel *already uses* the widening f32 FMA (`vfwmaccbf16`) and then explicitly rounds back to bf16 after every madd (`RVV10/GeneralBlockPanelKernel.h:203-227`). Keep the accumulator in f32 across k inside the existing kernel, narrow once in `acc()`. 2. **NEON bf16 kernel** (ARMv8.6, `__ARM_FEATURE_BF16_VECTOR_ARITHMETIC`): dedicated bf16 `gebp_kernel`/traits override using `BFDOT`/`BFMMLA` into f32 accumulators (Graviton3, Neoverse V1/N2, Apple M2+; acknowledged TODO at `NEON/PacketMath.h:4627-4630`). Model it on `MatrixProductMMAbfloat16.h`. 3. **AVX512BF16 kernel**: `_mm512_dpbf16_ps`-based override with f32 accumulators (2 bf16 MACs/lane/instr, ~2× fp32 peak; Cooper Lake, SPR+, Zen4/5). Companion conversion work: relax the GCC-only `cvtneps_pbh` guard, use `cvtne2ps_pbh` for pack/store. 4. **SME widening kernels**: extend the existing fp32 ZA-tile GEMM (`arch/SME/GeneralBlockPanelKernel.h`) with `FMOPA` (fp16→f32) and `BFMOPA` (bf16→f32) variants — the packers/predication machinery generalizes. 5. **GPU half contraction**: an f32-accumulating kernel (or converting mappers into the float kernel) to replace the generic shared-memory kernel that keeps half registers/shmem (`TensorContractionGpu.h:41-57,1364`) — matches cuBLAS HGEMM-with-f32-compute semantics. Same option for bf16 once GPU bf16 packets exist (§6.1). 6. **ThreadPool contractions**: an f32-accumulating 16-bit `TensorContractionKernel` specialization — the clean TF-frozen seam (its template signature, `TensorContractionInputMapper`, and `TensorContractionBlocking` signatures must not change). Flag that bf16/half contraction outputs change bitwise (improve) — TF golden values. 7. **OpenBLAS backend (#2936)**: the merged sbgemm path (!1909, `EIGEN_USE_OPENBLAS_BFLOAT16`) has two correctness bugs — `fbeta=1` against an *uninitialized* temp (`GeneralMatrixMatrix_BLAS.h:149,189-191`) and overwrite-instead-of-accumulate violating the `res += αAB` contract (`:196`) — and is untested/undocumented/unwired. Fix + test first, then add the requested sbgemv following the same pattern. 8. **Local correctness fix that rides along**: `TensorContractionSubMapper::loadPacketPartial` ignores its count argument (`TensorContractionMapper.h:412-418`) — OOB hazard for the Power bf16 pack path through tensor mappers; implement real partial semantics (kernel-local fix, no framework change). Where no widening hardware exists (x86 half without AVX512FP16, pre-P10 bf16 already handled), the *generic* remedy is convert-on-pack — that stays in A.2 because it needs shared pack/traits infrastructure. ### A.2 Deferred — the generic GEBP framework redesign 1. **`AccPacket`/`ResPacket` decoupling in generic `gebp_traits`** with a converting `acc()` epilogue — benefits every backend without dedicated kernels, but requires redesigning the register-blocking state machine; coordinate with #3024 (DoublePacket/predux_half overloadability, ChipKerchner's stated refactor intent). 2. **Convert-on-pack infrastructure**: f32-emitting `gemm_pack_lhs/rhs` for 16-bit inputs + f32 kernels with converting store — the generic fallback for backends without widening ops (in-tree precedent: AltiVec's VSX `gemmbfloat16`, `MatrixProduct.h:3127`). 3. **Mixed-precision user API** (bf16×bf16→f32 as a first-class product/contraction result): `ScalarBinaryOpTraits` specializations + contraction output-type override — today `Matrix<half>*Matrix<float>` doesn't compile and `(A*B).cast<float>()` is a documentation trap. API design question as much as kernel work. 4. **Cost-model integration**: contraction cost model assumes native-FMA throughput for 16-bit (under-provisions threads, over-coarsens grains) — fix via a gebp-level emulation/throughput trait, *not* NumTraits AddCost (which feeds Core unrolling everywhere). Related 16-bit accumulation-policy items best decided with the framework: shard-by-inner-dim partial buffers typed on an accumulator scalar, TensorConvolution per-tap accumulation, Tensor scan/cumsum. ## Appendix B — deliberate decisions to record (docs material) - Via-float per-op rounding **is** the accuracy contract for 16-bit element-wise ops (provably exact per op); native-fp16 chains (AVX512FP16, NEON+fp16, RVV zvfh) match scalar semantics per op but differ in reductions — the AVX512FP16 predux items in §7 exist to close that gap; product-accumulation semantics are Appendix A territory (A.1 hardware-local now, A.2 generic framework deferred). - `dummy_precision` = 10.2×eps (half) / 6.4×eps (bf16) is much tighter relative to eps than float's 84×; it's simultaneously too tight for accumulation-heavy checks and too loose to catch 1-ulp regressions. Resolution: keep the values (test-ecosystem-wide impact), add the exhaustive-ULP harness (E1) for exactness, and document the ~6-10-ulp meaning of `isApprox` at 16 bits. - bf16 `is_iec559=true` (`BFloat16.h:188`) is technically wrong (not an IEC interchange format; `std::bfloat16_t` implementations say false) — low priority, breaking-change-flavored; revisit at the C++17 bump. - Tensor uniform RNG is [0,1) step 2⁻¹⁰/2⁻⁷; core `random()` is [-1,1) step 2⁻⁹/2⁻⁶ — divergent by design; document, don't unify. - AltiVec bf16 subnormal handling is compiled out by default (`SUPPORT_BF16_SUBNORMALS` knob) — documented, deliberate; leave.