Ozaki 2.5: Engineering the Deconstruction Path of fp64-Emulated Dense Matrix Multiplication on FP8 Tensor Cores∗ Satoshi Matsuoka†
arXiv:2609.09095v1 [cs.MS] 8 Sep 2026
Director, RIKEN Center for Computational Science (R-CCS) Kobe, Hyogo, Japan
Draft 27 of September 3, 2026
Abstract FP8 Ozaki II emulates fp64 matrix multiplication with a fixed schedule of low-precision tensor-core products over a CRT residue system. Converting the fp64 operands into residue planes—the deconstruction charged by the fourth term of the Tensor–Memory Equilibrium (TME) model of Part 1 [11], identified in the NVIDIA technical review [2]—consumes integerpipe and memory resources before any tensor instruction issues. This paper engineers that path. It is a companion to Part 1 but an independent contribution: Part 1 argues a platformcoverage thesis and prices deconstruction as one cost term among four, never opening the conversion path itself, whereas every technical result below—the split closures, the twolimb kernels, the modulus/encoding codesign, the route menu and its crossover algebra, the reconstruction-aware engine, and the closed-form reach floor with its hardware asks— appears only here. We make four contributions, all conditional reduced-model projections pending the measurement campaign we specify; no new kernel is measured here. (1) A deconstruction-aware performance model, with the honest ceiling stated up front. Deconstruction competes with the low-precision stream through the harmonic mean n̄ = 2mn/(m+n) of the output dimensions: below a size crossover the emulated rate rises with n̄ (the branch through NVIDIA’s preliminary ∼ 200-TFLOPS “Emulated DGEMM” at modelinferred n̄ ≈ 500–900) and approaches the raw arithmetic roof Pfp8 /(3r+1) (≈ 473 TFLOPS at Rubin’s 17.5-PFLOPS rate, r=12) only within a single thread-block cluster (a ≤2562 output). A larger output is re-split on the fly, once per cluster, holding the achieved rate at the deconstruction-λ floor, ≈ 235 TFLOPS on Rubin for cluster-aligned large squares—half that 473-TFLOPS roof, the one-half being a derived ratio of three design integers (cluster reach, pipe-provisioning ratio, plane-formation cost), not a fitted efficiency—and the regime a DGEMM benchmark such as HPL runs in (≈1.9 EFLOPS fp64 on a 10,000-GPU cluster). That figure is route- and shape-conditional and we state the condition up front: it is the codesigned hybrid set E at NB ≈ 1024, the upper end of the cited 892–1024 GPU-HPL panel-width range; the published set S, which carries the round-to-nearest Ozaki II theorem, floors at ≈ 182 TFLOPS (0.38 of its 473-TFLOPS roof, ≈1.4–1.5 EFLOPS). Ragged, non-aligned sizes dip further, to a disclosed rigid-schedule worst case of ≈0.36 of that same ∗ “Ozaki 2.5” names an implementation-level proposal: the mathematics of Ozaki Scheme II is unchanged on the published modulus set (accuracy inheritance is proved for the round-to-nearest variant; the implemented truncation rule is a stated obligation, and codesigned variants carry per-set proof obligations, Appendix A), and the terminology will be coordinated with the Ozaki-II/FP8-Ozaki-II authors. All performance numbers in this paper are conditional reduced-model projections; no new kernel is measured here. A companion paper, “FP8 is All You Need” (Part 1, arXiv:2606.06510), covers the shared tensor–memory equilibrium model and hardware co-design proposal. † Correspondence: [email protected]
1
roof—the only minimum we construct. A ragged-edge-aware (predicated) schedule is modeled to recover ≈0.54 there; we log that as a Part-3 validation target, not a second guaranteed floor. Crossover and floor are one curve, drawn together (Fig. 3), not a late reveal. The positive corollary is that the floor is a large-square phenomenon: the tall/skinny and small-batch shapes that dominate real solvers (block-Krylov, batched GEMMs, panel factorisations) re-split their big operand only once (λA =1, bounded λeff ≈1.25–1.5) and stay near the crossover—only lightly clipped (0.89–0.94×), never at the floor—so Ozaki 2.5 is already worth ≈1.6–1.9× (≈2× on block-Krylov shapes) over simple deconstruction on Rubin today, with no hardware change (Table 11). (2) The Ozaki 2.5 method. Leaving the Ozaki II reconstruction framework unchanged, we engineer the deconstruction path around the actual 9–11-bit FP8 modulus set: convert-once residue-plane workspaces, an exact two-limb constant-reduction GEMM on integer tensor pipes (or its pure-SIMT dp4a realisation), a ≈6.3-instruction SIMT residual, and conversion pipelined behind the MMA stream—moving the Rubin crossover from ≈1,211 to ≈480–730 and modeled throughput at n̄=512 from ≈200 to 364–438 TFLOPS (a conditional 1.8–2.2× convert-once envelope). (3) Modulus/encoding codesign. A script-checked study shows the modulus set is itself a performance parameter: an all-byte system (15 coprime moduli ≤ 256) legalises a onepass reduction and is projected faster below n̄ ≈ 540 despite a 20% lower roof; a hybrid (squares ≤ 332 plus a byte tail) concedes only 7.5% asymptotically while dominating below n̄ ≈ 620; and two exhaustive supply bounds (squares cap at 94.1 bits, sub-64 systems at 89.9, both short of the 111.8 required) close the space—motivating a regime-switched modulus dispatch. (4) Potential hardware co-design. Because the floor is the closed-form R Pint /(cq r)— a per-element deconstruction load cq r/reach—it names its own escape: we specify cluster reach, TMEM capacity, and L2 service bandwidth for dense DGEMM (with a target table), and the deconstruction-cost (cq ) datapath of Part 1—minimally an in-flight-convert copy engine, which also un-binds the conversion-bound sparse kernels—each a concrete, measurable ask that lifts the emulated rate toward the arithmetic roof, none yet present in silicon. None of this is Rubin-specific: instantiated at Blackwell Ultra (GB300) rates the same machinery yields a 135-TFLOPS roof already at its own floor (roof-bound), a ∼166-TOPS residual-int8 cap that reshuffles the route ranking (the all-byte set leading via pure-SIMT dp4a), and, against a ≈1.4-TFLOPS native fp64 pipe, every modeled route above native for n̄ ≥ 32—so the method is present-tense on shipping hardware. Measured GEMM/SYRK call-shape traces across four application classes (LOBPCG, multifrontal LU, CCSD, blocked LAPACK) ground the shape analysis. Every rate is flagged achieved today—meaning modeled on shipping silicon with no hardware change, a reduced-model projection, not a measured kernel—or convert-once envelope (clipping-limited) where it appears; four falsifiable predictions with a decisive size sweep, baselines (GEMMul8 [26], fused-kernel work [10]), and negative controls define the test, and all scripts, traces, and the parameter file generating every figure and table are supplied with this submission as a reproduction artifact (arXiv ancillary files); an archival Zenodo DOI and commit hash will be added in the first arXiv revision.
Keywords: FP64 emulation; Ozaki Scheme II; Ozaki 2.5; FP8 tensor cores; deconstruction cost; dense matrix multiplication; hardware co-design; Tensor–Memory Equilibrium model; NVIDIA Blackwell (GB300); NVIDIA Rubin.
1
Introduction
Double-precision dense matrix multiplication is entering its emulated era. On NVIDIA’s Rubin generation, native fp64 matrix throughput is no longer the headline; instead the official specifications list “Emulated DGEMM” as a first-class column—∼ 200 TFLOPS per GPU [17, 9]—a preliminary, “up to” specification whose algorithm, mode, and benchmark dimensions are not published; the natural candidate is the Ozaki II scheme [23] on fp8 tensor cores [27, 13], which 2
NVIDIA’s public material describes as the direction of its emulation work [14]. NVIDIA thus positions emulated DGEMM as the primary high-throughput fp64 matrix path, and its projected rate is a product specification. That specification invites a comparison the companion Part 1 paper [11] makes precise. At the official 17.5-PFLOPS dense fp8 rate of the DGX Rubin NVL8 specification [18], the analytic ceiling of fp8 Ozaki II at r=12 moduli—3r+1 = 37 fp8 MMA passes per fp64 product—is Pfp8 17,500 = ≈ 473 TFLOPS, 3r + 1 37 so the announced specification is ≈ 42% of the raw arithmetic quotient. That 473-TFLOPS quotient is a peak-operation roof, not an attainable ceiling—tile efficiency, launches, reconstruction, and resource contention all take their share, and B200 calibration suggests a sustained-fp8 efficiency well below unity at large sizes [27]—so a delivered efficiency of ∼40% (a ∼60% shortfall) could be unremarkable. But part of the margin could also have a specific, identifiable, and largely removable cause, with a size-dependent signature that the alternative lacks. The answer, up front. The removable part is the deconstruction cost, and honesty requires stating its ceiling before the constructions that chip at it. For a small output the emulated rate rises with n̄ and can approach the roof; but a large dense output—the regime a DGEMM benchmark such as HPL runs in—spans many thread-block clusters, on-the-fly conversion repeats once per cluster, and the achieved rate settles at the deconstruction-λ floor— ≈235 TFLOPS on Rubin for cluster-aligned sizes, half the 473-TFLOP roof, with ragged sizes dipping to a disclosed rigid-schedule worst case (Fig. 3, drawn with that floor so it is visible early, with the crossover, not only at the end). Ozaki 2.5 thus does not deliver the raw roof for large DGEMM on today’s silicon; it delivers about half—still ≈8× native fp64, consistent with NVIDIA’s ∼200-TFLOP figure—and, because the floor is set by the load cq r/reach, a co-design coordinate we lift back to the roof (§5.2). The equally important positive is that the floor is a large-square phenomenon: the tall/skinny and small-batch matrices that dominate real solvers (block-Krylov, batched GEMMs, panel factorisations) re-split their big operand only once (λA =1; bounded λeff ≈1.25–1.5), stay near the crossover (a mild 0.89–0.94×, not the floor), and are already worth ≈1.6–1.9× (≈2× on block-Krylov shapes) over simple deconstruction with no hardware change (Table 11). Every performance number below is flagged achieved today or convert-once envelope (clipping-limited) where it appears—where achieved today means modeled to be attainable on today’s silicon with no hardware change, a reduced-model projection like every rate here, not a measured kernel. The rival explanation for the margin—a size-independent efficiency artefact—carries an orthogonal signature that a single size sweep separates.
1.1
Prior work
The Ozaki scheme writes a high-precision product as a sum of low-precision inner products over a splitting/residue system; the error-free int8 slicing of the Ozaki I lineage [22] and the CRT-based Scheme II of Ozaki et al. [23] are the two realisations relevant here. Its use for fp64 emulation on fp8 tensor cores has been developed and measured by Uchino et al. [27] and Mukunoki et al. [13], with accuracy guarantees proved under round-to-nearest hypotheses [24, 25]; GEMMul8 [26] and fused-kernel emulation [10] are the public baselines. NVIDIA’s material positions Ozaki-style emulation as the direction of its “Emulated DGEMM” path [14], and its technical review [2] first isolated the per-input deconstruction cost that the companion Part 1 [11] folds into the TME performance model as a fourth term. This body of work establishes the scheme, proves its accuracy, and reports achieved throughput. What it does not do—and what this paper adds—is engineer the deconstruction path around the actual FP8 modulus set (the accuracy contract inherits only for the round-to-nearest variant; the shipping truncation rule and the codesigned sets carry stated obligations, Appendix A), treat the modulus set as a runtime performance variable, and localise the large-DGEMM ceiling as the closed-form R Pint /(cq r) 3
floor that names its own hardware escape. The technical recap the rest of the paper builds on—the Ozaki schemes and the int8→fp8 substrate transition that created the deconstruction cost—is §2.
1.2
Contributions and roadmap
The name is chosen deliberately: the method adds no new mathematics to Ozaki II—the moduli and the reconstruction framework are untouched—and stops short of the hardware datapaths whose evaluation belongs to follow-up co-design work. It is the missing half-step, Ozaki II with the conversion path engineered as deliberately as the multiplication path always has been, an engineering need that did not arise before the substrate transition. Concretely: • A deconstruction-aware model with an honest ceiling (§3–§4). We recall the fourterm TME model of Part 1 and show that its fourth term—the per-input deconstruction cost cq , counted at cq =16 instructions per modulus per element on the shipping path [2] (Appendix B)—gates dense GEMM behind a size crossover n∗ = cq rPfp8 /(αPint ) ≈ 1,211 on Rubin, and that a large multi-cluster output is pinned at the deconstruction-λ floor (0.50 of the 473-TFLOPS roof—a derived ratio of three design integers, Eq. (15), not a fitted efficiency). The reduced model passes through the published 200-TFLOPS figure at model-inferred n̄ ≈ 500–900; we state this as a falsifiable hypothesis, not a finding. • The Ozaki 2.5 method (§5). Ozaki II unchanged in its mathematics, with the deconstruction path engineered around the actual modulus set (convert-once workspaces, an exact two-limb integer-tensor reduction, a ≈6.3-instruction SIMT residual, pipelined conversion; Algorithm 1), moving the Rubin crossover from ≈1,211 to ≈480–730. • Modulus/encoding codesign (§5.1). A script-checked study treating the modulus set as a design variable: an all-byte system that legalises the one-pass reduction, a hybrid that concedes only 7.5% asymptotically, and two exhaustive supply bounds—motivating a regime-switched modulus dispatch. • Potential hardware co-design (§5.2). Because the floor is the closed-form R Pint /(cq r), it names its escape: cluster reach, TMEM capacity, and L2 bandwidth for dense DGEMM, and a cq -removing copy-engine datapath that also un-binds the sparse kernels—each a concrete, measurable ask. Relation to Part 1, and why this is a separate paper. Part 1 [11] argues a platform thesis: that an FP8-dominated datapath, priced by the four-term TME model, can serve the scientific dwarfs at large. Ozaki emulation enters there as one workload among many, and the deconstruction term enters as a cost to be priced, not a path to be engineered. This paper takes the opposite cut: it holds the workload fixed—fp64 GEMM emulation—and engineers the single term Part 1 could only price. Every technical result below is new to this paper and appears in no FP8-platform or Ozaki-scheme publication we know of: the carry-corrected redundant split that closes the six nonsquare S-tail moduli 487–511 in e4m3 (Lemma 1); the exact two-limb byte-identity reduction kernels and their pure-SIMT dp4a realisation; the modulus/encoding codesign study—sets E/A/D with per-set proof obligations (Appendix A)—and its two exhaustive supply bounds; the route menu with per-route numerical contracts and the k≈6777 route-crossover algebra; the reconstruction-aware per-record engine with Nacc accounting, grounded by measured LD_PRELOAD call-shape traces; and the closed-form reach floor of Eq. (14) with its quantified hardware asks. From Part 1 we import exactly two things, both by citation: the TME cost model (§3 recalls it) and the memo-derived constants of Appendix B. Nothing here duplicates Part 1’s coverage analysis, and nothing there anticipates the modulus,
4
kernel, or floor results; the two papers share a cost model the way two instruction-set studies share an ISA manual. Every projection is instantiated for the Blackwell generation (in particular GB300) as well as Rubin; because GB300 pairs a ≈1.4-TFLOPS native fp64 pipe with a 135-TFLOPS emulated roof, the method is a present-tense proposition on shipping hardware, not one that waits for Rubin. We close (§6–§7) with measured call-shape traces that ground the shape analysis and four falsifiable predictions whose decisive experiment is a DGEMM size sweep, with a slicingbased int8 kernel on B200 as the control arm—identifying the responsible term whether the branch hypothesis survives measurement or falls to it, in the same two-way spirit as the TME model itself [11]. Because every result below carries one of a small set of epistemic labels, we fix them here in Table 1; each results-table caption repeats its label, and the detailed ledger is in §7. Table 1: Claim-status labels used throughout the paper. label
what carries it
proved algebraically
the harmonic-mean crossover (reduced model); the two-limb byte identity; the square-modulus and Karatsuba plane identities and the accumulator floor they imply; the e4m3 layout envelope (m ≤ 321 nonsquare under a canonical split, m ≤ 577 under a redundant one, m ≤ 1089 square); the carry-corrected three-plane layout of Lemma 1 and its plane bounds; the exact-accumulation bound Kslab max |term| ≤ 224 ; the operand-bandwidth ratio pstore /α and its dsm counterpart (1−1/c) pstore /α supply bounds; coprimality and exact CRT products; centred digit maps and e4m3 representability of all A/E/D planes and the six square S moduli. The six nonsquare S-tail moduli (487–511) are closed: no canonical balanced split places them in e4m3 at any base, but the carry-corrected redundant split of Lemma 1 does, verified exhaustively over all m2 ordered centred pairs per modulus (verify/stail_layout.py) the fp8/int8/fp64 platform rates and the ∼200-TFLOPS emulatedDGEMM figure cq =16 and the Pint ≈75-TOPS normalisation (Appendix B) every route knee, envelope, bracket, and composite speedup (Tables 3–11) the LD_PRELOAD call-shape traces (geometry only)—and nothing else
script-checked
vendor-announced memo-derived modeled projection measured
2
Background: the Ozaki Schemes and the Substrate Transition
From error-free splitting to modular arithmetic. The original Ozaki scheme [22]— Ozaki I in the present numbering—computes a high-precision matrix product by error-free splitting: each operand is split into a short sum of lower-precision slices such that every pairwise slice product is exact in the target arithmetic, and the exact partial products are summed. Its cost grows quadratically in the slice count, which is what makes it attractive on int8 (few wide slices) and prohibitive on fp8 (many narrow ones). Ozaki Scheme II [23] replaces the splitting by residue arithmetic: one exact integer product is computed per modulus of a Chinese Remainder Theorem (CRT) system, so the cost grows linearly in the number of moduli. Because the integer core is exact, the scheme delivers componentwise fp64-grade accuracy—the “Grade A” contract in Part 1’s terminology, i.e., fp64-equivalent componentwise error bounds— with the error analysis and the automatic precision-escalation machinery (ESC/ADP) supplied by [24, 27, 25]. The scheme in equations. Let C = AB with A ∈ Rm×k , B ∈ Rk×n in fp64. (i) Scale to integers. Choose per-row exponents σ1 ..σm for A and per-column exponents τ1 ..τn for B,
5
keeping t significant bits, and set A′ = diag(2σ ) A ,
B ′ = B diag(2τ ) ,
C ≈ diag(2−σ ) (A′ B ′ ) diag(2−τ ),
(1)
where ⌊·⌉ denotes the integer quantiser selected by the named variant (Table 13): round-tonearest for Ozaki-II-RN, under which the cited error theorem [24] holds within its hypotheses; truncation toward zero for Ozaki-II-TZ, the rule the shipping FP8 implementation uses [27] and the one our projections assume. The round-to-nearest theorem is not claimed for the truncation rule: TZ does not automatically inherit the RN bound, and TZ-grade accuracy remains an explicit validation obligation rather than an established result. The integer product C ′ = A′ B ′ is computed exactly under either quantiser, so floating-point rounding enters only in the scaling of Eq. (1) and in the single final conversion of each output to fp64. (ii) Choose the residue system. The entries of C ′ are bounded, |c′ij | ≤ k 22t =: µ; choose pairwise coprime moduli m1 , . . . , mr (coprimality, not primality, is what CRT requires; the shipping FP8 implementation uses 9–11bit square and near-prime moduli, 1089, 1024, 961, 841, 625, 529, 511, 509, 503, 499, 491, 487 at r=12 [27, 26]) with M =
r Y
mi > 2µ,
(2)
i=1
so that C ′ is uniquely determined by its residues, taken in the half-open symmetric range [−M/2, M/2). The same convention is used for every per-modulus residue in this paper: intervals [−x/2, x/2), with the tie map stated once and for all as “for even x the endpoint −x/2 is included and +x/2 excluded.” (iii) One exact low-precision product per modulus. For each modulus form the residue operands and their exact products, A(i) = A′ mod mi ,
B (i) = B ′ mod mi ,
C (i) = A(i) B (i) mod mi ,
(3)
where every entry of A(i) , B (i) lies in [0, mi ) (or, in the implementations, in the centred range [−mi /2, mi /2), which is what the digit decompositions below require); exactness of the residue products is obtained through the low-precision decomposition of the next paragraph, not by representing full residue products directly, and rounding-free FP32 accumulation holds for k ≤ 216 , with k-blocking (and its recombination cost) beyond that bound [27]. This is where the tensor cores do all the heavy work. (iv) Reconstruct by CRT in Garner form. With precomputed constants wj = (m1 m2 · · · mj−1 )−1 mod mj [4, 6], the mixed-radix digits of each output element and the element itself are v1 = c(1) ,
vj = c(j) − v1 + v2 m1 + · · · + vj−1 m1 · · · mj−2 c′ = v1 + m1 v2 + m2 (v3 + · · · + mr−1 vr ) ,
wj mod mj ,
(4) (5)
evaluated in bounded, fixed-width multi-limb int32 arithmetic (Horner form), lifted to the symmetric range, unscaled per Eq. (1), and rounded once to fp64. Reconstruction touches only the mn outputs and is amortised over the inner dimension k; deconstruction—the residue formation on the left of Eq. (3)—touches every input. That asymmetry is what this paper is about. The fp8 variant. On fp8 tensor cores a residue modulo a 9–11-bit modulus does not fit a single e4m3 significand. The source method therefore represents each residue by two limbs at a piece width b, d = d1 2b + d0 , 0 ≤ d0 < 2b , (6) with b chosen so that every pairwise limb product is exact in the fp8 MMA datapath. Two limbs give three products, and the two modulus kinds reach them differently. For a square modulus m = s2 the d1 e1 term is annihilated modulo s2 , so the two stored planes suffice directly, with 6
epilogue coefficients {1, s, s}; for a nonsquare modulus a third, sum plane d0 +d1 is stored and the Karatsuba identity is evaluated [27]. Either way the count is three plane-product MMAs per modulus, so the stored-plane load pstore (two or three bytes per modulus per scalar) and the issued-MMA count pmma = 3r are different quantities and are kept apart throughout this paper; the companion Part 1 freezes the resulting kernel, including the accumulator ledger these identities imply. In the accurate mode one additional fp8 GEMM estimates an input magnitude bound (it is not a Karatsuba pass), and fast/accurate modes trade the modulus count r at comparable accuracy [27]. The accurate-mode schedule thus costs α = 3r + 1 fp8 MMA passes per fp64 product
(= 37 at r=12),
(7)
the working point that suffices for fp64-equivalent accuracy on well-scaled data; the source analysis reports, under its definition, ≈55 effective significand bits at r=12, ≈59 at 13, and ≈64 at 14, with the precise mode- and set-specific condition given there [27, 13]. Dividing the official dense fp8 rates by 37 gives the emulation ceilings used throughout: ≈ 135 TFLOPS on GB300 (5 PFLOPS fp8) and ≈ 473 TFLOPS on Rubin (17.5 PFLOPS; the primary HGX specification also lists 250-TOPS dense int8, 33-TFLOPS fp64, and the 200-TFLOPS emulated-DGEMM figure itself [20, 18]). The GB300 rates used throughout are NVL72-derived per-GPU values— NVIDIA’s NVL72 rack specification with the sparsity convention unpacked: 720 PFLOPS sparse fp8 and 24 POPS sparse int8 over 72 GPUs give 5 PFLOPS and 166.7 TOPS dense per GPU, and 100 TFLOPS rack fp64 gives ≈1.4 per GPU [19]; the air-cooled HGX B300 SKU carries different per-GPU ratings, so “GB300” in this paper always means the NVL72-derived per-GPU operating point. The substrate transition, and why it matters here. Blackwell-class emulation ran on the int8 tensor substrate: B200 supplies 4,500 TOPS of dense int8, and the shipping cuBLAS emulated DGEMM delivered ≈ 150 TFLOPS on it [14]. The int8 substrate is deconstructionfriendly in two distinct ways. In the error-free-slicing realisation of the Ozaki I lineage [22], deconstruction is nearly free: slicing a scaled integer into 8-bit slices is byte extraction—a scale, a round, and a type-pun, a few SIMT instructions per element in total, with no per-modulus arithmetic at all. Even the CRT realisation is cheaper per modulus on int8: the centred residue of a byte modulus (m ≤ 256) fits an int8 operand directly, so the Karatsuba piece split (and its converts) disappears, leaving cq ≈ 13 against the fp8 route’s 16. Per modulus only, however, and that qualification is not cosmetic: an int8-native set must be all-byte, fourteen byte moduli cannot reach the published set’s 111.8 bits (§5.1), and so ri = 15— candidate A, 117.8 bits. The deconstruction load is then cq ri ≈ 195 against the fp8 route’s 16 · 12 = 192, a wash. What the int8 substrate actually buys is a smaller divisor in the arithmetic roof P/α: αint8 = ri +1 = 16 against 3r+1 = 37, a 2.31× advantage per unit of substrate rate; the deconstruction advantage cancels itself. From the Blackwell-Ultra generation onward, however, silicon area has been redirected to low-precision floating point: the int8 tensor rate falls to a residual ∼ 166–250 TOPS while fp8 scales to 5–17.5 PFLOPS, so the only high-throughput emulation route on Rubin-class parts runs through the CRT/fp8 variant [11]— whose deconstruction is not byte extraction: every streamed element must be reduced modulo each of r moduli and split into three fp8 pieces before a single MMA can issue. The cost of that step is the subject of this paper. (Algorithm attribution of the shipping path requires care: the released cuBLAS GEMM emulation is publicly described as Scheme-I slicing [10], while NVIDIA’s cuEST emulation stack now exposes both schemes—slice-count and modulus-count controls, selected by compute capability [16]. The memo-counted conversion SASS [2] exhibits per-modulus reduction structure, i.e. a CRT path; whether that path is the one behind the published B200 figure is a provenance question the memo’s authors have been asked to confirm. Where this paper reads the B200 figure through CRT constants, the reading is a sensitivity
7
illustration conditional on that attribution; the slicing realisation matters below as the clean experimental control, §7.)
3
The Updated TME Model and the Fourth Term
Part 1’s Tensor–Memory Equilibrium model [11], as corrected by the NVIDIA technical review [2], bounds the execution time of an emulated kernel by four terms: Tserial = max
α Wmma β Q0 cq r φ nin , , Plow Bmem Pint
+ γ nout ,
(8)
the tensor time (α = 3r+1), the memory time (β ≥ 1 the bandwidth multiplier on the native reference bytes Q0 ), the deconstruction time—cq integer-pipe instructions per modulus for each of the φnin converted input elements—and the reconstruction latency γ per output. This displayed quantity is Tserial , a service-tail endpoint: the tensor, memory, and deconstruction maxima (already mutually overlapped) followed by an unoverlapped reconstruction tail γ nout . It is not fully serialised—the genuine no-overlap endpoint is the sum Tfp8 +Tmem +Tdec +γnout — and, being a max, it does not by itself upper-bound measured time. We keep it distinct from the ideal service-overlap time that the tables actually use,
Tsvc = max Tfp8 , Tmem ,
I +N I Ndec N S +N S rec , decPS rec PI
,
(9)
in which deconstruction and reconstruction are charged to their host pipes (I the integer/tensor deconstruction pipe of rate PI , S the SIMT reconstruction pipe of rate PS ) and overlapped with the fp8 and memory times. These two are model endpoints, not rigorous bounds on measured time: with perfect service overlap the time is Tsvc , and imperfect overlap moves Tmeas toward— and, absent the assumed service-independence, potentially past—the tail estimate. Table 11 and the companion engine (_compute_rate) evaluate Tsvc , while Eq. (8) is retained as the service-tail sensitivity Tserial ; a measured schedule is what settles where in the band Tmeas falls. The fourth term is the review’s first-order correction: reconstruction is charged per output, but deconstruction is charged per input, and for memory-bound kernels inputs outnumber outputs by orders of magnitude. On the shipping cuBLAS conversion path the review counts cq = 16 instructions per modulus per element (13 counted from the shipping conversion SASS plus 3 for the fp8 Karatsuba piece converts), against an integer-pipe budget of Pint ≈ 75 TOPS counted on GB300 (instruction throughput for the specific conversion sequence; vendor “TOPS” conventions differ, and the normalisation must accompany any fitted value) and presumed unchanged on Rubin pending the pipe census of the follow-up measurement work [2, 11]. Two consequences of Part 1 frame what follows. First, a streamed kernel that saturates HBM consumes elements at Bmem /8 per second, so the SIMT pipe can spend at most 8Pint /(rBmem ) instructions per modulus per element without falling behind: 6.25 on GB300 but only 2.3 on Rubin—the budget shrinks with the generation, because bandwidth grew 2.75× while the integer pipe did not. Second, dissecting the cq = 16 into its six stages (Figure 1) shows the bulk of the cost—the per-modulus reduction, ∼8–10 of the 13 counted instructions—is linear P over the byte planes of the scaled integer, x mod m = k bk (28k mod m) mod m, and therefore executable as matrix arithmetic on the tensor side; what must remain on SIMT (scale/truncate, the final reduction to the centred residue interval, and the structural Karatsuba split) sums to a proposed SIMT-residual count of cres ≈ q
4 r |{z}
scale/truncate
+
1 |{z}
limb recomb.
+ |{z} 2 + final mod
3 |{z}
≈ 6.3
instr. per modulus at r=12,
piece split
(10) where the limb-recombination term is required by the exact two-limb encoding of §5 (the actual moduli exceed one byte), and each per-instruction count is an instruction-ledger projection to 8
be verified at SASS level, not a measured value. (The lazy-reduction variant of Part 1 lowers this to ≈4.3 for the published set—≈3.3 on the one-pass all-byte set—but doubles the tensor cost to α′ ≈ 6r+1; benign for memory-bound kernels, it would halve the dense ceiling and is therefore excluded throughout this paper.) The reconstruction term, quantified. The fourth term γ nout cannot be waved away with an “O(r2 /k), negligible” aside, and we give it a rate bound. Garner reconstruction is Nγ = r(r−1)/2 + 2r ≈ 90 narrow integer operations per output at r = 12. Hosted on the SIMT pipe, its service time as a fraction of the emulated GEMM’s is Tγ Tmma
=
Nγ /2 · Pemu 45 Pemu = , k Pint k Pint
(11)
which at k = 1000 is 8.1% on GB300 and 28.4% on Rubin— not < 1%—and falls below 1% only at k ≳ 8.1×103 (GB300) / 2.84×104 (Rubin). The review is right that SIMT-hosted Garner is a small-k wall. We therefore adopt the SIMT-Garner cost as the paper’s conservative reconstruction cost ceiling—the most expensive implementable reconstruction, and hence a performance floor within the model—a lower bound on the modeled rate, not a measured one, since every other term in the composite remains a projection—and report every trace composite (Table 11) at it. The ≈r-op int8 hosting is now presented as a codesign target, not a delivered count: an exact reconstruction of the ≈112-bit residue state (the supply lemma needs 111.8 bits) requires either a multi-limb integer recombination (≈rL byte-MACs, L ≈ 14; by op-count rL≈168 exceeds Garner’s ∼90 narrow ops, so the comparison is of host-pipe times, rL/PbyteMAC vs Nγ /Pnarrow , which the engine finds comparable to or below Garner—no decisive gain), or a floating-point fractional-CRT recombination (≈r fp32 MACs) whose fp64-exactness is a Part-3 validation obligation. Until that validation the ≈r thresholds (k ≳ 977 (GB300) / 2270 (Rubin)) are reported as a target/sensitivity, and the conservative thresholds are the SIMT-Garner values above. Both hostings are carried in the per-record engine (params.py: _compute_rate, dispatching the reconstruction host like the storage mode); the trace composites of Table 11 are computed at the conservative SIMT-Garner cost floor, applied to both the emulated envelope and the shipping baseline. Relative to the ≈r codesign target they are unchanged for the computeadvantaged records (k ≫ threshold) but lower for the small-k records: some shipping-relative ratios move between the target and guaranteed modes—the guaranteed values being those tabulated in Table 11 (e.g. UMFPACK GB300 1.00 → 1.19; dense-QR Rubin 2.06 → 1.60)—while the vs-native multipliers for the small-k multifrontal and dense-LU fronts on Rubin fall more (there ≈4 → 1 and ≈9 → 3); the codesign target would restore these, pending Part-3 validation (§6.1). The formal fourth term and the companion’s per-record engine now use the same overlap semantics (Tsvc ): reconstruction demand is added to its host pipe’s service time and overlapped with the fp8 and memory times (_compute_rate), with the serialised Tserial (Eq. (8)) retained only as the worst-case upper bound.
4
The Dense Puzzle: a Size Crossover, and Where 200 TFLOPS Falls
Amortisation, quantified. For dense GEMM the folk intuition is that deconstruction cannot matter: O(n2 ) operand elements against O(n3 ) multiply work. The intuition is asymptotically right and quantitatively misleading. For a m×k by k×n product with each operand element deconstructed once and its residue planes reused thereafter, the conversion time is cq r(mk+kn)/Pint against a tensor time α 2mkn/Pfp8 ; conversion stops binding when n̄ :=
2mn cq r Pfp8 ≥ n∗ = , m+n α Pint 9
(12)
(a) Anatomy of the deconstruction cost per streamed element (memo-counted stages) fp64 element
in HBM
1 scale & truncate
∼ 4 /elem (4/r /mod)
2 byte-plane marshal
no arithmetic; LSU/SMEM b/w
3 reduction Σbk(28k mod m)
∼ 8–10 /mod on shipping SASS migrates: pair of K = 8 limb GEMMs (INT8 TC) or packed dp4a
4 final mod to [−m/2, m/2)
∼ 2 /mod
5 Karatsuba split
∼ 3 /mod (3 stored planes; squares: 2)
6 ESC / max pass
amortised
stored FP8 planes
30 (S) / 44 (A) / 33 (E)
SIMT
data movement
linear → migrates
SIMT
SIMT: structural
linear
workspace
cq ≈ 16 per modulus per element on the shipping path (memo-counted estimate, not a throughput measurement); stages 3–5 are what Ozaki 2.5 re-engineers
(b) The cq ladder, its keep-up budgets, and the knees each rung implies Rubin 2.27
per-modulus SIMT keep-up budget (r = 12)
GB300 6.25
cq = 16 S0 shipping (memo-counted)
knees: Rubin 1211 / GB300 346
16 cq ≈ 10.3 L1 direct: pure-SIMT dp4a, published set
knees: 782 / 223
10.33 cq ≈ 7.5 L1 target (reuse schedule unverified)
knees: 568 / 162 (aspirational)
7.5 cqres ≈ 6.3 two-limb SIMT residual (tensor-migrated)
eff. knees (INT8 capacity): 666 / 287
6.3 cq ≈ 7.3 A at dp4a grain (single-limb; 3-plane split)
knees: 553 / 158
7.27 0
2
4
6
r 0 = 15 budgets: GB300 5.00, Rubin 1.82
8
10
12
14
16
SIMT instructions per modulus per streamed element
Figure 1: Anatomy of the deconstruction cost cq per streamed element (reproduced from Part 1 [11]). Panel (a): the six-stage conversion pipeline, colour-coded by disposition—red stages are elementwisenonlinear and stay on SIMT, green stages are linear and migrate to the tensor pipes, amber is data movement. Panel (b): the five audited cq rungs, each drawn against the per-modulus SIMT keep-up budget 8Pint /(rBmem ) that applies to that rung (the paired ticks; the strip above names the two columns). The budget is not one number, because r is not: the four r=12 rungs are charged against 6.25 instructions per modulus on GB300 and 2.27 on Rubin, while the all-byte A rung runs at r′ =15 and is charged against 5.00 and 1.82, which is why its ticks stand to the left of the others. Every bar ends to the right of both of its own ticks; the narrowest miss is the two-limb residual, whose 6.3 all but coincides with GB300’s 6.25. Ozaki 2.5 is the engineering of this pipeline toward its SIMT-residual count, Eq. (10).
the harmonic mean of the two output dimensions clearing a critical size—one independent of the inner dimension k. At the memo-counted cq =16 (r=12, α=37, Pint ≈75 TOPS): n∗ ≈ 346 on GB300
n∗ ≈ 1,211 on Rubin,
10
the knee moving 3.5× outward in one generation because Pfp8 grew while Pint did not. Below the crossover, delivered throughput under conversion–MMA overlap is Pfp8 n̄ Pint PDGEMM (n̄) = min , 3r + 1 cq r
!
(13)
,
and with conversion fully serialised the ceiling is divided by (1+n∗ /n̄); real kernels land between the brackets. The convert-once premise, and its reach. Eq. (13) assumes each operand element is deconstructed once and its planes reused thereafter (λ=1). That is exact only while the output fits one thread-block cluster—an output edge we call the reach R, equal to 256 for the shipping 4×4 cluster of 642 tiles and given in closed form as a function of cluster shape by Eq. (19) in §5.2. A larger output spans several clusters, and on-chip planes cannot be shared across them, so each operand is re-split λA = ⌈n/R⌉, λB = ⌈m/R⌉ times, and the per-element deconstruction load cq r λ/n̄ stops falling—it pins at cq r/R. Substituting n̄ → R in Eq. (13) therefore gives the plateau in closed form, and the substitution is worth doing explicitly because of what cancels: floor Pdec
PI8TC /2 ISIMT = R · min , res aint8 cq r
!
,
capped at
Pfp8 , 3r + 1
(14)
in the symbols of Eq. (16): aint8 the useful int8 MACs per streamed element, PI8TC /2 the int8 tensor MAC rate, ISIMT the SIMT instruction rate, and cres q the SIMT residual of Table 3. The fp8 peak has cancelled. Below the cap, the floor is a product of exactly two things: the reach R, which is a geometry the hardware fixes, and a per-modulus cost the modulus system fixes. Raising Pfp8 moves it not at all. That is the whole co-design argument in one line, and it is the reason §5.2 spends its effort on reach rather than on peak. The two terms bind on different routes—S and E are int8-capacity bound (0.71 and 0.92 TFLOPS per unit of reach), A is SIMTresidual bound (0.95)—so the ranking of routes is itself a function of which host pipe is scarce. Eq. (14) reproduces the engine exactly at every rung of Table 6 (verify/cluster_map.py). The convert-once crossover is therefore an upper envelope, realised up to n̄≈R; beyond it the achieved dense-square rate plateaus at this deconstruction-λ floor for cluster-aligned sizes (a sawtooth in ⌈E/R⌉, not a hard constant; the ragged worst case is reported below). So the arithmetic roof (473 TFLOPS on Rubin) at n̄ ≳ n∗ is attained only in the convert-once (singlecluster, or materialised) regime, and the honest large-DGEMM ceiling is that floor—which §5.2 quantifies (≈0.50 of that same 473-TFLOPS roof for the codesigned routes), shows to be a codesign coordinate, and lifts back to the roof (Fig. 6). We flag it here so the crossover of Eq. (13) and the floor of §5.2 read as one deconstruction curve, not two claims. Note that Eq. (14) is a Pdec-only statement (§5.1): it charges no reconstruction, and at R=256 it gives 243 TFLOPS for A against the 235 that the reconstruction-aware Psvc reports for E as the winning route at k=4096, the HPL-relevant depth of Table 8 (the ordering inverts beyond k≈6777; Table 2). Why the floor lands at one half, exactly—and why minor hardware recovers all of it. The 0.50 reads like a fitted efficiency, the kind of number a benchmark produces; it is not. Write the floor of Eq. (14) for the winning hybrid route E, whose int8-tensor term binds: P floor = R · PI8TC-MAC /aint8 = 256 × 125 T/136 = 235 TFLOPS. Divide by the 473-TFLOPS roof Pfp8 /(3r+1) and every rate cancels into a ratio of three design integers: P floor PI8TC-MAC 3r + 1 256 × 37 = R × × = = 0.497. roof P P a 140 × 136 |{z} {z } | int8 | fp8-MAC {z } 256
1/140
37/136
11
(15)
The reach R=256 is cluster geometry (Eq. (19)); 140 is Rubin’s provisioning ratio between the fp8 and int8-tensor MAC pipes; 136 is route E’s two-limb plane-formation cost in int8 MACs per streamed element. Nothing is fitted, and that the product lands within 0.3 percentage points of one half is arithmetic coincidence: against route E’s own 437.5-TFLOPS roof the same floor is 0.54, and the published set S floors at 0.38 of its 473-TFLOPS roof. The number moves exactly as the integers move—doubling the int8-tensor provisioning lifts the formula to 471 TFLOPS, where route E’s own roof caps it at 0.92 of the common 473-TFLOPS roof; halving it drops the floor to 0.25; and between cluster-aligned sizes the achieved value follows the ⌈E/R⌉ sawtooth. The same three integers explain why Blackwell Ultra already has full recovery: at GB300’s provisioning ratio (≈60, not 140) the formula gives 156 TFLOPS—above route E’s 125-TFLOPS roof there—so the min of Eq. (14) is taken by the roof and the floor never binds (the engine reports E attaining its GB300 roof exactly; the 135 TFLOPS GB300 floor the tables quote is set S’s α=37 roof, which the dp4a route L1d attains at the memo’s Pint —see the sensitivity in Appendix B). Rubin’s one half is thus a provisioning statement, not a method limit: one generation widened the fp8:int8-tensor ratio from ≈60 to 140 while plane formation stayed charged to the narrow pipe. Figure 2 draws the data movement behind that ledger. The floor exists only because plane formation is charged to the arithmetic pipes: an output larger than one cluster’s reach re-splits its operands (λA =⌈n/R⌉, λB =⌈m/R⌉) because on-chip planes cannot be shared across clusters, and materialising them through memory instead is feed-starved: at the 64 × 64 output tile a materialised plane stream costs p/TILE bytes per useful flop, capping route E at 171 TFLOPS from L2 and 43 from HBM—both below the 235 the on-the-fly schedule holds (§6.2). The minor datapath addition of §5.2 (Option C, an in-flight-convert copy engine) moves formation onto the copy path at stream rate: deconstruction leaves the arithmetic ledger entirely, the λ term vanishes rather than shrinks, and the rate returns to min(roof, memory bound)—the roof, at unchanged reach—while the added deposit traffic is a fraction of an operand read the kernel already performs (§5.2). Full recovery from a minor addition, because the obstruction was an accounting assignment, not a bandwidth shortage. Where the published figure falls. Figure 3 plots Eq. (13). Reading it at cq = 16: the modeled rate is 200 TFLOPS at n̄ = 512 overlapped, and 200 TFLOPS at n̄ ≈ 900 serialised. NVIDIA’s published Emulated-DGEMM figure is reproduced by the deconstruction term alone, with the fp8 tensor pipes idle 58% of the time, for benchmark sizes anywhere in the n̄ ≈ 500–900 window. We state the reading precisely: Hypothesis: the published ∼200-TFLOPS Rubin Emulated-DGEMM figure lies on the deconstruction-limited branch of Eq. (13)—a deconstruction-bound, not tensorbound, operating point at a model-inferred size. This is a hypothesis, not a finding, and the alternative is real: a flat ∼40% delivered efficiency (tile efficiency, sustained-versus-boost clocks, workspace traffic) explains the same single number if the benchmark was large. The B200 figure cuts both ways, and we read it as a sensitivity illustration, not an anchor. Its int8-substrate emulated DGEMM delivered ≈ 150 of a 281TFLOPS ceiling (≈ 53%)—a ratio so similar to Rubin’s that a common, generic library-efficiency margin is the parsimonious first reading. A CRT accounting can nonetheless be laid over it: at cq ≈ 13 (the memo’s Karatsuba-free count) and αint8 = ri +1 = 16 on the all-byte set (ri =15; §2)—with the Pint ≈ 75-TOPS normalisation transferred from the GB300 memo, not measured on B200—the B200 path would have a knee of its own at n∗ ≈ 13 · 15 · 4,500/(16 · 75) ≈ 730, and its published 150 TFLOPS would fall on its serialised branch at n̄ ≈ 840: the two generations’ margins are then jointly consistent with the fourth-term reading at a single, common benchmark size in the upper half of that window, n̄ ≈ 840–900. This illustrative joint reading rests on transferred constants and is conditional on the B200 figure’s algorithm 12
(a) Plane formation is charged to the arithmetic pipes, once per reach window HBM FP64 operands
plane formation on SIMT / INT8-TC pipes: ai8 = 136 MACs/elem (E), charged once per window (λA = ⌈n/R⌉, λB = ⌈m/R⌉)
SMEM/TMEM planes per cluster; not shareable
FP8 tensor cores (one reach window)
materialise planes in L2/HBM? feed-starved: caps route E at 171 (L2) / 43 (HBM) TF
OPTION C: copy engine forms planes IN FLIGHT; deposit = T/R = 25% of the operand read
λ off the ledger: roof, reach unchanged
output: ⌈E/R⌉2 windows, R = 256
(b) The ledger: a ratio of three design integers, not a fitted efficiency 37 1 P floor/P roof = R ⋅ (PI8TC /PFP8 ) ⋅ (3r+1)/ai8 = 256 × 140 × 136 = 0.497
(reach ⋅ provisioning ⋅ formation cost)
GB300's provisioning ratio is 60, not 140: there the formula clears the roof, so the floor never binds. route-E roof 437.5
common roof 473
0.50 of the 473 roof
Rubin TODAY: the λ floor 0.25
Rubin, INT8-TC ×12 (dial down) Rubin, INT8-TC ×2 (dial up)
0.92 (own-roof capped)
Rubin + OPTION C (λ off the ledger)
roof, reach unchanged
GB300 today: formula 156 > roof 125 ⇒ roof-bound: full recovery already
1.00 of its roof 0
100
200
300
400
500
modeled rate at the large-square floor (TFLOPS)
Figure 2: Why the large-square floor is one half on Rubin—and why a minor datapath addition recovers all of it. (a) Data movement at the floor: each reach-sized output window re-splits its operand panels into residue planes on the arithmetic pipes (the λ charges), because planes cannot be shared across clusters and streaming materialised planes from L2/HBM is feed-starved, capping the rate below the floor itself; the dashed Option-C path forms planes on the copy engine instead, removing the charge. (b) The ledger: Rubin’s floor as the derived product of Eq. (15) against the two roofs; the provisioning dial (halved/doubled int8-tensor rate) moving the plateau; GB300, where the formula exceeds the roof and the floor never binds; and Option C returning Rubin to the roof at unchanged reach. All quantities are computed by the artifact engine (render_lambda_floor.py); no value is fitted.
being the CRT path (§2 provenance note); if the shipping path is Scheme-I slicing [10], its deconstruction is byte extraction, the 53% margin is generic library efficiency, and only the Rubin side of the joint reading survives—it validates nothing either way. The two explanations thus cannot be separated by the published numbers alone—but their signatures are orthogonal. The fourth term predicts a delivered rate linear in n̄ up to a knee whose position moves with each path’s cq (≈730 on B200, ≈1,211 on Rubin), and predicts no knee at practical sizes for an error-free-slicing int8 kernel, whose per-element byte extraction puts its crossover below n̄ ≈ 25; efficiency artifacts are flat in n̄ on every path. A DGEMM size sweep on Rubin and B200, with a slicing-based int8 kernel as the control arm, decides the question and fits each path’s cq from its knee position in the same run (§7). With the ceiling stated, the rest of the paper is comprehensive about how far each deconstruction algorithm gets—and honest about where each stops. Table 2 is the menu the switched dispatch chooses from: §5–§5.1 design its entries (the modulus sets and reduction routes that bring cq from the shipping 16 toward ≈5–6 and set both the small-size rate and the floor’s height), and §5.2 collects the limitations that bound them (Table 5) and the preferred co-design target that removes the dominant one. Read together, the two tables are the whole design space on one page. 13
modeled TFLOPS (upper envelope)
(a) Rubin (PFP8 = 17.5 PFLOPS, ISIMT = 75 T/s, PI8TC = 250 TOPS) 500
Ozaki II ceiling PFP8/(3r+1) ≈ 473 TFLOPS n * = 1211
single-cluster reach n ̄ = 256: curves achievable to its left; larger outputs → floor
400
published spec: S0 branch, n ̄ ≈ 512–900
software routes: 1.8–2.2 × (cond.)
300
achieved large-output ceiling: deconstruction-λ floor (0.50 roof)
200
illustrative native FP64 ( ≈ 30 TF assumed) S0: shipping path, cq = 16 Ozaki 2.5 two-limb, INT8-capacity-aware L1 direct: pure-SIMT dp4a, cq ≈ 10.3 A at dp4a grain, cq ≈ 7.3 (single-limb; 3-plane split) L1 target cq ≈ 7.5 (unverified; excluded from envelopes) S0 serialised bracket
published “Emulated DGEMM” ≈ 200 TFLOPS
100
0 32
64
128
256
512
1024
2048
4096
8192
modeled TFLOPS (upper envelope)
(b) GB300 (PFP8 = 5 PFLOPS, ISIMT = 75 T/s, PI8TC = 166 TOPS) 140
achieved large-output ceiling: deconstruction-λ floor (1.00 roof)
Ozaki II ceiling PFP8/(3r+1) ≈ 135 TFLOPS
120 100 single-cluster reach n ̄ = 256: curves achievable to its left; larger outputs → floor
80 60
native FP64 ( ≈ 1.4 TF) S0: shipping path, cq = 16 Ozaki 2.5 two-limb, INT8-capacity-aware L1 direct: pure-SIMT dp4a, cq ≈ 10.3 A at dp4a grain, cq ≈ 7.3 (single-limb; 3-plane split) L1 target cq ≈ 7.5 (unverified; excluded from envelopes) S0 serialised bracket
40 20 output-width-bounded band (n ̄ ≤ 2nb)
0 32
64
128
256
512
1024
2048
4096
8192
outer-dimension harmonic mean n ̄ = 2mn/(m+n)
Figure 3: Modeled emulated-DGEMM upper envelopes versus problem size on Rubin (a) and GB300 (b), Eq. (13), at the memo-counted reference constants and under the Ozaki 2.5 ladder (conservative dp4a counts; generated from params.py). The rising curves are the convert-once upper envelope, achievable only up to the single-cluster reach (green dashed, n̄=256); a larger multi-cluster output is re-split on the fly, so its achieved rate is capped at the solid brick deconstruction-λ floor (≈235 TFLOPS, 0.50 of set S’s 473-TFLOPS roof on Rubin; equal to GB300’s 135-TFLOPS set-S roof, via dp4a: roof-bound), and the shaded gap between floor and roof is unreachable without the co-design of §5.2. We draw the floor here, with the crossover, so the honest large-output ceiling is visible at the outset rather than disclosed late. The shaded band is the overlap bracket of the shipping path (S0, cq =16); the published ∼200-TFLOPS figure lies on the S0 branch at n̄=512 (overlapped) to n̄≈900 (serialised)—a model-inferred size, not a published one. The two-limb route is drawn at its int8-capacity-aware effective knee (≈670 on Rubin, ≈290 on GB300, both at ηred =1); L1 direct (cq ≈10.3) and A at dp4a grain (cq ≈7.3) use no tensor capacity; the aspirational L1 target (cq ≈7.5, unverified reuse schedule) is the thin dotted curve, excluded from all envelopes. The flat dash-dotted line is an illustrative native-fp64 reference (Rubin: ≈ 30 TFLOPS, the 33-TFLOPS vector specification at ∼0.9 blocked efficiency; GB300: ≈1.4 TFLOPS)—the fallback path emulation does not accelerate. The shaded vertical band marks the n̄ range of output-width-bounded kernels (block-vector GEMMs and width-bounded panels, n̄ ≲ 2nb for widths nb = 64–256); the measured traces of §6.1 place LOBPCG and QR-panel work in this band and multifrontal fronts at n̄ ≈ 400–900. In the band the switched dispatch of Figure 5 is worth ≈1.9–2.4× over S0 (platform-dependent; band dispatch—distinct from the schedule-filtered trace composites of Table 11) (independent tensor service assumed; the serialised bracket and the int8-tensor-contention-free dp4a floor of §6 and the storage-mode accounting of §6.2 apply). On GB300 every plotted route exceeds the native reference throughout the evaluated range n̄ ≥ 32.
14
Table 2: The Ozaki 2.5 route menu (Rubin), dispatched by size and platform. “pipe” = where deconstruction runs; “floor” = the 16-CTA large-output deconstruction-λ ceiling (Fig. 3) as a fraction of that route’s own roof; “roof@C” = cluster size at which this route reaches its own roof. Tensor routes carry a high per-element deconstruction on the fast INT8 pipe (higher floors); dp4a routes carry a tiny count on the slower SIMT pipe (lower floors, but zero integer-tensor demand—decisive on GB300 and deployable on any CUDA core today). No route reaches the raw roof for large DGEMM without the co-design of §5.2. The floor entries below are evaluated at (m, n, k)=(4096, 4096, 4096), the HPL-relevant panel depth of Table 8 (reported GPU HPL practice is NB≈892–1024). The floor is a large-k quantity, and this reorders the routes relative to the single-cluster crossover of Table 10: at HPL-scale reduction depth the Garner reconstruction epilogue is not yet amortised, so the all-byte A route—which peaks at 243 TF at n̄=256 (its crossover) and rises to that only as k→∞—delivers a floor of 231 (its 0.61 column, = 231/380), below the hybrid E floor of 235; E is therefore the best floor in this regime, while A keeps the best roof-fraction. This ordering is itself finite-k: route A’s finite-k rate keeps rising past route E’s plateau and overtakes it at k≈6777 (Rubin; the crossover is invariant to m=n at this scale, including the 4096 this table’s floors use), roughly 6–7× beyond any panel width this paper reports as practiced, so E is the better choice at every HPL-relevant depth but not as k → ∞. The n̄=256 crossover value and the large-square floor are thus distinct numbers for a route and must not be conflated, and neither should be read as holding for arbitrarily large k. route set (r)
5
pipe
roof floor roof@ where it wins / what limits it (TF) /roof C
S0 S2L
S (12) S (12)
SIMT INT8-tens.
473 473
0.21 0.38
– 128
Et
E (13)
INT8-tens.
438
0.54
64
At
A (15)
INT8-tens.
380
0.61
64
Dt
D (17) INT8-tens.
337
–
–
L1d
S (12)
SIMT dp4a
473
0.32
256
Adp
A (15) SIMT dp4a
380
0.45
128
shipping baseline; high cq =16 highest roof, full FP64; deconstruction-heavy best floor at HPL-relevant k (235 TF); hybrid moduli best roof-fraction; one-pass reduction; lower roof not carried: lowest roof; needs nblk =3 (288>256 KiB at 2) zero int-tensor, any CUDA core today; SIMT-bound zero int-tensor; GB300 leader; SIMT-bound
The Ozaki 2.5 Method
Definition 1 (Ozaki 2.5). Ozaki 2.5-S is the Ozaki II algorithm on the published modulus set—MMA schedule and reconstruction mathematics unchanged; the accuracy contract inherits as proved for the round-to-nearest variant, the implemented truncation rule being a stated obligation (Appendix A)—with the deconstruction path engineered to its proposed SIMT residual and overlapped behind the tensor stream; the codesigned variants (A/E/D, §5.1) retain the reconstruction framework but carry the per-set obligations of Appendix A. We write “Ozaki 2.5” for the family. The engineering discipline is: (O1) each operand element is converted once where the storage mode permits (λ=1; fused modes admit λ ≥ 1, §6.2) and its residue planes reused thereafter; (O2) the linear per-modulus byte-plane reduction leaves scalar code for dot-product grain—a small constant GEMM on an integer tensor pipe, or packed dp4a on SIMT (route L1; Eq. (17))—a route choice dispatched per platform; (O3) the SIMT residual (scale/truncate, final reduction, Karatsuba split) is implemented at DP4A/narrow-integer grain, at the floor of Eq. (10); (O4) conversion of tile t+1 is pipelined behind the MMA of tile t through the asynchronous copy path; and (O5) the lazy-reduction trade is not taken, preserving the full Pfp8 /(3r+1) ceiling. The name marks its place: mathematically it is Ozaki II (hence not “III”), but as an implementation discipline it is the missing half-step between the shipping realisation and the hardwareassisted deconstruction datapaths whose evaluation is follow-up codesign work. Its net effect 15
Algorithm 1 Ozaki25-Dgemm(A, B, C, r) — Ozaki II with the deconstruction path at its proposed SIMT residual. Per-stage costs in comments; stages tagged [SIMT], [TENSOR], [TMA]. 1: input: A, B, C, r; storage mode ∈ {F, M-L2, M-HBM, P}; conversion multiplicities λA , λB ;
workspace location and lifetime
▷ per §6.2
2: setup: coprime moduli m1 ..mr (the published 9–11-bit set—six squares ≤ 332 plus a near-29
tail—or a codesigned set of §5.1); constant R ∈ Z8×r , Rki = 28k mod mi ; Garner constants; ESC exponent scan of A, B ▷ once per call; amortised 3: allocate plane workspaces by storage mode (stored planes: 30 S / 44 A / 33 E): F/splitk—ephemeral CTA/cluster SMEM only, no HBM object, conversion multiplicities λA , λB charged; M-L2/M-HBM—call-scoped workspace of intended residency (L2 or HBM), written once and read back; P—reuse the existing operand-stationary object, its one-time setup charged separately; hybrid—materialise the smaller operand’s planes, stream the larger ▷ O1: Nconv = λA mk + λB kn; ledger and byte ledger kept separate, §6.2 4: for each k-panel p async, double-buffered do ▷ O4: overlaps MMA of panel p−1 5: [SIMT] S ← trunc(DA Ap ) ▷ scale + truncate (source convention): ∼ 4 instr/element 6: [TMA] marshal byte planes B[0..7] ← bytes(S) ▷ data movement, no arithmetic 7: [TENSOR/int8 or SIMT-dp4a] U ← B × R ▷ O2 route choice: two-limb reduction GEMMs (16r ≈ 192 MACs/elem, exact int32 accum. + recomb.)—or route L1, packed dp4a (Eq. (17)); dispatched per §5.1 8: [SIMT] Vi ← Ui − mi ⌊Ui /mi ⌋ ▷ final mod to centred [−mi /2, mi /2): ∼2 instr/mod 9: [SIMT] WA ← karatsuba3(Vi ) ▷ 3 fp8 pieces per nonsquare modulus (2/square): ∼ 3 instr/modulus 10: same for Bp → WB ▷ SIMT residual incl. limb recomb.: res cq ≈ 4/r + 1 + 2 + 3 ≈ 6.3/modulus 11: for i = 1..r; piece-pairs of the (3r+1)-pass schedule of [27] do
[TENSOR/fp8] Cbi += MMAfp8 (WAi , WBi ) ▷ exact: e4m3-exact Karatsuba pieces, bounded accumulation; the compute roof, Pfp8 /(3r+1) = 473 TFLOPS on Rubin b1..r : ≈ r MACs/output as a 13: [int8-TENSOR + SIMT] mixed-radix reconstruction of C length-r recombination on the int8 tensor pipe (the codesigned hosting) plus a few SIMT carry/normalise ops; single rounding to fp64 ▷ γ per output, tax ∝ 1/k (Eq. (11)); int8-hosted < 1% at k ≳ 977/2270 (GB300/Rubin); SIMT-only Garner (90 ops) < 1% only at k ≳ 8.1×103 /2.84×104 14: [SIMT] ESC check; native-fp64 fallback for out-of-range rows ▷ accuracy contract of [27, 25] per variant status, Appendix A 12:
on the model constants is cq : 16 −→ ≈ 6.3 (two-limb SIMT residual), hence ∗
n (Rubin) : ≈ 1,211 −→ ≈ 480–730, n∗ (GB300) : ≈ 346 −→ ≈ 130–290, the ranges spanning the SIMT-residual and int8-capacity limits (GB300’s residual int8 tensor rate, ∼166 TOPS, caps harder than Rubin’s ∼250). Figure 4 illustrates the dataflow and the pipelining; Algorithm 1 specifies the method. O1: convert-once residue-plane workspace. Deconstruction is charged per converted element, and the conversion count is schedule-dependent: Nconv = λA mk + λB kn, with peroperand multiplicities λA , λB ≥ 1 set by the storage mode of §6.2—the materialised and persistent modes (M/P) achieve λ = 1 but pay plane traffic, while the fused modes (F) generate no 16
(a) Ozaki 2.5 dataflow for one k-panel (Algorithm 1) executed on
FP64 panel Ap / Bp
in HBM
HBM / TMA
1
scale & truncate
∼ 4 instr/elem
SIMT
2
byte planes b0. . b7
marshal only
TMA / SMEM
192 MACs/elem, exact INT32; +1/mod recomb.
INT8 tensor
3 pair of K = 8 limb GEMMs [R (0)], [R (1)], recomb. ×256
the 8–10 instr/mod reduction: migrated OFF the SIMT pipe
∼ 2 + ∼ 3 instr/mod
SIMT (DP4A-grain)
3r FP8 residue planes WA, WB
convert once, reuse ∀ MMA
workspace
3r+1 MMA passes
̂ C1. . r exact
FP8 tensor
Garner → FP64 C
INT32, per output
SIMT
4+5
final mod, Karatsuba split
proposed SIMT residual: cqres ≈ 4/r + 1 + 2 + 3 ≈ 6.3 per modulus (to be verified at SASS level) SIMT (nonlinear stages)
tensor-migrated (linear)
data movement
FP8 MMA (Ozaki II core)
(b) Double-buffered pipeline: conversion of panel p+1 hidden behind the MMA of panel p convert (SIMT+INT8 TC) MMA (FP8 TC)
conv p
conv p+1
MMA p−1 (3r+1 passes)
MMA p (3r+1 passes)
conv p+2
MMA p+1 (3r+1 passes) time
̄ ≈ 1/3 of stream rate at n ̄ ≈ 500 required convert rate = PFP8/(37 n): ⇒ conversion fully hidden for n ̄ ≳ n */3
Figure 4: The Ozaki 2.5 method. (a) Dataflow for one k-panel (Algorithm 1): the fp64 panel is scaled and truncated to integers on SIMT, its byte planes are marshaled, the per-modulus reduction runs as a pair of exact 8×r constant limb GEMMs (recombined with weight 256) on the int8 tensor pipe (the 8–10 instructions per modulus that dominate the memo-counted cq , migrated off the SIMT pipe—or, as route L1, kept on SIMT as packed dp4a dot products, Eq. (17), when the integer-tensor rate is the binding resource, as on GB300), the final mod and the Karatsuba split—the irreducible nonlinearities, the cres q ≈ 6.3 proposed SIMT residual (incl. limb recombination)—remain on SIMT, and the resulting stored fp8 residue planes (30 for the published set) are written once per operand lifetime (call-scoped in mode M; persistent in mode P; never written in modes F/split-k, §6.2) and reused by every MMA pass. (b) Double-buffered pipelining hides the conversion of panel p+1 behind the (3r+1)-pass MMA of panel p; the required conversion rate is Pfp8 /(37n̄), about a third of the full stream rate at n̄ ≈ 500.
plane traffic but admit λ > 1. The conversion-instruction ledger and the HBM byte ledger are therefore separate ledgers, related only under a specific schedule; neither “convert once” nor “no plane traffic” holds unconditionally. Each k-panel of A and B is converted into its stored fp8 operand planes—30 for the published set (two per square modulus, three per nonsquare [27]; 44 for candidate A, 33 for E), a ≈3.8× transient footprint over the fp64 panel for S, bounded by the panel size, not the matrix—and the planes are reused by every MMA that touches the panel. 17
For operand-stationary workloads—Krylov and block-Krylov solvers with a fixed A, repeated application in iterative refinement—the planes persist across calls, amortising the conversion to zero over the iteration count. Retained preprocessing of this kind is already exposed by the public GEMMul8 interface [26]; O1 adopts it as a baseline requirement (panel-granular and cross-call persistent) rather than claiming it as new. The discipline’s storage cost and conversion multiplicity depend on where the planes live—fused in the owning CTA, transiently materialised, or persistent across calls—a coupling with the traffic model made explicit, with its mode boundaries, in §6.2. O2: the byte-plane reduction as tensor work—the exact two-limb encoding. The dominant stage of the memo-derived cq is the per-modulus reduction, and it is linear over P the byte planes: x mod mi = k bk (28k mod mi ) mod mi . A one-pass int8 GEMM against Rki = 28k mod mi is not exact for the actual moduli: with mi up to 1089, the constants reach 826, far outside both signed and unsigned 8-bit range. The exact realisation splits each constant (1) into two unsigned byte limbs, R = R(0) +256 R(1) with Rki ≤ 4, and computes two int8→int32 GEMMs (unsigned bytes times unsigned limbs; all dot products bounded by 8 · 255 · 255 < 219 , hence exact), recombined as U = U (0) + 256 U (1) at one SIMT instruction per modulus per element (the +1 of Eq. (10)). The cost is 2 · 8r ≈ 192 MACs per element; a signed-input correction (the scaled operands are two’s-complement) subtracts the precomputed 264 mod mi — equivalently, adds the stored constant (−264 ) mod mi —in the final reduction, with residues carried in the centred interval [−mi /2, mi /2) (for even mi the endpoint −mi /2 is included and +mi /2 excluded; for m=256 this represents −128); the sign and the boundary cases join the exhaustive-per-modulus residue unit test of §7. A one-pass single-limb variant becomes exact only under a sub-256 modulus system—a modulus/encoding codesign (more moduli, cheaper formation); §5.1 takes this up and produces concrete candidate systems. Against the 37n̄ fp8 operations per element of the main MMA stream this is nominally small, but the reduction GEMM must run on an integer datapath (bytes are not exactly representable in e4m3), and its service demand is first-order in the conversion-bound regime: at the SIMT-residual knee the two-limb route demands ∼380 TOPS (two-operation convention), exceeding the preliminary 250-TOPS Rubin int8 tensor rate. Confined to that pipe, the reduction itself limits the knee to n̄ ≈ 730 rather than ≈480; splitting the reduction between the int8 tensor and the SIMT dp4a route below (the GEMMul8 baseline machinery [26]) is a scheduling question the measurements must settle. Logical-shape efficiency is a further open cost: the exact logical shapes are a pair of K=8, N =r GEMMs (or one K=8, N =2r product with the two limb outputs kept separate, then recombined with weight 256—concatenating the limbs along K would sum them unweighted and is not the identity); such narrow shapes do not map to dense MMA tiles at peak, and the vendor int8 figure is a dense-tile rate, not a proven concurrent capacity beside near-peak fp8 work (the two-endpoint concurrency bracket of §6 quantifies both extremes of that uncertainty). Writing ηred (M, K, N ) ≤ 1 for the issued-versus-useful efficiency of this kernel, the capacity knees below assume ηred = 1; we therefore quote knee and throughput ranges, not single values. Formally, the tensor-migrated routes obey n∗tensor
res
cq r Proof aint8 Proof = max , , ISIMT ηred PI8TC /2
Proof =
Pfp8 , α
(16)
with aint8 the useful reduction MACs per element and PI8TC /2 the MAC rate (the vendor TOPS figure counts two operations per MAC); Eq. (16) instantiates the fourth term of the TME model for these routes and is the form evaluated throughout (params.py). The same dot product on SIMT: the dp4a instruction and route L1. The tensor pipe is not the only dot-product engine on the die. Since Pascal, the SIMT integer pipe exposes
18
(PTX dp4a) dp4a(a, b, c) = c +
3 X
(17)
aj bj ,
j=0
one instruction that reads a and b as four packed bytes each and accumulates their exact dot product into int32: four MACs per issue slot—precisely the grain of the stage-3 sum P 8k mod m ). A worked instance at m=487, four planes: the constants (1, 256, 278, 66) i k bk (2 split into byte limbs (1, 0, 22, 66)+256·(0, 1, 1, 0); for x = 0xDEADBEEF, planes (239, 190, 173, 222), two dp4a return 18,697 and 363, and 18,697 + 256 · 363 = 111,625 ≡ 102 = x mod 487. Eight planes cost two dp4a per limb, and a dp4a returns one scalar dot product per instruction, so narrow high-limb coefficients do not make the second limb free. The honest direct count for the published set is therefore ≈ 4/r + |{z} 2 + |{z} 2 + cL1 q |{z} scale
R(0)
R(1)
1 |{z}
recombine
+ |{z} 2 + |{z} 3 ≈ 10.3, mod
split
and this is the reference L1 count used in every table and figure. The high-limb vectors are narrow (entries ≤ 4) and partially shared—509 and 487 have identical ones—so an explicit common-subexpression schedule could plausibly reach cq ≈ 7.5; we carry that figure only as an unverified optimisation target, excluded from the dispatch envelope until a reuse schedule and SASS count exist. The distinction dissolves for byte-range constants: for candidate A of §5.1 every constant fits one byte, the second limb and its recombination vanish, and cq ≈ 4/15 + 2 + 2 + 3 ≈ 7.3 with no unverified term—the split is charged at three, because three operand planes (d0 , d1 , d0 +d1 ) are stored per modulus, matching α′ = 3r′ +1 (an accounting consistency caught in external review: two digits do not imply two stored operands)—the fastest rigorously-counted pure-SIMT route in this paper. L1 and its A-grain sibling are thus O2’s architectural alternatives, not rivals in mathematics: the same limb identities, executed on the SIMT pipe instead of the tensor pipe, trading int8-tensor capacity for SIMT issue slots; which trade wins is a platform question §5.1 resolves by rate ratio. O3: the proposed SIMT residual. What remains on SIMT is the elementwise-nonlinear residue of Figure 1: diagonal scale and truncate-to-integer (∼ 4 per element, amortised over the r moduli), the final reduction of the GEMM-produced partial residues to the centred interval [−mi /2, mi /2) (∼ 2 per modulus), and the Karatsuba split into three fp8 pieces (∼ 3 per modulus, structural for the fp8 substrate). Together with the limb recombination of O2: Eq. (10), cres q ≈ 6.3—an instruction-ledger projection to be verified at SASS level (constant division compiles to multiply-high sequences, not single instructions), implemented at dp4a/narrow-integer grain (Eq. (17)) where the ISA allows. O3′ : closing the nonsquare S-tail layout. One line of that ledger—“the Karatsuba split into three fp8 pieces”—was until now an assumption for six of the twelve published moduli, and it is the line the whole roof rests on: three stored planes per nonsquare modulus is what makes α = 3r+1 = 37 and hence 473 TF, rather than α = 43 and 407 TF. Under the canonical balanced split of Eq. (6) the six nonsquare S-tail moduli 487–511 have no e4m3-exact layout: their balanced digit sums d0 +d1 reach ±22 and pass through ±17, ±19, ±21, ±23, none of which is representable. Earlier drafts of this paper recorded that as a proved impossibility and carried the six layouts as an open SKIP. That was too strong a reading of its own negative result—it proved only that the canonical split fails—and Lemma 1 closes the gap constructively. Lemma 1 (Carry-corrected three-plane layout). Let E denote the integers exactly representable in e4m3, i.e., those of the form M ·2e with |M | ≤ 15 and magnitude at most 448. Let m be a nonsquare modulus with m ≤ 513 and let d range over the centred residues [−m/2, m/2). 19
Take the canonical balanced base-16 split d0 ≡ d (mod 16), d0 ∈ [−8, 7], d1 = (d − d0 )/16; if d0 +d1 ∈ / E, replace (d0 , d1 ) by (d0 − 16ς, d1 + ς) where ς = sign(d0 ). Then d = d0 + 16 d1 still holds identically, and all three stored planes d0 , d1 , d0 +d1 lie in E, with |d0 | ≤ 14, |d1 | ≤ 16, |d0 +d1 | ≤ 22. Every pairwise plane product is therefore bounded by 222 = 484, so the Karatsuba epilogue d e = (1−16) d0 e0 + 16 (d0 +d1 )(e0 +e1 ) + (162 −16) d1 e1 (18) accumulates exactly in fp32 to depth K ≤ ⌊224 /484⌋ = 34 663—i.e., at Kslab = 4096 with Nacc = 3, with an order of magnitude to spare. The bound m ≤ 513 is sharp for this rule: m = 514 fails. The correction is a single predicated add-pair, taken by 28 of 511 residues at m=511; it costs ≈1 SIMT instruction per modulus and is already inside the cres q ≈ 6.3 ledger, because the split was charged at three from the outset. The construction is redundant rather than canonical—|d0 | may exceed 8—which is precisely why it escapes the canonical envelope; the price is that d0 is no longer the unique base-16 digit of d, which nothing downstream requires. Verification is exhaustive, not sampled: for each of the six moduli we check Eq. (18) over all m2 ordered centred pairs (237 169 to 261 121 per modulus), obtaining zero mismatches and confirming plane representability at every residue (verify/stail_layout.py). Two consequences are worth stating plainly. The published set’s roof of 473 TF no longer rests on an unexamined property of a third party’s stored-plane construction. And the “m ≤ 289 nonsquare” e4m3 envelope that earlier drafts quoted in Table 1 was doubly wrong: it is not tight even for canonical splits (the true canonical bound is m ≤ 321, at base 17), and it does not bound redundant ones at all (the reachable optimum is m ≤ 577, at base 18). The table now carries both corrected numbers. Algorithms 2 and 3 write out the two stages that Algorithm 1 compresses into comments. They are given separately because they are the two stages the cost model actually charges—the first sets the deconstruction floor, the second the fourth term γ nout of Eq. (8)—and because the per-line instruction ledger is the object a SASS-level audit has to reproduce. O4: pipelining. The conversion of panel t+1 runs concurrently with the MMA passes over panel t, double-buffered through the TMA/asynchronous-copy path. The required conversion element rate is Pfp8 /(37n̄)—at n̄ ≈ 500 about one third of the full HBM stream rate—so the SMEM marshaling budget that binds fully streamed kernels [11] is not binding here; overlap is what moves a real kernel from the serialised toward the overlapped bracket of Eq. (13). O5: what Ozaki 2.5 deliberately does not do. The lazy-reduction trade (deferring the final mod into the accumulation at α′ ≈ 6r+1 ≈ 73) buys cq ≈ 4.3 on the published set (≈3.3 on the one-pass all-byte set) but halves the dense ceiling to ≈240 TFLOPS on Rubin; it is the right rung for memory-bound streamed kernels and the wrong one for dense GEMM, and Ozaki 2.5 excludes it.
5.1
Modulus/Encoding Codesign: Sub-256 Systems
The two-limb encoding of O2 accepts the published modulus set and pays for it. The converse design question—posed by the representability problem itself—is whether a different modulus system can make the one-pass reduction exact. We report a script-checked design study (constraints: pairwise coprimality; CRT range ≥ the current set’s 111.8 bits; every residue digit and every Karatsuba sum exactly representable in e4m3 under the centred-residue convention (residues carried in the symmetric range [−m/2, m/2), as the symmetric CRT lift already requires); reduction constants 28k mod mi ≤ 255 for a one-pass GEMM). Two supply bounds, computed by exhaustive dynamic programming over prime-factor masks with no cardinality cap (script and result manifest to be released with the paper; an earlier draft reported greedy constructions as bounds, an error caught in external review), close the space. (i) Pairwise-coprime 20
Algorithm 2 Deconstruct(X, {mi }, DX ) — one operand panel to its stored fp8 planes; lines 5–11 of Algorithm 1 written out. This is the stage the deconstruction floor charges for: the right-hand comments are the per-element instruction ledger that sums to cq (Eq. (10)). Charges are per element per modulus unless marked “/el”. 1: input: fp64 panel X; moduli m1 ..mr ; ESC diagonal scale DX ; byte-plane constants pre(0)
(1)
(1)
split into unsigned limbs, Rki + 256 Rki = 28k mod mi with Rki ≤ 4; the signed-input correction σi = (−264 ) mod mi 2: output: fp8 operand planes—2 per square modulus, 3 per nonsquare; p = 30 stored bytes per scalar for set S (44 for A, 33 for E) 3: for all elements x of X in parallel do 4: [SIMT] z ← trunc(DX x); B[0..7] ← bytes(z) ▷ scale + truncate: ∼4/el, i.e., 4/r per modulus 5: for all moduli mi , i = 1..r do P P (0) (0) (1) (1) 6: [TENSOR/int8 or SIMT-dp4a] Ui ← k B[k] Rki , Ui ← k B[k] Rki ▷ 19 O2. Route S2L: 2·8 int8 MACs, both partials < 2 , hence exact. Route L1: 2+2 dp4a (Eq. (17)) (0) (1) 7: [SIMT] Ui ← Ui + 256 Ui − [ z < 0 ] σi ▷ limb recombination + two’s-complement correction: 1 8: [SIMT] d ← Ui − mi ⌊Ui /mi ⌋, re-centred to [−mi /2, mi /2) ▷ O3, final reduction: ∼2 (multiply-high, not a divide) 9: if mi = s2 is square then 10: [SIMT] (d1 , d0 ) ← balanced base-s split of d; store (d0 , d1 ) ▷ 2 planes; d1 e1 dies mod s2 ; epilogue {1, s, s}: ∼2 11: else 12: [SIMT] d0 ← balanced d mod 16 ∈ [−8, 7]; d1 ← (d − d0 )/16 13: if d0 + d1 ∈ / E then ▷ fails only for d0 +d1 ∈ {±17, ±19, ±21, ±23}: 28 of 511 residues at m=511 14: [SIMT] ς ← sign(d0 ); d0 − =16ς; d1 + =ς ▷ one predicated carry, Lemma 1 15: [SIMT] store (d0 , d1 , d0 +d1 ) ▷ 3 planes; epilogue {−15, 16, 240} (Eq. (18)): ∼3 incl. the carry 16: invariant: every stored plane lies in E and every pairwise product is ≤484, so fp32 accumulation is exact to K = 34 663 ≫ Kslab ▷ Nacc = 3 for the nonsquare tail 17: ledger: ∼4/r + 1 + 2 + 3 ≈ 6.3 SIMT instructions per element per modulus = cres q ; the reduction of line 6 is the term that migrates off SIMT ▷ S row refines this to 5.8, §5.1
21
Algorithm 3 Reconstruct({Cbi }, {mi }, D, E) — mixed-radix CRT recovery of one output tile; line 16 of Algorithm 1 written out. This is the fourth term γ nout of Eq. (8): unlike deconstruction it is charged per output element, so its tax falls as 1/k (Eq. (11)) and it is a small-k wall, not a large-k one. exact per-modulus accumulators Cbi , i = 1..r; moduli m1 ..mr ; ESC scales D, E; precomputed inverses µij = (mi )−1 mod mj for i < j 2: choose host ∈ {int8-tensor, simt} ▷ dispatched per §6; both hostings are carried per record 3: if host = int8-tensor then P 4: [TENSOR/int8] y ← i wi Cbi as a length-r recombination against the precomputed mixed-radix weight vector w ▷ ≈r MACs per output; the codesigned hosting 5: [SIMT] carry-normalise y against the mi ladder ▷ a few ops; γ small 6: else 7: [SIMT] v1 ← Cb1 8: for j = 2..r do P Q 9: vj ← Cbj − i<j vi l<i ml µj−1,j mod mj ▷ Garner: r(r−1)/2 multiply–reduce 1: input:
▷ 2r further ops; Nγ = r(r−1)/2 + 2r ≈ 90 at r=12 l<j ml 11: [SIMT] rescale by D −1 E −1 and round once to fp64 ▷ single rounding: the accuracy 10:
y←
P
j vj
Q
contract of [27] 12: [SIMT] ESC range check; re-issue out-of-range rows on the native-fp64 path ▷ Appendix A 13: cost: int8-hosted, < 1% of service time at k ≳ 977 (GB300) / 2270 (Rubin); SIMT-only Garner, < 1% only at k ≳ 8.1×103 / 2.84×104 ▷ Table 11 is priced at the conservative SIMT-Garner floor
squares under the digit constraint (s ≤ 33) supply at most 94.1 bits (optimum, attained at r=11): squares alone cannot span the fp64 range’s 111.8 bits, so some tail is mandatory and the codesign question is only which tail. (ii) Moduli ≤ 64 supply at most 89.9 bits (optimum at r=18; single-digit moduli ≤ 16 at most 19.5): one-pass-per-modulus schedules of that granularity are unreachable at fp64 accuracy—CRT supply, not representability, is the hard wall. The word exhaustive is scoped in two parts: the two supply bounds are exhaustive—feasibility upper bounds over their precisely defined admissible universes (pairwise-coprime squares with s ≤ 33; moduli ≤ 64)—whereas the chosen S/A/E/D performance designs are not exhaustive: they are selected points in the far larger, unswept space of mixed moduli, schedules, and encodings. Table 3 collects those selected points and their dp4a reference routes on both platforms; it is the rate table the rest of the paper reads from. Accounting note: the S row of Table 3 refines the uniform 6.3/192 accounting used elsewhere by counting per modulus (squares split at ∼2 instructions, 1024 reduces by masking), giving 5.8/176 and knees ≈440/666; the uniform and refined figures bracket the same conclusion. Three plane counts, not one. The single symbol “number of planes” is used in the literature for three quantities that this design separates, because on Rubin they take three different values and each prices a different resource. Let r be the number of moduli and let a modulus be square when m = s2 . pstore stored planes per scalar: what the operand path must carry from L2 or across the DSM crossbar. A square modulus stores two limbs d = d1 s + d0 ; a nonsquare modulus stores three (d0 , d1 , and the Karatsuba sum plane d0 +d1 ). It prices bandwidth. pmma = 3r operand feeds into the tensor pipe. A square modulus has weights {1, s, s}, so its 22
Table 3: Candidate modulus systems and dp4a reference routes on both platforms (reduced model, ηfp8 =1; script-generated conditional projections, produced—with every figure—from the params.py, supplied in the artifact). “mod.@n̄” = modeled upper-envelope TFLOPS of useful fp64 work, evaluated with Pdec-only (n̄)—the deconstruction-limited envelope at k → ∞, carrying no reconstruction charge (§5.1). Every rate column of this table uses that one function; none of its numbers may be read against a Psvc figure quoted elsewhere, and the two rank routes differently at small n̄. S = published set with the two-limb reduction; A = all-byte; E = hybrid (square moduli s2 with s ≤ 33, plus byte tail); D = 7-bit. Italicised rows are the pure-SIMT dp4a realisations of §5. Values assume independent fp8/int8 tensor service; the serialised bracket and the int8-tensor-contention-free dp4a floor are given in the text (θ-bracket, §6). All route values are ideal-overlap envelopes; schedule eligibility (storage mode, shape, capacity) per §6.2; serialisation endpoints in the text. Claim status: modeled projection. system
r′
α′ roof (TF)
reduction
MACs/el
cres q
knee n∗ mod.@256 mod.@512
Rubin (17.5-PFLOPS fp8; int8 cap 250 TOPS; native ≈30 TF illustrative) S: published + two-limb 12 37 473 two-limb 176 5.8 A: all-byte 15 46 380 one-pass 112 5.3 E: hybrid squares+byte 13 40 438 mixed 136 5.3 D: 7-bit 17 52 337 one-pass 128 5.2 S + dp4a (L1 direct) 12 37 473 pure SIMT — 10.3† A + dp4a (single-limb) 15 46 380 pure SIMT — 7.3†
666 401 476 400 782 553
182 243 235 216 155 176
364 380 438 337 310 352
GB300 (5-PFLOPS fp8; int8 cap 166 TOPS; native ≈1.4 TF) S: published + two-limb 12 37 135 two-limb 176 A: all-byte 15 46 109 one-pass 112 E: hybrid squares+byte 13 40 125 mixed 136 D: 7-bit 17 52 96 one-pass 128 S + dp4a (L1 direct) 12 37 135 pure SIMT — A + dp4a (single-limb) 15 46 109 pure SIMT —
287 147 205 148 223 158
121 109 125 96 135 109
135 109 125 96 135 109
5.8 5.3 5.3 5.2 10.3† 7.3†
† Total c , entire conversion on SIMT. L1 direct = 10.3: 4/r + 2+2 dp4a (two limbs) q
+ 1 recombine + 2 final mod + 3 split; an explicit high-limb reuse schedule could lower this toward ≈7.5, an unverified target (§5). A’s byte constants need no second limb; its 7.3 (three-plane split charged, matching α′ = 3r′ +1) carries no unverified term. Knees are model outputs at stated rates; read them as ≈two-significant-figure bands.
high limb is fed twice; the count is therefore three per modulus regardless of how many distinct planes were stored, and pmma > pstore exactly on the squares. With the single reduction pass this gives the familiar α = 3r + 1. It prices issue. Nacc fp32 accumulator tiles that must stay resident across a whole k-slab, hence in tmem. It is 2 per modulus wherever the merging −2b pre-scale is taken and 3 where it is not, and it is not pstore , pmma , or 2r in general. Whether the merge is taken is a per-modulus decision with two gates—the pre-scaled plane must be substrate-exact, and the merged term must leave accumulator margin at Kslab —and on S’s 9-bit tail it is the second gate, not the first, that we decline (below). It prices capacity. Table 4 gives the three counts for the three tensor routes. The S row is the one that resists a mnemonic: its six 9-bit nonsquare tail moduli (487–511) stay at Nacc =3, so S totals 6(2) + 6(3) = 30, not 2r = 24, while E and A do reach 2r, at 26 and 30. Earlier revisions justified that exception by int8 range—the pre-scaled plane −2b d1 leaves [−128, 127]—which is true of int8 and beside the point here, because every Ozaki product pass in the menu runs on the fp8 substrate. On e4m3 the merge is in fact available: under Lemma 1 the pre-scaled plane reaches exactly −256 = −28 ∈ E on all six, and regrouping (18) as d e = (1−16) [ d0 e0 − 16 d1 e1 ] + 16 (d0 +d1 )(e0 +e1 ) reproduces d e over all m2 ordered centred pairs of all six with zero mismatches. What the merge costs is accumulator headroom: the merged term reaches 4095, so exact fp32 accumulation holds only to K ≤ ⌊224 /4095⌋ = 4097, and the schedule is exact at Kslab = 4096 by a single term—against 34 663 for the unmerged three planes of Lemma 1 and 15 420 for the merged byte-modulus accumulator. We decline a margin of one and keep 23
the tail at three (verify/stail_merge.py). The exception is therefore a design choice with a stated margin, not an impossibility; either way the consequence is that S and A make the same tmem ask from opposite directions, and that Nacc cannot be inferred from r alone—it has to be counted per modulus (verify/nacc_table.py). Table 4: The three plane counts and what each prices, for the three tensor routes at T =64 on the fp8 product substrate. “squares” counts moduli of the form s2 ; pstore is 2 per square and 3 per nonsquare; pmma = 3r always; Nacc is counted per modulus—2 wherever the merging pre-scale is taken, 3 on S’s six 9-bit tail moduli, where it is available but declined for accumulator margin (see text). “op. ratio” is the operand demand at the route’s own roof, pstore /α B/FLOP-per-B/FLOP, which is tile- and platformindependent; the corresponding remote (dsm) share is (1−1/c) times it, 0.75 at c=4, and is not tabulated separately. “tmem nblk =2” is the live tile count under modulus blocking (eq. 20), at 16 KiB per 642 fp32 tile, against a 256 KiB budget. Claim status: counted, script-checked. route
r
S (published) 12 E (hybrid) 13 A (all-byte) 15
squares pstore pmma 6 6 1
30 33 44
36 39 45
α
Nacc
37 6(2)+6(3) = 30 40 13(2) = 26 46 15(2) = 30
op. ratio
tmem nblk =2
0.81 0.83 0.96
15 tiles, 240 KiB 14 tiles, 224 KiB 16 tiles, 256 KiB
Read across the rows, the table is also the case against A as a design point, and it is worth being blunt about it because A is the route this paper’s own deconstruction analysis makes look best. A holds exactly one square in fifteen moduli (256 = 162 ), against six in twelve for S and six in thirteen for E, so it pays pstore = 44 against pmma = 45: one modulus in fifteen is amortised on the operand path and the other fourteen are not. Its operand ratio 0.96 leaves four percent of margin at its own roof, where S and E leave nineteen and seventeen. And its blocked tmem footprint is 256 KiB against a 256 KiB budget—exactly zero headroom, so any co-resident use of tmem, any tile larger than 642 , or any less-than-perfect allocator makes A infeasible rather than merely slow. And its roof is the lowest of the three to begin with, 380 TF against S’s 473. Three independent margins—operand bandwidth, tmem capacity, and the roof itself—close on the same route, and not one of them is the deconstruction cost that makes A look attractive in the first place. Two performance functions, kept apart. Route comparisons in this section and the next are reported through two distinct functions, and a good deal of confusion—including in earlier drafts of this paper and its companion—comes from reading a number computed by one as though it were the other. They are: Pdec-only (n̄) the deconstruction-limited envelope: the rate at which residue formation and the MMA schedule can supply a problem of edge n̄ in the limit k → ∞, with no reconstruction charge. This is the right function for asking which modulus set deconstructs most cheaply, and it is what the knees, crossovers and envelope figures of this section report. Psvc (m, n, k) the reconstruction-aware service rate at finite k: the same schedule with Garner reconstruction charged against the mn outputs it actually touches. This is the right function for asking what a library call would deliver, and it is what the floor, the rung table, and every “delivered” figure report. The two agree only as k/n̄ → ∞, and they can rank routes differently: at n̄=256 the all-byte set A leads on Pdec-only and trails on cubic Psvc , for the reason given below. No table row in this paper mixes them, and each is labelled with the function it uses. Candidate A (all-byte). Fifteen coprime moduli {256, 255, 253, 251, 247, 241, 239, 233, 229, 227, 223, 217, 211, 199, 197} give 117.8 bits; the largest reduction constant is 243, so the one-pass int8 GEMM is exact as written (112 MACs per element; 8 · 255 · 243 < 219 ). Byte 24
moduli also make Karatsuba representability clean, provided residues are carried centred: two balanced base-16 digits span |d1 · 16 + d0 | ≤ 136, which covers the centred range |x| ≤ 128 of every byte modulus (it would not cover uncentred residues up to 255—the digits would overflow). Under the canonical tie rule d0 ∈ [−8, 7], d1 = (x − d0 )/16 the decomposition is unique, the digits lie in [−8, 8], the stored sum plane d0 +d1 lies in [−16, 15], and every integer in these ranges is e4m3-exact (verify_candidates.py)—unlike the published nonsquare tail, whose balanced digit sums reach ±22, where 17, 19, 21 are not representable. We should be precise about how much this buys, because Lemma 1 has narrowed it: the tail is no longer unlayoutable, only no longer canonically layoutable. A’s advantage over S here is one predicated add-pair on ≈5% of residues, not a feasibility gap; the honest case for byte moduli rests on the capacity and deconstruction terms below, not on representability. The price is the pass count: α′ = 46 (+24%), roof 473 → 380 TFLOPS (−20%), and three stored planes per modulus (the split charged at three accordingly). The return is the deconstruction side improving substantially, the effective knee falling 666 → 401 (now SIMT-residual-bound; the int8-capacity knee is 341); below n̄ ≈ 540, A is projected faster than the published set despite the lower roof (243 vs 182 TFLOPS at n̄=256, +34%). That last comparison is a statement about Pdec-only , the deconstruction-limited envelope at k → ∞, and it does not survive transfer to the reconstruction-aware service rate—so we make the distinction explicit here rather than let the reader carry the wrong number forward. A’s advantage is real where k is deep: at n̄=256 and k=4096, Psvc gives 231 for A against 182 for S, +27%. In the cubic regime k=n̄ it inverts, and sharply: 131 for A against 168 for S, −22%. The mechanism is the modulus count A buys its exactness with—Garner reconstruction grows like r2 , so r=15 carries a residual that k=n̄=256 cannot amortise. A is therefore a deep-k candidate, not a small-square one, and the two functions are tabulated separately throughout (§5.1); no row of this paper mixes them. Candidate E (hybrid)—the projected sweet spot. Keep the six square moduli (1089 down to 529: squares of s ≤ 33, where the square shortcut needs no third plane and, again centred, two balanced base-s digits with |di | ≤ 16 span 16s+16 ≥ (s2 −1)/2—exactly tight at s=33, which is presumably why 332 tops the published set), and replace only the nonsquare tail with seven byte moduli {251, 247, 241, 239, 233, 229, 227}: r′ = 13, α′ = 40. The roof concedes 7.5% (473 → 438 TF), the knee falls to 476, modeled throughput at n̄=512 rises 364 → 438, and on the deconstruction envelope E dominates S everywhere below n̄ ≈ 620 (above which the two-limb S passes E’s roof). The reconstruction-aware picture is again different, and in E’s case it is the more favourable of the two at scale: under Psvc at k=4096, E leads S at every n̄ we evaluate (235 vs 182 at n̄ ≥ 256), and in the cubic regime E leads S from n̄ ≈ 768 upward while S leads below. E is thus the route we recommend for large problems on either accounting, and the one place S retains an edge—small cubic shapes—is precisely where the regime-switched dispatch below sends the work to S anyway. Two-limb work survives for five moduli only (1024 reduces by bit masking, free). Regime-switched dispatch. The modulus set and the reduction engine are runtime parameters, so a library need not choose: dispatch by n̄, recovering the sub-crossover region at ≤ 7.5% asymptotic cost with the Ozaki-II reconstruction framework unchanged (CRT range matched or exceeded in every candidate; per-set proof obligations in Appendix A). The dispatch menu comprises the tensor-migrated routes (two-limb S; one-pass A; mixed E) and the pure-SIMT dp4a routes of §5 at their conservative counts (L1 direct, cq ≈ 10.3; A at dp4a grain, cq ≈ 7.3 with the three-plane split charged). On Rubin the roomy 250-TOPS int8 tensor rate keeps the tensor-migrated routes ahead everywhere: A below n̄ ≈ 410, E to ≈620, and the two-limb S above, on the roof from ≈670; the dp4a routes are dominated there (L1 direct trails twolimb, and A-dp4a trails one-pass tensor A). On GB300 the residual int8 rate (∼166 TOPS) 25
throttles every tensor-migrated reduction (effective knees ≈290/150/210 for S/A/E), and the small-size band compresses: tensor A leads to its 109-TFLOPS roof (capacity knee ≈150), with its pure-SIMT dp4a realisation within ≈7% of it (knee ≈160; no int8 tensor use at all) as the capacity-free fallback; E carries ≈180–210; L1 direct takes the envelope to the full 135-TFLOPS roof at n̄ ≈ 220; and the two-limb S ties on the roof from ≈290, freeing SIMT slots there. Modulus codesign thus pays on both platforms—the all-byte set leads both small-size bands—and on GB300 the dp4a realisation offers nearly the same throughput with zero int8-tensor demand. Should the L1 reuse schedule reach its ≈7.5 target (§5), L1 would take the entire GB300 sub-roof envelope and the Rubin mid-band above ≈530; nothing below rests on that. Caveats as elsewhere in this paper: greedy set selection is not proven optimal, the accuracy contract must be re-verified per set against the source analysis, the 3-pass schedule is assumed for balanced-digit Karatsuba, and ηred , tile shape, and int8/fp8 concurrency remain measurement questions; the candidate sets therefore join the validation sweep of §7. Reverting to int8 as the substrate is not advantageous on GB300: at αint8 = ri +1 = 16 (the all-byte set’s ri =15) the residual rate caps emulation at ∼10 TFLOPS, so int8 serves the reduction GEMM only. Figure 5 shows the resulting dispatch on both platforms: the switched envelope is at least as fast as any single route everywhere, and the arithmetic roof is conceded nowhere.
5.2
Hardware Co-Design: Lifting the Deconstruction Floor to the Roof
The preceding two sections leave a single, sharp question, and it is a co-design question. One term—deconstruction—governs the whole story, and it does so through one quantity: the peruseful-FLOP conversion load ≈ cq r λ/n̄. Inside a single cluster λ = 1, so as n̄ grows this load falls—the crossover of §4, along which the arithmetic roof is approached. But the load cannot keep falling. Growing n̄ past the cluster’s 256-element reach forces λ = n̄/256, at which point cq r λ/n̄ = cq r/256: constant. The crossover curve and the “0.50-of-roof floor” (235 of Rubin’s 473 TFLOPS) are therefore not two claims but one curve—the switched envelope rises to ≈235 TFLOPS at n̄=256 and then plateaus (Fig. 6a). Its λ=1 extrapolation would reach the roof near n̄≈700 (and the shipping cq =16 path’s near n̄≈1200, §4), but that extrapolation is unreachable: any n̄>256 is multi-cluster, so λ>1 and the curve plateaus instead. The practical effect is that a real large dense-DGEMM benchmark on Rubin sees ≈235 TFLOPS, half the 473-TFLOP arithmetic roof—while the traced applications, governed by their tall/skinny, smallk, or sub-cluster shapes, are largely spared (§6.1), and the whole question vanishes on GB300, whose 135-TFLOP roof already sits at or below its own floor. Before acting on that floor, Table 5 lays out every limitation identified for Ozaki 2.5—what each binds, its magnitude, and its fix—so the design space is explicit, and so the reader can see that the deconstruction-λ floor is the one limitation that binds large dense DGEMM, the rest being small-size, platform-, or design-space-specific. Read with the route menu (Table 2), it is the whole picture on two pages. The floor is R Pint /(cq r)—Eq. (13) evaluated at n̄ → R, which is Eq. (14)—so what pins is a per-element deconstruction load Ldec = cq r/R, and the rate is Pint /Ldec under the cap and the host-pipe min(·) of Eq. (14). We write it with r explicit because the load is per modulus as well as per element: the shorthand “cq /reach” used in earlier drafts drops a factor r and is dimensionally incomplete, though it never entered the engine, whose Eq. (13) second term has always been n̄Pint /(cq r). Here the cluster’s reach R is the harmonic mean of the two output edges a cm ×cn cluster of T ×T accumulator tiles spans (Eq. (19) below; T the accumulator-tile side, cm cn = C the CTAs per cluster), so the floor is a co-design coordinate, not a wall, along three hardware axes and one algorithmic (Fig. 6b; the concrete targets are collected in Table 7). These are potential co-design points—none exists on shipping or preliminary-Rubin silicon—and the value of the model is that it localises them precisely: the closed form R Pint /(cq r) names exactly the four quantities a change can move—R and Pint hardware, cq the deconstruction datapath, r the modulus set—and no other. 26
modeled TFLOPS (upper envelope)
(a) Rubin: dispatch A → E → S at n ̄ ≈ 414 / 616 500
dispatch A (tensor)
roof 473
E
dispatch S (two-limb: on roof, frees SIMT)
roof 438
400
roof 380
300
dp4a routes dominated on Rubin (omitted); shipping path falls below native for n ̄ ≲ 80
200
switched envelope (dispatch by n)̄ S: published set + two-limb (tensor) A: all-byte, one-pass (tensor) E: hybrid squares + byte tail (tensor) original shipping path, cq = 16 illustrative native FP64 ( ≈ 30 TF assumed)
100
0 32
64
128
256
512
1024
2048
4096
8192
modeled TFLOPS (upper envelope)
(b) GB300: dispatch A → E → L1 → S at n ̄ ≈ 178 / 207 / 287
140
roof 135
dispatch A (tensor)
E
L1
dispatch S (two-limb: on roof, frees SIMT)
roof 125
120
roof 109
166-TOPS residual INT8 cap throttles the tensor routes; A (tensor) leads the small-size band, A (dp4a) within ≈ 7%
100 80
switched envelope (dispatch by n)̄ S: published set + two-limb (tensor) A: all-byte, one-pass (tensor) E: hybrid squares + byte tail (tensor) L1: published set, pure-SIMT dp4a (cq ≈ 10.3) A at dp4a grain, pure SIMT (cq ≈ 7.3) original shipping path, cq = 16 native FP64 ( ≈ 1.4 TF)
60 40 20 0 32
64
128
256
512
1024
2048
4096
8192
outer-dimension harmonic mean n ̄ = 2mn/(m+n)
Figure 5: The regime-switched dispatch on both platforms (reduced model, int8-capacity-aware effective rates, conservative dp4a counts; modeled convert-once upper envelopes, not delivered performance— realised only up to the single-cluster reach (n̄≲256); a large multi-cluster DGEMM is capped at the deconstruction-λ floor (0.50 of the 473-TFLOPS roof on Rubin), Fig. 3 and §5.2). (a) Rubin: below n̄ ≈ 410 the all-byte set A is fastest (one-pass exact reduction, steepest effective slope); the hybrid E carries the middle band to ≈620; the two-limb S takes the dispatch above, on the roof from ≈670. The pure-SIMT dp4a routes are dominated on Rubin and omitted for clarity. (b) GB300: the ∼166-TOPS residual int8 cap throttles the tensor-migrated reductions and compresses the small-size band: tensor A leads to its 109-TFLOPS roof (its dp4a realisation, shown dotted, within ≈7% with no int8-tensor use), E carries ≈180–210, L1 direct reaches the full 135-TFLOPS roof at n̄ ≈ 220, and the two-limb S ties from ≈290 (freeing SIMT). The shaded envelope is the dispatched performance: no single-route curve exceeds it anywhere. For comparison, the original shipping path (cq =16, dotted) and the native-fp64 reference lines (illustrative, flat) are shown: on Rubin the shipping path falls below native for n̄ ≲ 80 while the switched envelope clears it from n̄ ≈ 32 (≈35 against the official 33-TFLOPS figure [20]); on GB300 every plotted route exceeds the native reference throughout the evaluated range n̄ ≥ 32.
Reach is the largest-headroom axis, with two hardware routes to the same on-the-fly reach— raising C or raising T . Reach is an integer-geometry quantity, and we state it as one rather √ than through the continuous idealisation T C used in earlier drafts: a cm ×cn cluster of square T ×T tiles reaches 2 T cm cn R = , (19) cm + cn 27
Table 5: Every limitation identified for Ozaki 2.5: what it binds, its magnitude, and what lifts it (software / hardware). The deconstruction-λ floor (row 2) is the one that binds large dense DGEMM; the rest are small-size, platform-, or design-space-specific. limitation
binds
magnitude
lifted by (software / hardware)
deconstruction cq (4th TME term) deconstructionλ floor
everything
16→5–6→0
large dense DGEMM (HPL) small-k kernels
0.50 of 473 (Et); 0.21–0.61 of own roof
modulus codesign (SW); in-flight convert (HW) reach↑; L2+blocking; cq →0
reconstruction (Garner) tax TMEM capacity (256 KiB)
plane-feed / materialise cluster reach (256) SIMT–fp8 overlap ρ INT8 residual cap moduli supply (111.8 bit) small-n̄ knee (crossover)
schedule choice
SIMT 111%@256; INT8 2.2%@1024 at nblk =2: 240/224/256 S/E/A (A exact); D 288 needs nblk =3 needs 222 TB/s; L2 88
the λ itself dp4a routes
R=2T cm cn /(cm +cn ) serialised ≈0.5×
GB300 tensor routes modulus-set design small outputs
compresses small-size band squares 94, sub-64 90 (short) n∗ ≈480–1211
route set + reach
INT8 mixed-radix hosting (SW); large NB larger TMEM (HW)
on-the-fly (SW); L2+blocking (HW) larger cluster / TMEM (HW) measure (Part 3); scheduling (SW) all-byte dp4a (SW) hybrid/all-byte-tail dispatch (SW) lower cq
√ the harmonic mean of its two output edges, which coincides with T C exactly when cm =cn and falls short otherwise (a 2:1 grid loses 5.7%). “Realizable” here carries two constraints, not one. Besides cm cn ≤ C, the grid must be placeable: CUDA requires each grid dimension to be divisible by the corresponding cluster dimension [15], so the admissible (cm , cn ) are those dividing the CTA-grid extents ⌈m/T ⌉ × ⌈n/T ⌉ of the launch. A 40962 output at T =64 is a 64×64 CTA grid, and 64 is divisible by 4, 8 and 16 but by none of 5, 6 or 11. The distinction is load-bearing rather than pedantic, because the product constraint alone admits better reaches than the ones we tabulate: 5×6 gives R = 349.09 against the 8×4 rung’s 341.33 at C=32, and 11×11 gives R = 704 against 16×8’s 682.67 at C=128. Neither divides a 64×64 grid. A shape whose CTA-grid extents happened to be divisible by 5, 6 or 11 could use them; power-of-two cluster edges are the ones that stay admissible across the whole shape sweep, which is why those are what we tabulate, and on the shapes this paper reports the rungs of Table 6 are therefore the admissible optima rather than convenient round numbers. Every entry there is the engine’s own floor at that reach, not an interpolation. The two routes are not physically interchangeable: more CTAs per cluster (C=16 today—Blackwell allows 8 portable, 16 optin—to 64–128) requires a cross-GPC distributed-shared-memory fabric the current within-GPC crossbar does not provide, while a larger accumulator tile T is bounded by the fixed 256-KiB per-CTA TMEM SRAM. What is being asked for on the C axis is fabric reach, not a new instruction. At today’s C=16 the sharing is expressible in shipping PTX: the producing CTA publishes its converted planes into its own shared::cta window and each consumer pulls with cp.async.bulk.shared:: cluster.shared::cta.mbarrier::complete_tx::bytes, addressing the peer window through mapa and completing on a .cluster-scope mbarrier; all three primitives were introduced in PTX ISA 8.0 and require target sm_90 or later [21]. Because each consumer’s MMA then issues against its own shared memory, the plain .cta_group=1 descriptor suffices and no multi-CTA operand descriptor is invoked. The C=64/128 rungs of Table 6 therefore ask the fabric to carry the same instruction sequence further, across GPC boundaries—a physical-reach change, which 28
is why we book it as hardware co-design rather than as software the reader could write today. Part 1 [11] gives the per-CTA byte accounting for this exchange. That TMEM bound deserves a precise statement, because there is a tempting argument that it is not a bound at all. Under the modulus blocking that Part 1 freezes, the r moduli are cut into nblk blocks and only one block is live at a time, so the resident tile count falls from P Nacc = i w(i) to live Nacc (nblk ) = min max B
B∈B
X
w(i) ≥
(20)
Nacc /nblk ,
i∈B
the minimum over partitions B of the modulus set into nblk blocks. The inequality is generally strict: blocks partition moduli, not accumulators, so a modulus contributes its whole width w(i) ∈ {2, 3} to whichever block holds it. At nblk =2 the true counts are 15/14/16 tiles for S/E/A, not the 15/13/15 the ceiling suggests, because r is odd for E and A (verify/blocking.py); an earlier version of this argument quoted the bound as though it were the count. One might then raise nblk until any tile fits. That does not work, and the reason is worth recording: each extra block re-traverses the k-slab and re-reads the fp64 operand at 8/T B per useful FLOP, so the l2 bill is (nblk −1) (8/T )×rate—and the T in the denominator is cancelled by the rate the larger tile unlocks, while the nblk needed to fit a fixed TMEM grows like T 2 . At T =96 the fit needs nblk =5 and l2 affords 4; at T =128 it needs 7–9 and affords 4; at T =192 no nblk ≤ r fits at all, and beyond r blocking stops being operand-free. So the area ask is real, and at the nblk =2 operating point it is larger than earlier drafts claimed: a 1282 tile is a 3.5–4.0× TMEM ask (not ≈3×) and a 1922 tile 7.9–9.0× (not ≈6.75×), the range running over S/E/A. The tile route does still cut the materialised plane feed to p/T (222 → 111 TB/s at 1282 ), so of the two it remains the more leveraged axis—but it is bought with SRAM area, not with scheduling. Table 6: Reach rungs on both escape axes, from the realizable integer geometries of Eq. (19). Rates are the engine’s Rubin reconstruction-aware floor Psvc at m=n=k=4096 evaluated at that reach, against the 473-TFLOPS roof; they are modeled projections, not measurements. The cluster axis holds T =64; the tile axis holds the shipping C=16=4×4. TMEM is the nblk =2 live-accumulator footprint over S/E/A against the 256-KiB per-CTA budget. Reproduced by verify/cluster_map.py and verify/blocking.py. axis
configuration
reach R
TFLOPS
of roof
TMEM ask
cluster C
4×4 (shipping, opt-in) 8×4 8×8 16×8
256.00 341.33 512.00 682.67
235 314 438 473
0.50 0.66 0.92 1.00
0.9–1.0× 0.9–1.0× 0.9–1.0× 0.9–1.0×
tile T
642 (shipping) 962 1282 1922
256.00 384.00 512.00 768.00
235 342 438 473
0.50 0.72 0.92 1.00
0.9–1.0× 2.0–2.2× 3.5–4.0× 7.9–9.0×
L2 service bandwidth is the second axis and the weakest: fed from L2 at the fp8 rate a materialised (λ=1) schedule reaches the roof, but that needs ≈10×HBM, and each set attains its own roof at a different multiple, because the requirement is BL2 = roof · pstore /TILE and pstore differs: set S reaches its 473 TFLOPS at 10.08× (221.7 TB/s), set E its 438 at 10.25× (225.6 TB/s), and set A its 380 at 11.89× (261.5 TB/s). That is against the ≈4× assumed, and it is subject to the L2 capacity predicate p k(m+n) ≤ CL2 (§6.2): the large multi-cluster outputs the floor governs need ∼1 GB of resident planes at k=4096, 8× over the 126-MB L2, so the planes spill to HBM and the ceiling collapses to ≈47 TF. The L2 escape thus reaches the roof only in the moderate-k window (N ≲ 1449 square, set S) where the floor scarcely bites—the least useful of the three hardware asks. cq , the algorithmic-and-silicon axis, is in fact the highest-leverage single lever: it scales the whole curve at every reach, and an Option-C-class datapath (§6.5.2 of [11]) that drives cq → 0 removes the deconstruction term outright—lifting the dense floor to 29
Table 7: Potential hardware co-design targets for dense emulated DGEMM (Rubin); none exists on shipping silicon. The floor is R Pint /(cq r) (Eq. (14)), so each axis is a concrete, measurable ask, and each row here moves exactly one of its four quantities with the other three held at the memo-derived reference constants of Table 3: the two reach rows move R; the L2 row moves the supply the deconstruction load is served from, by making a materialised (λ=1) schedule feasible; the Part 1 datapath row moves cq ; and r is not a hardware knob at all but a modulus-set choice (§5.1). Reach has two routes (cluster size or TMEM tile) to the same on-the-fly value but they hit different physical limits (cross-GPC DSM fabric vs. fixed TMEM SRAM); L2 is gated by an additional capacity predicate; cq —developed in Part 1 as the sparse-critical lever—is the highest-leverage as it removes the term at every reach. Reach targets are integer-geometry quantities (see text): the tabulated C=64 (8×8) and both tile rungs are exact square grids, while C=128 must be 16×8 and so reaches 683, not 724—still above the reach 666 at which the floor meets the roof, so the roof entry holds. axis
knob
reach
CTAs/cluster C 16
reach
TMEM tile T
bandwidth L2 service algorithm
cq
today
target
dense-DGEMM effect
floor 0.50 → 0.92 / roof; needs cross-GPC DSM 642 1282 (3.5–4×) / 1922 (7.9–9×) floor 0.50 → 0.92 / roof; also feed 222 → 111 TB/s 4×HBM 10–12×HBM roof iff planes L2-resident (N ≲ 1449); else → 47 TF memo-count options A/B/C highest leverage: removes term at every reach; sparse lever 64 / 128
the roof and un-binding the conversion-bound sparse kernels, which meet cq directly with no re-split floor to escape. Because it is the sparse-critical lever, the companion Part 1 develops it in full (the three options A/B/C: a rebalanced narrow-integer issue, a cvt.fp64.residues instruction, and a residue-decompose copy-engine mode); this paper cross-references rather than repeats it. Table 7 collects the four axes as concrete, measurable asks, and Fig. 6b plots the two reach routes and the L2 escape against the floor they have to clear. Every carried dispatch route is retained (route D is not carried), because the knob that rescues each differs (Table 2). The all-byte and hybrid sets, with their lower cq , already stand at 0.54–0.61 of their own (lower) roofs today—438 TFLOPS for E, 380 for A— and reach them at a 64-CTA cluster; the high-roof published two-limb and pure-SIMT routes sit lower (0.32– 0.38 of their own, higher, 473) and reach that roof only at a 128-CTA cluster (via the switched envelope; a pure-SIMT dp4a route on its own needs the next cluster step)—so the “best” route is reach-dependent, and a design that can switch modulus sets by shape (§5.1) can also switch by the available reach. The sparse and streaming kernels face the same cq term from the other side—below the threshold intensity they are conversion-bound rather than re-split-bound—and their lever is deferred/persistent deconstruction (convert the stationary operand once, reuse the planes across the iteration), analysed in Part 1 [11]. On shipping and preliminary-Rubin silicon none of the hardware knobs exists yet, so the projections of §6 compare software routes at the 16-CTA reach; the point here is that the 2× dense-DGEMM gap is named and closable, and which improvement closes it is a Part-3 measurable (the per-route floors and the cluster each needs are Table 2). A worked example: HPL on a Rubin cluster. To make the ceiling and its escape concrete, take High-Performance Linpack—the fp64 benchmark whose whole reputation is a high roof fraction. HPL is dominated (≳90% of its FLOPs) by the rank-NB trailing update C ← C −A[M, NB] B[NB, N ], whose output M ×N is the local trailing submatrix under the twodimensional block-cyclic distribution: it fills the GPU and is therefore large and multi-cluster— exactly the regime the deconstruction-λ floor governs. So HPL does not see the 473-TFLOP roof. It is worth being precise about why the panel width NB matters, because it does two 30
(b) the floor is a co-design coordinate, not a wall
(a) crossover and floor are ONE deconstruction curve 600
one curve rises, then plateaus
400
co-design headroom (2 × )
300
200 FP8 arithmetic roof 473 TF λ = 1 ideal: convert once (unreachable past 256) achieved on-the-fly (this design) deconstruction-λ floor 235 TF
100
achievable large-DGEMM rate (TF)
emulated FP64 rate (TF)
500
0
128-CTA 16x8 or 1922 tile 473 TF (1.00)
500 L2 10x → 469
400
300
64-CTA 8x8 438 TF (0.92)
32-CTA 8x4 314 TF (0.66) L2 6x → 282 today: 16-CTA 4x4 235 TF (0.50)
200
100
roof 473 TF on-the-fly floor vs cluster reach materialise ceiling at higher BL2
0 64
128
256
512
1k
2k
4k
200
square output edge = n,̄ k = 4096
300
400
500
600
700
cluster reach R = 2Tcmcn/(cm+cn), integer cm×cn grids
Figure 6: (a) The crossover and the floor are one deconstruction curve: the λ=1 ideal (convert-once, the arithmetic-roof crossover) is unreachable past the 256-element cluster reach, and the achieved on-the-fly rate plateaus at the deconstruction-λ floor (≈235 TFLOPS, 0.50 of the 473-TFLOPS roof; the gap is co-design headroom). (b) The floor is a co-design coordinate: raising the cluster reach (bigger cluster, or bigger TMEM tile) lifts the envelope floor to the roof (64-CTA → 0.92, 128-CTA or 1922 tile → roof); a ≈10×-HBM L2 reaches it through materialisation instead.
unrelated jobs and an earlier version of this passage credited it with the wrong one. As the depth k of the trailing update, NB does not move the deconstruction floor at all: the deconstruction load is cq r(λA mk + λB kn) against 2mnk useful FLOPs, and k cancels. What it does is amortise reconstruction, which is charged per output element and therefore falls as 1/NB. Separately, and in its second role as the block size of the block-cyclic layout, NB sets the granularity of the local edges M, N , so an NB that is a multiple of the reach R=256 keeps those edges cluster-aligned. That second effect is real but small at HPL’s scale: at a local edge of 8192 the misalignment penalty for a ragged M =N =8000 is ≈2%, not the sharp ragged worst case that afflicts small matrices. The first effect is the large one, and Table 8 quantifies it. Table 8: Sensitivity of the modeled per-GPU HPL trailing-update rate to the panel width NB, at a local output edge M =N =8192 on Rubin. Values are Psvc (reconstruction-aware, §5.1) in TFLOPS of useful fp64; bold marks the better route at that NB. Reported GPU HPL practice is NB ≈ 892–1024 [3, 5], i.e. near the NB=1024 column of this table. Claim status: modeled projection from params.py, not measurement. route S (published set, RN contract) E (hybrid, codesigned)
NB=128
256
512
1024
2048
4096
186 169
186 159
182 202
182 234
182 235
182 235
Two things in that table deserve to be said plainly rather than left for a reader to notice. First, the 0.50-of-roof headline (235 of Rubin’s 473 TFLOPS) is an NB-conditional claim: E reaches 235 only for NB ≳ 1024, and at NB=256 it delivers 159, below the published set. The claim survives because reported GPU HPL practice sits at NB ≈ 892–1024 [3, 5], but it survives on an empirical convention, not on a theorem. Second, the two routes carry different accuracy contracts, and the faster one carries the weaker: S is covered by the round-to-nearest Ozaki II theorem, while E’s guarantees are the per-set obligations of Appendix A. The box 31
below therefore quotes both, and a reader who wants the theorem should read the S row— 182 TFLOPS, 0.38 of the 473-TFLOPS roof—as the contract-backed number. The box is an explicit order-of-magnitude reference, not a submission: Top500 does not currently accept Ozaki-emulated DGEMM for fp64 Linpack, so the numbers below indicate only what the per-GPU floor implies at scale. HPL per Rubin GPU (emulated fp64, today’s 16-CTA silicon, NB=1024, Psvc at a local edge of 8192), quoted for both accuracy contracts: route per GPU of 473-TF roof accuracy contract S (published set) ≈182 TF 0.38 round-to-nearest Ozaki II theorem E (hybrid, codesigned) ≈235 TF 0.50 per-set obligations, App. A Both fractions are taken against the common 473-TFLOPS arithmetic roof Pfp8 /37, so that the two routes are directly comparable on one scale. Table 2 instead reports each route against its own roof; there the E floor reads 0.54 of set E’s 438 TFLOPS, which is the same 235 TF. That is ≈6× and ≈8× the ≈30-TFLOP native fp64, and both bracket NVIDIA’s announced ∼200-TFLOP “Emulated DGEMM.” The S contract is componentwise fp64 with exact residuedomain accumulation, for which HPL’s residual check is routinely passed on GPUs. 10,000-GPU cluster (ballpark): DGEMM peak 1.82 EFLOPS (S) or 2.35 EFLOPS (E); at the HPL efficiency empirically established across a decade of large GPU clusters (75–85% over panel factorisation, pivoting, look-ahead, and communication) Rmax ≈ 1.4–1.5 EFLOPS (S) or 1.8–2.0 EFLOPS (E) fp64 sustained—reference figures, no bespoke system model claimed. With the preferred co-design target below (conditional on the Option C requirement table, Table 9, being met): each route reaches its own roof—not a common one. Because Option C removes the deconstruction term cq outright, the ranking is then set by α alone, and the theorembacked published set S (α=37) becomes the fastest at ≈473 TFLOPS, ahead of the hybrid E (α=40, ≈438) and the all-byte A (α=46, ≈380). On the S row the cluster DGEMM peak becomes 4.73 EFLOPS and Rmax ≈ 3.6–4.0 EFLOPS: HPL roughly doubles against today’s E row on the same 10,000 GPUs and reaches ≈2.6× today’s S row—while gaining the round-to-nearest contract rather than trading it away.
The preferred co-design target to reach the roof. Constrain the problem the way a vendor would—minimise silicon and new datapaths, and let the proposed software carry the complexity—and one option dominates under this model; area, power, routing, and concurrency remain a Part-3 obligation (Table 9), so we prioritise rather than establish minimality: Option C, a residue-decompose transform mode in the TMA/asynchronous copy engine (§6.5.2 of Part 1 [11]; illustrated in Fig. 7). It converts each fp64 tile to residue planes in flight, HBM→SMEM, so deconstruction leaves the compute budget entirely—cq →0 on the SIMT/integer pipes—and the deconstruction load cq r/R that sets the floor vanishes at every reach: the emulated rate reaches the roof at today’s 16-CTA cluster, with no larger cluster, TMEM tile, or L2 bandwidth, because re-splitting (λ) is now free (the copy engine re-converts at stream rate). This is conditional on the planedeposit path: Option C zeroes the compute-pipe cq but not conversion time, which becomes Tcvt = max(Tinput , Nscalar NMAC /Ptransform , Tdeposit , Ttensor-read ); “reaches the roof” holds once the transform sustains the stream rate and the deposit sustains its (pstore /8)× share. Those are provisioning requirements on a vendor, not properties we derive, so we state them as a bill of materials: Table 9. At HPL’s panel widths the model puts the result at the roof to within the ≈1–4% Garner reconstruction epilogue, which overlaps on the residual int8 pipe (§3). What Option C actually costs, stated as a requirement. Table 9 is the ask. Three of its rows deserve comment because two of them are smaller than a reader would guess and one is larger. The deposit is small, and exactly so: the SMEM plane-write rate is pstore /R bytes per useful fp64 FLOP while the operand read the MMA already performs is pstore /T , so the new 32
(a) Software path today: deconstruction runs on the compute pipe, paid λ times registers / SIMT lanes
HBM FP64 tile
SIMT/INT convert cq ⋅ λ ops/elem
SMEM residue planes
tensor cores FP8 MMAs
SMEM round-trip
the deconstruction-λ FLOOR: 235 TF, which is 0.50 of the 473 TF envelope roof (Rubin); conversion competes on the SIMT/INT pipe and is repeated once per cluster for a large output Option C removes the SIMT/INT convert and the SMEM round-trip; it does NOT remove conversion time
(b) OPTION C: residue-decompose transform in the TMA / async copy engine (in flight) HBM FP64 tile
convert in flight
COPY ENGINE residue-decompose transform 166–325 TMAC/s aggregate
SMEM planes in MMA layout
tensor cores FP8 MMAs
cq → 0 on the COMPUTE PIPE only. Conversion time does not vanish: Tcvt = max(Tin, NscalNMAC/Pxform, Tdep, Tread). REQUIREMENT on the vendor, each route at ITS OWN roof (S 473, E 438, A 380 TF, not a common 473): transform 166–325 TMAC/s, i.e. 1.3–2.6 times the platform's ENTIRE residual INT8 tensor throughput; input 12–15 TB/s (about 17% of L2 service, not an HBM ask); plane deposit 55–65 TB/s of SMEM writes = T/R = 64/256 = 25% of the operand read the MMA already performs, an identity, not an estimate. Meet these and each route reaches its own roof, and the same datapath un-binds the conversion-bound sparse kernels; datapath width, ports, area and power are a Part-3 obligation, not a delivered count.
Figure 7: The central hardware co-design target, illustrated. (a) Today the residue conversion runs on the SIMT/integer pipe, competes with the fp8 MMA stream, round-trips through SMEM, and—for a large output—is repeated once per cluster: this is the deconstruction-λ floor (0.50 of the 473-TFLOPS roof on Rubin). (b) Option C performs the residue-decompose transform in the async copy engine, converting each fp64 tile in flight and depositing planes in MMA operand layout; deconstruction leaves the compute budget (cq →0 on the SIMT/INT pipe), re-splitting becomes free (λ-free), and—once the requirements of Table 9 are met—the emulated rate reaches each route’s own roof at today’s cluster, while the same datapath un-binds the conversion-bound sparse kernels. Option C does not zero conversion time, only the compute-pipe term: the transform must sustain 166–325 TMAC/s (1.3–2.6× the residual int8 pipe) and the deposit 55–65 TB/s of SMEM writes. The codesigned modulus sets (byte extraction plus a few mod-by-invariant multiply–truncates) keep it a fixed-function residue block rather than a general converter; everything else stays in software.
write traffic is exactly T /R = 64/256 = 14 of a read that already happens—an identity, not an estimate, and the strongest honest argument for Option C on the bandwidth axis. The input side is an L2 ask, not an HBM ask: 8/R bytes per FLOP is 12–15 TB/s, ≈17% of L2 service, because the slab panels are L2-resident under the fused schedule. The transform datapath is the large ask: at each route’s own roof it must sustain 166–325 TMAC/s, which is 1.3–2.6× the platform’s entire residual int8 tensor throughput. Earlier drafts described this block as “≈300 narrow MACs per copy engine”; that is a per-engine width, and it is not derivable without a clock and a copy-engine count, neither of which is a published quantity. We therefore state the aggregate rate, which is derivable, and leave the width, ports, area and power to the vendor—a Part-3 obligation, not a delivered count. What the software does buy is that the datapath is narrow: because our codesigned modulus sets reduce conversion to byte extraction plus a handful of mod-by-invariant multiply–truncates, it is a fixed-function residue block, not a general converter, and the moduli, the two-limb reduction, and the exact reconstruction all stay in software. Two further properties make it the preferred co-design choice under the present performance model. It is precedented: the copy engine already performs layout transforms (TMA swizzle, Blackwell decompress-on-copy). It also fixes sparse: the conversion-bound SpMV/GEMV kernels of Part 1 33
meet cq directly, so the same datapath un-binds them—one change closes both regimes, which no other row of Table 9 does. This is Plan A: the option prioritised by the present model, the only lever that reaches the roof at unchanged reach, and the only one that closes the sparse side as well. Table 9: Option C as a vendor requirement, not as a result. Every row is what a vendor must provision, evaluated per Rubin GPU at each route’s own roof under the 4×4 cluster of 642 tiles (reach R=256); none of it exists on shipping silicon. Script: verify/optionc_req.py. The three sets do not share a roof: with cq removed the ranking is set by α alone, so S (α=37) leads at 473, E (α=40) at 438, A (α=46) at 380 TF. Correction to earlier drafts: the deposit band was printed as 55–83 TB/s. That 83 folded two errors: it priced set A at set S’s 473-TF roof, which set A cannot reach, and it used pstore (A) = 45, which counted no squares in a set whose first modulus is 256 = 162 ; the correct value is 44 (§5.1, Table 4). Evaluated route-consistently at the corrected plane count the band is 55–65 TB/s. The two tmem-tile alternatives are likewise corrected: the multipliers here are what a vendor must provision, i.e. the worst case over the three sets (4× and 9×, both set A’s), against the per-route ranges 3.5–4.0× and 7.9–9.0× of §5.2; earlier drafts printed 3× and 6.75×, which are neither. Claim status: modeled requirement derived from params.py; no measurement, and no vendor commitment, is implied. requirement
what it is
S
E
A
each route’s own PFP8 /α 473 438 380 unique fp64 scalars the 1848 1709 1486 transform consumes, P/R as bandwidth (TB/s) 8/R B per FLOP, served from 14.8 13.7 11.9 L2 (88 TB/s) transforms per scalar (MAC) frozen-kernel narrow 176 136 112 multiply–adds transform rate (TMAC/s) the datapath ask 325 232 166 vs residual int8 × the platform’s 2.60× 1.86× 1.33× 125-TMAC/s pipe plane deposit (TB/s) SMEM writes, pstore /R B per 55.4 56.4 65.4 FLOP as a share of the read = T /R exactly, on every route 25% tensor-side operand read (TB/s) pstore /T B per FLOP 222 226 262 (unchanged by C) of which remote (DSM) (1−1/c)× the above 166 169 196 staging buffer (KiB/tile) planes of one 642 fp64 tile 120 132 176 (32 KiB in) latency / ports drain the staging buffer at the Part-3 obligation deposit rate area / power fixed-function residue block in Part-3 obligation the copy engine
route roof (TF) input scalar rate (Gscalar/s)
the alternatives, for the same dense effect C=64 (8×8) cluster cross-GPC DSM at 4× the fabric reach C=128 (16×8) cluster cross-GPC DSM, reach 683 TMEM tile 1282 TMEM SRAM 256 → 1024 KiB (4×) TMEM tile 1922 TMEM SRAM 256 → 2304 KiB (9×) faster L2 (Plan B) L2 service 4×HBM → 10–12×HBM
floor 0.50 → 0.92 floor → roof floor 0.50 → 0.92 floor → roof roof, dense only
Plan B, independent of the conversion datapath: a faster L2 with software-blocked materialisation. If a programmable copy-engine transform is judged too invasive, the datapath-independent fallback scales an existing structure—L2 service bandwidth—and lets software carry the complexity. Materialise the residue planes, but blocked: tile the output so each block’s planes fit the 126-MB L2 (at NB=2048 a 10242 block’s planes are almost exactly 34
126 MB), keep the row-panel’s planes L2-resident, and stream both operands to the tensor cores at the fp8 rate (λ=1 within the block). The materialised ceiling is then BL2 /(p/TILE): 0.60 of the 473-TFLOPS roof at 6×HBM, 0.79 at 8×, and the roof at ≈10–12×. The lone hardware change is L2 bandwidth—no new datapath, instruction, or fabric—and the blocking that satisfies the L2 capacity predicate is pure software, so this is the hardware-lightest path second to Option C. Two honest caveats: ≈10×-HBM L2 is a real bandwidth investment (today ∼3–5×), and, unlike Option C, it is dense-only—sparse kernels have no operand reuse to amortise the materialised feed, so it does nothing for SpMV/GEMV. Two further “scale an existing structure” options round out the menu: a larger TMEM tile (a 3.5–4.0× SRAM capacity bump → 0.92 on-the-fly, no materialisation, and it would let D run at nblk =2 rather than 3—though that halves D’s re-read charge without touching its roof, which is why D stays out of the menu either way), and a pure panel (NB) increase—free, but it reaches only the 0.50 floor, never the roof. The whole ladder is Table 5.
6
Projected Performance
Table 10 evaluates Eq. (13) across the Ozaki 2.5 ladder—S0, the shipping path; L1, the dp4a route of §5; the two-limb tensor route; and the codesigned sets of §5.1—on Rubin and GB300. Three readings summarise it. (i) The announced operating region—and where the gains need no hardware change. The n̄≤256 column is achieved today, single-cluster, with no hardware change: the switched dispatch reaches 243 TFLOPS at n̄=256 (≈2.4× the shipping path). Above 256, the reduced model projects n̄=512 rising from ≈200 to 364 (two-limb on the published set), 380 (all-byte), or 438 TFLOPS (hybrid/switched)—a conditional 1.8–2.2×. These larger-n̄ figures are convertonce envelopes: they are realised today for the tall/skinny and small-batch shapes of real applications (Table 11), where the switched dispatch is worth ≈1.6–1.9× (≈2× on the block-Krylov rows) over simple deconstruction with no hardware—but for a large square output they require the co-design of §5.2 and are otherwise clipped to the ≈235-TF floor. The practical reading is therefore positive: on Rubin, Ozaki 2.5 is already useful—near the crossover, only lightly clipped (0.89–0.94×, not the 0.50 floor)—for the block-Krylov, batched, and panel work that dominates real solvers; only the large-square-DGEMM headline waits on hardware. These are projections at ηfp8 =1; every loss the roof ignores subtracts from both numerator and baseline in ways only measurement can apportion. On GB300 the switched dispatch runs tensor A in the block-width band (its dp4a realisation within ≈7%), E and L1 through the middle, and the two-limb route on the roof (from n̄ ≈ 290): the 135-TFLOPS roof is reached from n̄ ≈ 220, with 95 modeled at n̄=128 (1.9× the shipping path there)—against GB300’s 1.4-TFLOPS native reference, a ≈17–96× band throughout the evaluated range n̄ ≥ 32, which is the Part 1 argument in its starkest form. (ii) One-shot large DGEMM. This is precisely the regime the deconstruction-λ floor governs (§5.2): a single large square output is multi-cluster, so the convert-once crossover of Eq. (13)—roof approached above n̄≈1200 for the shipping path—is realised only as an upper envelope, and the achieved rate plateaus at ≈0.50 of the 473-TFLOPS roof (Rubin), not at the roof. That convert-once reading holds only inside a single cluster or under materialisation; for the deployable on-the-fly schedule the binding term is the re-split floor, and lifting it to the roof (larger cluster, TMEM tile, or L2 bandwidth) is the co-design question of §5.2. (iii) Structurally sub-crossover workloads are where the routes change the picture qualitatively; they are treated in §6.1. Tensor-service concurrency: the two-endpoint bracket. The tensor-migrated routes charge their reduction GEMMs to the int8 tensor rate as if it were a concurrent capacity beside near-peak fp8 work; as §5 notes, that independence is not a proven property, and the honest statement is a bracket between two analysable endpoints. Under full independence the 35
Table 10: Reduced-model upper envelopes, software routes only (TFLOPS of useful fp64 2mnk work; ηfp8 =1; int8-capacity-aware effective rates; conditional projections, not delivered performance; generated from params.py). “Switched” dispatches the route by n̄ (Rubin: A below ≈410, E to ≈620, two-limb S above; GB300: tensor A below ≈180—its dp4a realisation within ≈7%—E to ≈210, L1 to the 135TFLOPS roof at ≈220, two-limb S from ≈290; Figure 5). The aspirational L1 count (≈7.5) is excluded from the envelope pending its reuse schedule. Hardware is deliberately absent: Option-C-class conversion hardware (Part 1) would remove the conversion term at every size, but on current GPUs these software routes are what is deployable. All route values are ideal-overlap envelopes; schedule eligibility (storage mode, shape, capacity) per §6.2; serialisation endpoints in the text. Claim status: modeled projection. n̄: achieved today (single cluster) n̄: convert-once envelope†
knee
32
256
512
1024
2048
n∗
Rubin (17.5-PFLOPS fp8; int8 cap 250 TOPS) S0: shipping path, cq =16 12 25 50 L1: dp4a direct, cq ≈10.3 19 39 77 Ozaki 2.5 on S (two-limb) 23 45 91 Ozaki 2.5 on A (all-byte) 30 61 121 Ozaki 2.5 on E (hybrid) 29 59 118 A at dp4a grain, cq ≈7.3 22 44 88 Switched (A/E/S by n̄) 30 61 121
100 155 182 243 235 176 243
200 310 364 380 438 352 438
400 473 473 380 438 380 473
473 473 473 380 438 380 473
1211 782 666 401 476 553 —
GB300 (5-PFLOPS fp8; int8 cap 166 TOPS) S0: shipping path, cq =16 12 25 50 L1: dp4a direct, cq ≈10.3 19 39 77 Ozaki 2.5 on S (two-limb) 15 30 60 Ozaki 2.5 on A (all-byte) 24 47 95 Ozaki 2.5 on E (hybrid) 20 39 78 A at dp4a grain, cq ≈7.3 22 44 88 Switched (A/E/L1/S by n̄) 24 47 95
100 135 121 109 125 109 135
135 135 135 109 125 109 135
135 135 135 109 125 109 135
135 135 135 109 125 109 135
346 223 287 147 205 158 —
software route
64
128
†
Clipping status. Rates are tabulated for square outputs, for which n̄=edge and single-cluster ⇔ edge≤256. Columns n̄≤256 thus fit one cluster and are achieved today (throughout, in the §1 sense: modeled on shipping silicon with no hardware change, not measured), irrespective of clipping (a non-square shape at the same n̄ is instead governed by the size-weighted λeff of §5.1). Columns n̄>256 are convert-once envelopes: achieved today for tall/skinny or multi-panel shapes, whose λeff stays O(1) (≈1–1.5—the real-application regime of Table 11, where they are realised), but for a large square output (λeff =edge/256) they require the co-design of §5.2 and are otherwise clipped to the deconstruction-λ floor (≈235 TFLOPS on Rubin, Table 2). GB300 is the exception: its floor equals its 135-TFLOP roof, so its n̄>256 cells are roof-bound (achieved today), not clipping-limited. Every Rubin rate in this paper is thus either fundamental (achieved today) or clipping-limited (a convert-once envelope that needs the hardware of §5.2 for large square DGEMM); the two are distinguished wherever a number appears.
tensor-migrated routes deliver the Table 3 values. Under complete serialisation on a shared tensor service the delivered rate is the harmonic composition Pser = (1/PF + 1/PI )−1 of the fp8-side and int8-side rates: at n̄=512 on Rubin this gives ≈206 (S), ≈228 (A), and ≈227 (E) TFLOPS. The dispatch’s floor, however, is int8-tensor-contention-free: the pure-dp4a routes use no tensor int8 capacity at all, so even a fully serialised tensor service leaves A-at-dp4a at ≈352 and L1 direct at ≈310 TFLOPS at n̄=512 on Rubin. The honest software-route bracket at n̄=512 is therefore [352, 438] for the switched dispatch and [310, 364] for the published set (dp4a versus two-limb), degrading the conditional uplift over the published 200 from ≈2.2× to a floor of ≈1.76×; at n̄=256 the bracket is [176, 243]. On GB300 the dp4a route reaches the full 135-TFLOPS ceiling from n̄ ≈ 223, so the GB300 conclusions are int8-tensor-contention-free. Between the endpoints we interpolate as Ttensor = max(TF , TI ) + θ min(TF , TI ), θ ∈ [0, 1] (θ=0 independent, θ=1 serialised); measuring θ(shape, occupancy) joins the validation plan of §7. SIMT–fp8 overlap: the second serialisation axis. Avoiding int8-tensor contention is not the whole story: the dp4a routes still assume their SIMT conversion work overlaps the fp8 MMA stream. We parameterise that overlap by ρsimt,fp8 ∈ [0, 1] (0 ideal overlap, 1 complete 36
serialisation on the shared issue path). At complete serialisation the Rubin n̄=512 endpoints are the harmonic compositions (1/352 + 1/380)−1 ≈ 183 TFLOPS for A-dp4a and (1/310 + 1/473)−1 ≈ 187 for L1. Read those numbers as an equal-overlap-degradation sensitivity, not an unconditional floor—it conditions on both paths sharing one overlap coefficient—and keep its three ingredients separate: (i) the route endpoints (ideal 352/310, serialised 183/187 at n̄=512); (ii) the announced 200, whose size and algorithm provenance are unknown; and (iii) the matched-degradation ratio. For (iii): complete SIMT–tensor serialisation collapses the shipping baseline too, because S0 embeds the same overlap assumption—the memo’s Pint normalisation and the shipping path’s own SIMT conversion presume it—so the consistently serialised S0 is (1/200 + 1/473)−1 ≈ 141 at n̄=512 and ≈83 at 256, and the like-for-like ratio is ≈1.30× (512) and 1.46× (256). Comparing serialised routes against the announced 200 mixes worlds and is not claimed. Structurally, a tensor MMA occupies a shared issue slot only once per many tensor-busy cycles, so the shared-issue serialisation mechanism is weak—an argument, not a measurement, and ρsimt,fp8 joins θ in the paired-rate measurements of §7.
6.1
Application-Geometry Scenarios: Measured Call-Shape Traces
The value of the switched dispatch is best seen against the workloads that actually occupy the sub-crossover region, and against the alternative they would otherwise fall back to: Rubin’s native fp64 path, illustratively ≈30 TFLOPS of blocked DGEMM (Figure 3, flat line; to be measured). The geometry matters: rank-nb trailing updates have large n̄ and were never the problem; the structurally bounded class is where one output dimension is a block width— supernodal/frontal update panels in SuperLU_DIST- and MUMPS-class sparse direct solvers [8, 1], panel-internal factorisation kernels, LOBPCG and block-Krylov eigensolvers [7] (block widths 16–128), and batched small GEMMs in tensor contraction. Under the switched dispatch the reduced model projects, at n̄ = 64/128/256/384: 61/122/243/364 TFLOPS on Rubin, versus 25/50/100/150 on the shipping path—a uniform ≈2.4× across the band, with the emulationversus-native crossover pushed from n̄ ≈ 80 down to n̄ ≈ 32 against the illustrative 30-TF baseline (≈35 against the official 33-TFLOPS specification [20])—at or below typical block widths. On GB300 the same band runs 47/95/135/135 versus 25/50/100/135 (≈1.9× at block widths 64–128, converging at the 135-TFLOPS roof beyond the shipping knee of ≈350)—and with native fp64 at ≈1.4 TFLOPS there is no crossover to push within the evaluated range: every plotted route exceeds the native reference for n̄ ≥ 32, and the dispatch question is only which route. Measured call-shape traces. Whole-application gains depend on each code’s mix of GEMM geometries. Rather than assume mixes, we measured them: using an LD_PRELOAD interposer (traces/blas_shim.c, supplied in the artifact with every runner script), we captured the (m, n, k) of each dgemm/dsyrk call from six instrumented configurations across four library/application classes, all running unmodified library implementations on author-constructed inputs—SciPy LOBPCG on a 483 Laplacian at block widths 64 and 128; multifrontal sparse LU (UMFPACK) on a 7002 convection–diffusion operator; coupled-cluster tensor contractions (PySCF CCSD on benzene in cc-pVDZ, capped at six iterations); and netlib LAPACK blocked LU and QR at n=4096. The fixed captured traces measure geometry for the stated software stack; replaying them through the route model produces GEMM/SYRK-portion scenario projections, not application benchmarks or executions of an Ozaki-2.5 kernel (Table 11, Figure 8; work-weighted harmonic composition over Eq. (13), the same machinery as Table 10). Three measured facts stand out. First, LOBPCG at b=64 places 100% of its GEMM work at n̄ ≤ 256 (median 128: the 2mn/(m+n) → 2b geometry, measured), and netlib blocked QR places half its BLAS-3 work at n̄ ≈ 64—the block-width-bounded class is real and load-bearing. Second, multifrontal fronts are fatter than commonly assumed: UMFPACK’s median GEMM sits at n̄ ≈ 670 (quartiles 380–850), squarely in the crossover band of §4—real sparse-direct 37
work lives exactly where the deconstruction term bites, which independently motivates the size sweep. Third, sequential supernodal SuperLU (SciPy’s splu) issues essentially no BLAS-3 at all (its panel updates are dgemv-based), a reminder that the BLAS-3 supernodal geometry this paper targets is the multifrontal/distributed class. Amdahl fractions for non-GEMM work further dilute whole-run numbers, and the mixes below are GEMM-portion composites only; the native columns compare against flat illustrative references and should be read as order-of-magnitude bands, not predictions. Composite model conventions. The geometry of Table 11 is measured; the composites are model projections over it, computed as follows. Each record is assigned its eligible storage schedule (split-k / fused-λ / L2-materialised / HBM-materialised planes at 2p bytes, §6.2) and its route rate is capped by the record’s memory roof at its additive traffic Q0 + Qplane with Q0 = 8{k(m+n) + mn} (Eq. (22)); the reconstruction tax of Eq. (11) is applied on the dispatched host; and the shipping baseline S0 is capped by the same memory roof and tax, so the comparison is like-for-like. dsyrk records are charged as one unique operand (the Gram factor A of C = AAT ): triangular work m n k, a single converted operand nuniq = mk in the conversion and plane ledgers, and n̄ = nC ; the transposed operand view is a marshalling, not a second conversion; the output C is charged a single write under the overwrite convention βBLAS =0 that Q0 uses throughout (a nonzero βBLAS would add a read of the triangular C, which the interposer does not currently record). (These traces contain no dsyrk-tagged records— the LOBPCG Gram products were interposed as GEMM—so this convention is stated for completeness and leaves the present composites unchanged.) Quantiles and composite speedups P P are FLOP-weighted, in the harmonic form P = W / W /P (n̄ ) with speedup comp i i i i i P P base route S= / ; per-record predictions ship as traces/sched_*.csv. i Wi /Pi i Wi /Pi The reading is not that every solver triples, but that the measured geometry confirms the structural claim: the workloads for which emulated fp64 was least usable—LOBPCG and QR panels, whose traced work is block-width-bounded—are exactly where the software-only switched dispatch is modeled to be worth multiples, on current GPUs, with no hardware change; and the traced multifrontal fronts populate the crossover band itself. Relative to the earlier draft, the schedule/memory-roof filter and the on-the-fly deconstruction-λ correction together move several composites materially—on Rubin, multifrontal 2.07 → 1.19× and QR 2.27 → 1.60×; on GB300, LOBPCG 1.90/1.64 → 1.67/1.59×, multifrontal → 1.19×, and dense LU → 0.87× (below unity: at n̄=3200 the rank-64 trailing update is bandwidth-bound and the residue reconstruction makes emulation marginally slower than the shipping baseline there)— because the multifrontal small-front tail and the GB300 tall-skinny block shapes (operational intensity ≈8) are bandwidth-bound, where the emulated and shipping paths nearly coincide at the same memory roof; the GB300 LOBPCG composites are capped by exactly that roof. The GB300 block sharpens the point from the other side: per-route gains are more modest (the roof is nearer), but with a ≈1.4-TFLOPS native reference every traced workload is an order-ofmagnitude argument for emulating at all—led, below the roof, by pure-SIMT routes that any CUDA kernel can implement today. The final three columns answer, per workload, the question the co-design of §5.2 raises: if the deconstruction floor were lifted all the way to the compute roof (473 TFLOPS on Rubin, 135 on GB300), what would each of these workloads then get—against the shipping baseline, and against the native fp64 unit a user would otherwise run on? The answer is bounded well short of the roof—28–74% of it on Rubin, 36–95% on GB300—because once deconstruction is free the binding constraint on this measured geometry is the memory roof, not the arithmetic one. The additional headroom the co-design buys over Ozaki 2.5 as modeled here is therefore ≈2.5–3.3× on Rubin, where the two roofs are far apart, and only ≈1.1–1.5× on GB300, where they are not. Two things follow. The co-design case is a Rubin-class argument, not a general one; and even granting the co-design in full, these traced applications would remain bandwidth-limited, 38
Table 11: Measured call-shape traces and their model-projected GEMM-portion composites under the switched dispatch (both platforms). Geometry is measured (BLAS-3 interposer; hardware-independent); the composites are model projections over that geometry under the per-record schedule-, reconstruction, and memory-aware model detailed in the preceding paragraph. “work n̄≤256” = fraction of traced GEMM flops at n̄ ≤ 256. Claim status: modeled projection over measured geometry. Composites are reported at the conservative SIMT-Garner reconstruction cost floor; the ≈r codesign target of §3 would raise the small-k Rubin native multipliers as noted. The final three columns are the counterfactual of §5.2: the same per-record schedule, and the same memory roof and plane-feed caps, but with the deconstruction floor lifted to the compute roof (473 TF on Rubin, 135 TF on GB300). It is a ceiling bought by hardware that does not exist, not a delivered rate. Its two multipliers are that ceiling over the shipping baseline and over the native fp64 unit—the same pair as the delivered columns to their left, so the two cases can be read against one another rather than only the co-design case against today’s software; delivered is the composite rate so obtained and the fraction of the roof it represents. Because both pairs are ratios of the same two rates to the same two references, each pair carries the identical headroom factor, which is the quantity the next sentence reports. No row reaches 100%—once deconstruction is free the traced geometry is bandwidth-bound, so the additional headroom the co-design buys over Ozaki 2.5 as modeled today is ≈2.5–3.3× on Rubin but only ≈1.1–1.5× on GB300, whose memory roof is much nearer its compute roof. as modeled today workload (traced config.)
med. n̄ (q25–q75)
work delivered n̄≤256 vs ship vs native vs ship vs native (% of roof) dominated by
Rubin (native ≈30 TF illustrative; roof 473 TF) 128 (64–128) 100% ≈1.94× LOBPCG b=64 (483 Laplacian) LOBPCG b=128 256 (128–256) 100% ≈1.83× (same operator) Multifrontal LU 665 (383–850) 16% ≈1.19× (UMFPACK, 7002 ) CCSD contractions 800 (353–1302) 19% ≈1.69× (benzene/cc-pVDZ) Dense LU (netlib, 3200 (2560–3712) 1% ≈1.14× n=4096) Dense QR (netlib, 64 (63–3232) 50% ≈1.60× nb =32) GB300 (native ≈1.4 TF; roof 135 TF) 128 (64–128) LOBPCG b=64 (483 Laplacian) LOBPCG b=128 256 (128–256) (same operator) Multifrontal LU 665 (383–850) (UMFPACK, 7002 ) CCSD contractions 800 (353–1302) (benzene/cc-pVDZ) Dense LU (netlib, 3200 (2560–3712) n=4096) Dense QR (netlib, 64 (63–3232) nb =32)
if co-design reaches roof
≈2×
≈6.00×
≈6×
176 TF (37%) block
≈4×
≈6.00×
≈12×
352 TF (74%) block
≈1×
≈3.80×
≈4×
132 TF (28%) frontal
≈4×
≈4.30×
≈10×
288 TF (61%) contrac-
≈3×
≈3.68×
≈10×
299 TF (63%) trailing
≈2×
≈5.19×
≈6×
172 TF (36%) panel
width width panels tions updates factors
100% ≈1.67×
≈35×
≈2.18×
≈46×
64 TF (47%) block
100% ≈1.59×
≈67×
≈2.18×
≈91×
128 TF (95%) block
16% ≈1.19×
≈29×
≈1.38×
≈34×
48 TF (36%) frontal
19% ≈1.41×
≈58×
≈1.62×
≈66×
93 TF (68%) contrac-
≈0.87×
≈50×
≈1.34×
≈78×
109 TF (80%) trailing
50% ≈1.41×
≈33×
≈1.89×
≈45×
63 TF (46%) panel
width width panels tions 1%
updates factors
which is the honest reason we book the deconstruction floor as the first target rather than the only one. These are scenario projections over measured geometry, not application benchmarks; replaying the to-be-released traces through measured kernels is step (vi) of §7.
39
share of GEMM work (%)
LOBPCG b = 64 (SciPy, 483 Laplacian) 100
737 calls 667 GF total ñ 1/2 = 128 100% at n̄ ≤ 256
80
Multifrontal LU (UMFPACK, 7002)
LOBPCG b = 128 737 calls 2.7 TF total ñ 1/2 = 256 100% at n̄ ≤ 256
17,120 calls 12 GF total ñ 1/2 = 665 16% at n̄ ≤ 256
1024 4096
16
60 40 20 0 16
64
256
1024 4096
16
share of GEMM work (%)
CCSD contractions (PySCF, benzene/cc-pVDZ) 100 80
40,611 calls 749 GF total ñ 1/2 = 800 19% at n̄ ≤ 256
64
256
64
256
1024 4096
Dense QR (netlib dgeqrf, nb = 32)
Dense LU (netlib dgetrf 4096) 4,095 calls 45 GF total ñ 1/2 = 3200 1% at n̄ ≤ 256
248 calls 91 GF total ñ 1/2 = 64 50% at n̄ ≤ 256
60 40 20 0 16
64
256
1024 4096
n ̄ bucket (lower edge)
16
64
256
1024 4096
n ̄ bucket (lower edge)
16
64
256
1024 4096
n ̄ bucket (lower edge)
Figure 8: Measured call-shape distributions (share of traced GEMM/SYRK work per log2 n̄ bucket) for the six instrumented configurations, from the interposer logs (traces/, in the artifact). The shaded band marks n̄ ∈ [128, 512]. LOBPCG and QR-panel work is block-width-bounded; multifrontal fronts and CCSD contractions straddle the crossover band; dense LU trailing updates sit far right. Geometry is hardware-independent; only the projections built on it are model-conditional.
6.2
Storage Modes: Where the Planes Live, and What Each Mode Charges
The projections above lean on two disciplines at once—O1’s convert-once accounting and a fused traffic model in which the residue planes generate no HBM traffic of their own—and these are not simultaneously free: a plane that is never written to memory may have to be converted more than once, and a plane converted exactly once must live somewhere. We make the coupling explicit. For an m×k by k×n product with nuniq = k(m+n) unique operand elements and p stored plane bytes per scalar (30/44/33 for S/A/E), the storage modes are: mode F (fused)—the planes live in the SMEM of the owning CTA or thread-block cluster, conversion multiplicity λ ≥ 1 (a tile shared by several clusters may be converted by each), no plane HBM traffic; mode M (transient materialisation)—the planes are staged through a workspace, L2-resident or in HBM, and converted exactly once; and mode P (persistent)—the planes are retained across calls (operand-stationary O1), each call still reading p bytes per touched scalar. The rest of this section derives, mode by mode, the bound each of these charges. Because those bounds are not evaluated in isolation but as a per-call selection—the engine tries every eligible schedule on every route and keeps the best—we state that selection once, as Algorithm 4, before deriving its ingredients. It is the exact content of rate_rec in params.py: every projection in this paper, every composite in Table 11, and every per-record tag in traces/sched_*.csv is an evaluation of it. That last sentence is a claim about the correspondence between a printed algorithm and a program, which is exactly the kind of claim this paper has already had to retract once, so it is a test and not a promise: verify/alg4_fidelity.py transcribes the algorithm as printed, independently of the engine, and requires the two to agree on both the 40
winning rate and the winning tag at every point of a grid that includes the tall/skinny corner where the ledger omissions below used to hide. Reading the algorithm alongside the derivations that follow is how those omissions became visible in the first place. Mode M through HBM: the plane-only obstruction bound. Transient materialisation through an HBM workspace mandates the plane write and at least one read: the plane traffic is Qplane = 2p nuniq bytes, so Puseful ≤
Bmem n̄ , 2p
n̄ =
2mn , m+n
(21)
against 2mnk useful flops. (An earlier draft charged only the read; the write is not optional, and the crossovers below are the corrected values, doubled accordingly.) The bound stops binding above the plane-only crossover n̄mat = 2p·(set roof)/Bmem : on Rubin (22 TB/s) 1290/1522/1312 for S/A/E, on GB300 (8 TB/s) 1014/1196/1031. This is a necessary obstruction bound only: below n̄mat HBM materialisation cannot support the set roof whatever the compute engines do, but above it the plane-only obstruction merely ceases—roof attainment is a stronger, additive condition stated next. The additive full-traffic bound and the roof-attainment crossover. The obstruction bound above isolates the plane term; the schedule must in fact move the plane traffic and the call’s own operand/output traffic across the same bus. Charging both additively—plane bytes Qplane plus Q0 = 8{k(m+n) + mn}—the memory roof is Bmem Puseful ≤ , 8+c 4 + n̄ k
(
c=
2p p
both operands materialised, one operand (fused-hybrid),
(22)
and roof attainment requires n̄ above the additive crossover n̄add = (8 + c + 4 ⊮sq ) Proof /Bmem (with Proof the set roof and ⊮sq =1 for square shapes). For full materialisation (c = 2p) this is 1462/1695/1472 (large-k) and 1548/1764/1551 (square) on Rubin, and 1149/1332/1156 and 1216/1386/1219 on GB300, for S/A/E. The fused-hybrid one-operand crossovers (c = p) are lower, at 817/917/815 (Rubin) and 642/720/641 (GB300), which is why the switched dispatch attains the roof earlier than the full-materialised threshold would suggest. These one-operand values use c = p, hence the additive term (8+p)/n̄, which is the square case; for a general shape the additive memory cost per useful FLOP is Q0 + 2p k s 8 4 p = + + , 2mnk n̄ k ℓ
s = min(m, n),
ℓ = max(m, n)
(the per-record engine uses the exact 2p k min(m, n) plane term), reducing to (8+p)/n̄ only when m = n. We report both: n̄mat (plane-only, necessary) and n̄add (additive, roof-attainment). A model defect this exposition exposed, and its repair. Writing the additive bound down carefully made it possible to audit the engine against it, and the audit failed. Earlier revisions of this paper asserted that the per-record engine had always capped each materialised mode by its additive traffic Bmem W/(Q0 +Qplane ), and that an audit of the schedule tags showed no record taking an HBM-materialised branch. Both statements were false, and we correct them here rather than quietly restate them. The HBM-materialisation branch charged only the kslab reduction traffic into its additive roof, together with a plane-feed cap of p/TILE bytes per useful FLOP; it never charged the c = 2p write-plus-read of Eq. (22) at all. Those two terms cross at n̄ = 2 TILE = 128, so on square shapes the feed cap is the binding one and the omission is invisible—which is why it survived—but on tall-and-skinny shapes it is not. 41
Algorithm 4 Select-Schedule(π, ν, m, n, k) — the per-call storage-mode and route selection. π is the platform, ν a route of the menu (§5.1). Every schedule is capped by its own compute rate and by its additive memory traffic Q0 + Qplane and by any bandwidth it must be fed at; the best surviving rate wins. The envelope quoted throughout the paper is maxν of this function. 1: input: platform π, route ν; shape m, n, k; set S(ν) with p stored plane bytes/scalar and
r moduli; B = Bmem (π), BL2 = 4B, CL2 = 126 MB, TILE = 64, C = 16 CTAs/cluster, R∗ = 256, Kslab = 4096 2: W ← 2mnk; Q0 ← 8{k(m+n) + mn}; nuniq ← k(m+n) ▷ syrk: W ←mnk, one operand, mn terms halved 3: Sk ← ⌈k/Kslab ⌉; Qred ← 8 r mn (Sk −1) ▷ split-k residue partials, written and read 4: Tm , Tn ← ⌈m/TILE⌉, ⌈n/TILE⌉; P ← ∅ ▷ P: (rate, tag) candidates ▷ mode F, one cluster: DSM-shared, λ = 1 ρ ← ComputeRate(ν, λ=1) 7: if Sk = 1 then ▷ Qred =0: one slab, no partials 8: add (min{ρ, BW/Q0 }, grid-fused-nosplit) to P 9: else if 4 r mn > CL2 then ▷ partials spill: charge them to HBM 10: add (min{ρ, BW/(Q0 +Qred )}, splitk) to P 11: else ▷ partials stay resident: charge them to L2 12: add (min{ρ, BW/Q0 , BL2 W/Qred }, splitk) to P 13: else ▷ multi-cluster: mode F pays λ, mode M pays traffic 14: λA , λB ← ⌈n/R∗ ⌉, ⌈m/R∗ ⌉; λeff ← (mλA + nλB )/(m+n) ▷ size-weighted multiplicity; R∗ from Eq. (19) 15: ρλ ← ComputeRate(ν, λeff ); ρ1 ← ComputeRate(ν, 1) 16: Qre ← 8{mk(λA −1) + kn(λB −1)} ▷ mode F: fp64 operand re-reads, not planes 17: add (min{ρλ , BW/(Q0 +Qred +Qre )}, on-the-fly) to P 18: if p k min(m, n) ≤ CL2 then ▷ hybrid: small operand in L2, large on the fly p 19: add (min{ρλ , BW/(Q0 +Qred ), BL2 /Λ(p k min(m, n), 2 TILE )}, hybrid-L2) to P 5: if Tm Tn ≤ C then 6:
if p nuniq ≤ CL2 then ▷ mode M in L2, Eq. (23): no plane HBM traffic p 21: add (min{ρ1 , BW/(Q0 +Qred ), BL2 /Λ(p nuniq , TILE )}, mat-L2) to P 22: else ▷ mode M in HBM: planes cross the same bus 23: Qplane ← Qred + 2p nuniq ▷ c = 2p, Eq. (22): write and read 24: add (min{ρ1 , BW/(Q0 +Qplane ), B TILE/p}, mat-HBM) to P ▷ both terms: the feed cap alone is looser below n̄=128 25: return max P, with the maximising tag ▷ the tag emitted to traces/sched_*.csv 20:
26: function Λ(plane bytes Z, feed ϕ) 27:
return Z/W + max{ϕ, Z/W } n̄=64
▷ L2 bytes/useful FLOP, L2-resident workspace ▷ write once, read ≥ once: ϕ is under one pass below
28: function ComputeRate(route ν, multiplicity λ)
Proof ← Pfp8 /α; n̄ ← 2mn/(m+n) ▷ roofs tabulated in Table 3 if co-design roof assumed then return Proof ▷ overlay: co-design column of Table 11 31: σ ← λ cres ▷ SIMT ops/FLOP: deconstruction residual + Garner q r/n̄ + Nγ (r)/2k 32: δ ← λ aint8 /n̄ ▷ int8 MACs/FLOP; aint8 , cres q from Table 3 33: return min{Proof , ISIMT /σ, (PI8TC /2)/δ} ▷ Eq. (16), Eq. (14); an unloaded pipe contributes ∞ 29:
30:
42
In the released trace set, 276 CCSD records (all Rubin, all route Et , all at n̄ ≈ 42, 6.4% of that application’s FLOPs) took the branch, and the worst of them, (m, n, k) = (1953, 21, 4325), was credited 38.2 TF where Eq. (22) permits 12.3 TF—a 3.1× escape from this paper’s own obstruction bound. A second instance of the same omission class sat in the L2-resident branch: §6.2 says the planes there are “written and read at cache rates”, but the engine charged the streaming feed alone, which for n̄ < 64 is less than a single pass over the workspace—6167 of 8121 L2-resident records. The engine now charges c = 2p additively on the HBM branch (in addition to, not instead of, the feed cap) and write-plus-at-least-one-pass on the L2 branch. The repair is contained and moves every affected number the conservative way. All 276 HBM-materialised records fall back to the on-the-fly schedule, so the tag histogram of traces/ sched_*.csv now genuinely contains no HBM-materialised record—a property that is asserted by test rather than by prose (verify/paper_numbers.py), precisely because asserting it by prose is what failed here. A further 68 records move from L2-resident materialisation to onthe-fly. Exactly one composite in Table 11 changes: CCSD on Rubin falls from ≈1.73× to ≈1.69× over shipping (116 → 114 TF delivered, co-design headroom ×2.49 → ×2.54). Every other application, both platforms, the projection and co-design tables, and the whole co-design column of Table 11 are unchanged to the printed precision. Mode M in L2: the capacity predicate. The workspace need not touch HBM at all. Whenever the whole plane workspace fits in L2, p k (m+n) ≤ CL2 ,
(23)
the planes are written and read at cache rates and generate no plane HBM traffic whatever. We assume CL2 = 126 MB—a stated, checkable parameter, not a measured one—under which the square-problem reach is k ≤ 1449/1196/1381 for S/A/E; L2-resident materialisation thus covers precisely the moderate-k region in which the HBM bound of Eq. (21) would otherwise bind. This “no HBM plane traffic” claim is an L2-capacity-eligible hypothesis, testable but unmeasured: it presumes usable (not merely nominal) capacity net of competing footprints, a producer–consumer lifetime that keeps the planes resident between write and read, and a write-back policy that does not spill them. As a safety-factor sensitivity, halving the effective capacity to Ceff = 12 CL2 lowers the S-set square reach from 1449 to 1024, keeping the L2-resident materialisation predicate honest under a 2× residency haircut (the single-cluster grid-fused reach is the separate, smaller Tm Tn ≤ C bound of §6.2, ≈256 square). Grid-tiled split-k (single-cluster, unmeasured). When the output fits one thread-block cluster—Tm Tn ≤ C output tiles, with Tm = ⌈m/64⌉, Tn = ⌈n/64⌉ and C the CTAs per cluster (8 portable, 16 opt-in on Blackwell), i.e. roughly max(m, n) ≲ 256 for square outputs—the cluster’s CTAs share each k-slab’s converted A/B residue planes through distributed shared memory (DSM): every operand element is converted exactly once (λ = 1), the sharing stays on-chip, and there is no HBM or L2 plane storage at all. This covers Gram-type products (C = ATA) and the small-front shapes of the traces. Its ledger is a proposed, resource-eligible, unmeasured schedule: with Sk k-slabs each held by a CTA, the per-slab partials are kept in the residue domain (exact, no rounding on reduction) and combined through an L2-resident reduction tree. Larger outputs span multiple clusters, which cannot share converted planes across clusters without either repeated conversion or an L2-materialised plane workspace. The repetition is cluster-granular, not tile-granular: a cm × cn cluster spans 64cm × 64cn output elements, which is 256×256 at the 4×4 grid of C=16 and 256×128 at the portable C=8 grid 4×2 (Eq. (19)), so an operand block is re-split once per cluster that consumes it— λA = ⌈n/64cn ⌉ (the A panel reused down the n axis), λB = ⌈m/64cm ⌉, i.e. ⌈n/256⌉ and ⌈m/256⌉ at C=16— for a size-weighted per-element multiplicity λeff = (mλA + nλB )/(m+n) (Appendix figures; params.py: _lam_eff). For a square output this is E/256; those shapes take the on-the-fly 43
(or, for the smaller operand, L2-materialised) schedule above, not this branch. The k-slab reduction (Sk > 1) is orthogonal to this spatial-cluster question and is charged under either tiling as L2 service traffic ≈ 2 · 4r mn (Sk −1) bytes, i.e. ≈ 4r/Kslab bytes per useful FLOP (useful work 2mnk with k = Sk Kslab ; the earlier ≈ 8r/Kslab dropped the factor of two in the work). This traffic is L2-resident and generates no HBM plane traffic, but it is not below 1% of Q0 : as a fraction of Q0 it is ≈ rn/(2Kslab ), e.g. ≈ 150% for a width-n=1024 square output at Kslab =4096. What must instead be checked is the L2 service rate: at the achieved emulation rates the worst traced record demands ≈ 2.44 TB/s of L2 service (the CCSD (441, 441, 4324) record on the E route; ≈ 6.9 TB/s were any record run at the full 473-TFLOP S roof). That BL2 = 4Bmem (4× HBM) is a model assumption: NVIDIA documents the 126-MB L2 capacity, not a 4× sustained reduction bandwidth. Under 1 × /2 × /4× L2:HBM ratios (22/44/88 TB/s on Rubin, 8/16/32 on GB300) the worst-case demand stays under the cap in every case, so no traced record is L2-service-bound within the assumed 4× service model—the conclusion does not depend on the boost; the companion engine charges this as an explicit L2 service cap (params.py: B_L2, K_SLAB). The engine gates the λ = 1 branch on the single-cluster predicate Tm Tn ≤ C (params.py: CLUSTER_CTAS=16) and tags Sk = 1 records grid-fused-nosplit to distinguish pure spatial tiling from actual k-splitting; the launched-CTA count is Tm Tn Sk , not merely the Tm Tn spatial tiles. The residue accumulators live in per-CTA TMEM (PTX: a 512-column × 128-lane logical array, 256 KiB), and they cost 4Nacc bytes per output, not 4r: Nacc is counted per modulus (§5.1) and is 30/26/30/34 for S/E/A/D, so an unblocked 642 tile asks 480/416/480/544 KiB and nothing fits. What fits is the blocked schedule of (20): at nblk =2 the live counts are 15/14/16/18 tiles = 240/224/256/288 KiB, so S, E and A fit (A with no headroom at all) and D does not. D fits only at nblk =3 (12 tiles, 192 KiB), which doubles the fp64 slab re-read charge from ≈1/3 to ≈2/3 of l2; that, together with its 337 TF roof—the lowest in the menu—is why D is not carried as a design point (verify/blocking.py). This distinguishes the MMA K-chunk Kmma =64 from the split slab Kslab ≥ 4096. Fused-hybrid with integer λ, and the surcharge predicate. For k-large shapes outside the L2 predicate with one small output dimension, the stream-side operand is converted in the owning cluster with integer multiplicity λ = ⌈min(m, n)/256⌉ (the 256-element square reach of a 16-CTA cluster), delivering roof · min 1, n̄/(λ · knee) under ideal overlap. The smallside operand’s planes are materialised once—in L2 if p k min(m, n) ≤ CL2 , else in HBM at 2p k min(m, n) bytes—and the HBM surcharge relative to the call’s full traffic Q0 = 8{k(m+n)+ mn} is 2p k min(m, n) δ = . (24) Q0 For square shapes δ = p/12 = 2.50/3.75/2.75 for S/A/E—the surcharge dwarfs the call traffic— and requiring δ ≤ 0.2 in the k ≫ n limit forces an aspect ratio m/n ≥ 36.5/55.25/40.25: the hybrid is a genuinely tall/skinny schedule, eligible only where δ ≤ ε or its small planes are L2-resident. An earlier draft claimed a blanket “≤10–20% of call traffic” surcharge; that claim is withdrawn for general shapes—it was a tall-skinny value. L2 traffic, plane-feeding, and why the design generates on the fly. The residue planes must reach the tensor cores somehow, and there are two ways. Materialise them—write the r planes to memory once and stream them back—or generate them on the fly, splitting each operand block into planes in registers/TMEM immediately before the MMA that consumes it (the fused deconstruction schedule of the companion Part 1 [11]). Materialisation looks cheaper (convert once) but is feed-starved: at the 64×64 TMEM-limited output tile the operand reuse is 64-fold, so streaming both materialised operands to the cores costs p/TILE bytes per useful FLOP, and to hold set S’s roof that feed must sustain p/TILE · 473 = 222 TB/s. The companion figures for E and A, 226 and 262 TB/s, are each taken at that set’s own 44
roof (438 and 380 TF), not at S’s 473: the three sets do not share a roof, and pricing one set’s plane count at another set’s roof is exactly the error corrected in Table 9. L2 at the assumed 4×HBM supplies 88 TB/s and HBM 22; a materialised schedule is therefore capped at BL2 /(p/TILE) = 188/171/128 TFLOPS (S/E/A, L2-resident) or BHBM /(p/TILE) = 47/43/32 (HBM-resident)—both below the on-the-fly rate derived next. These three are the feed term alone; the per-record engine additionally charges the one-off write of the planes into L2, so its realised L2-resident cap sits a little under them (hence 186/169 rather than 188/171 at the small-NB end of Table 8), which tightens the inequality in the same direction. So the design generates on the fly, and the residue planes never occupy L2 or HBM. One consequence is worth stating plainly, because it reverses a claim of our own earlier draft: the emulated throughput is now insensitive to the L2 service bandwidth (the planes are not L2-resident), so the BL2 assumption that a previous version made load-bearing no longer conditions the result (the trace composites move by ≤ 1 native-× across BL2 = 4/2/1×HBM). What on-the-fly generation does cost is re-splitting. A single 16-CTA cluster of 64 × 64 tiles covers a 256 × 256 output (Tm Tn ≤ C), inside which every operand block is split once and shared through distributed shared memory (λ = 1). A larger output spans several clusters, and—planes being on-chip and unshared across clusters—each operand block is re-split once per cluster that consumes it: A[m, k] (reused down the n axis) ⌈n/256⌉ times, B[k, n] ⌈m/256⌉ times, for a size-weighted multiplicity λeff = (m⌈n/256⌉ + n⌈m/256⌉)/(m+n). For a square output of edge E this is E/256; for a tall/skinny one it stays O(1)—bounded, not unity: the large operand is split once (λA =1) and the small one many times but contributes little, giving 1.25/1.50 at LOBPCG’s 110,592 × 64/128. Those values are for the square 4×4 dispatch the engine models, which is a conservative upper bound rather than a property of the shape: a cluster shaped to the output (cm cn ≤ C; here 16×1 and 8×2) keeps every CTA productive and lowers them to 1.06/1.25 at no hardware cost, so the tall/skinny entries below are quoted against the pessimistic dispatch. The split-k reduction service—the term a previous round flagged—is by contrast negligible: its worst demand across all trace records is 2.44 TB/s, far under 88. The engine charges all of this (params.py: rate_rec, on-the-fly vs. materialise, best wins). How the limit manifests, kernel by kernel. The re-splitting cost has a scale-invariant consequence for dense multiplication that Figure 9 makes concrete. For a square output the continuous idealisation λ = E/256 grows with the edge exactly as the work does (n̄ = E), so the per-FLOP deconstruction load is constant and the emulated rate settles at ≈235 TFLOPS on Rubin, 0.50 of set S’s 473-TFLOPS roof, for cluster-aligned large squares. Under the realizable integer schedule λ = ⌈E/256⌉ this is a sawtooth: its aligned upper edge is the 0.50 plateau, and its ragged worst case (a +1 one-cluster sliver) is the 0.36 rigid-schedule minimum reported below. (The floor is the deconstruction re-split load; it binds on the INT8/TMAC deconstruction-MAC pipe for the winning E- and S-routes, not on the SIMT residual.) An output that fits one cluster (≤ 2562 ) keeps λ = 1 and follows the n̄-crossover of §3; between the two, a narrow band is served by materialising the smaller operand in L2. Strikingly, the ≈235-TFLOPS largeDGEMM floor lands on NVIDIA’s preliminary ∼200-TFLOPS “Emulated DGEMM” figure from an independent direction to the n̄=500–900 crossover of contribution (1): a DGEMM benchmark runs large square shapes, exactly the regime the re-split floor governs. Table 12 carries this through the traced workloads and the streaming kernels. The picture is benign wherever real shapes live: block-Krylov (LOBPCG) GEMMs are tall/skinny, so λeff stays O(1) (1.25/1.50 at b=64/128 under square dispatch, 1.06/1.25 shape-matched) and the hit is mild (0.89–0.94×); multifrontal LU and the rank-64 trailing updates of blocked dense LU are memory-bound (small k), sitting below the floor and untouched (1.00×); dense QR has no large-square/large-k product (0.95×); only CCSD, whose contractions are genuinely large-k (k∼1953) and compute-bound, takes a real 0.82×. Sparse and streaming kernels (SpMV, SpMM, stencils) are exempt outright: each value is converted once and consumed immediately on-chip, with no reuse to re-split. The 45
(a) Rubin: dense DGEMM
500
single-cluster: roof-bound multi-cluster: deconstruction-λ bound
300
200
125 100
materialisation is feed-starved: both ceilings below on-the-fly
100
FP8 roof (S) 135 TF on-the-fly emulation (this design) deconstruction-λ floor 135 TF L2-materialise ceiling 68 TF HBM-materialise ceiling 17 TF
150
emulated FP64 rate (TF)
emulated FP64 rate (TF)
175
FP8 roof (S) 473 TF on-the-fly emulation (this design) deconstruction-λ floor 235 TF L2-materialise ceiling 188 TF HBM-materialise ceiling 47 TF
600
400
(b) GB300: dense DGEMM
75 50 25
0
0 64
128
256
512
1k
2k
4k
64
square output edge (elements), k = 4096
128
256
512
1k
2k
4k
square output edge (elements), k = 4096
Figure 9: Emulated dense DGEMM under on-the-fly generation. Outputs up to one cluster’s reach (2562 ) keep λ = 1; larger square outputs are pinned at the deconstruction-λ floor (≈235 TFLOPS Rubin, 0.50 of set S’s 473-TFLOP roof; the GB300 floor coincides with its lower 135-TFLOP roof, so GB300 is roof-bound throughout). Every materialised alternative (L2- or HBM-resident planes, dash-dot/dotted) is feed-starved below the on-the-fly curve, which is why the design generates on the fly. Table 12: How the deconstruction-λ limit manifests per kernel (Rubin). Only large dense square GEMMs with large k sit at the floor; the traced workloads are governed by their real shapes and are largely immune. “B/A” is the composite ratio of the honest on-the-fly schedule (B) to the materialised one (A). The tall/skinny λeff entries assume the square 4×4 cluster dispatch the engine models; a shape-matched rectangular cluster lowers them to 1.06/1.25, so that row is conservative. kernel / regime 2
Dense DGEMM ≤ 256 Dense DGEMM ≳ 5122 , lg. k LOBPCG (b=64/128) Multifrontal LU CCSD contractions Dense LU trailing Dense QR Sparse SpMV/SpMM, stencil
shape
mechanism
outcome
1 cluster multi-cluster tall-skinny n̄ med. 665 n̄ 441, k∼1953 rank-64 k ∈ {32, 3232} streaming
λ=1; SMEM-fed on-chip re-split λ=E/256; feed-starved λA =1; λeff ≈1.25/1.5 (bounded) memory-bound fronts, below floor large-k, compute-bound (53% work) small k=64 ⇒ memory-bound no large-square/large-k product converted once, consumed on-chip
crossover model ≈235 TF (0.50 of 473) mild (0.89–0.94×) immune (1.00×) hit (0.82×) immune (1.00×) near-immune (0.95×) exempt
limit, in short, is a large-dense-DGEMM phenomenon, not an application-throughput one. Coverage, restated: a deconstruction-λ bound, not an L2 one. The earlier draft’s “no matrix-size regime is left uncovered” and its blanket ≥ 0.83 figure are withdrawn; so, one round later, is the intermediate “0.79, L2-bandwidth-conditional” figure of the draft that first charged plane-feeding—it charged the feed to a materialised schedule the design does not use, and (through a 1024-element cluster reach in place of the correct 256) undercounted the resplit multiplicity that actually binds. Both are replaced by the on-the-fly statement. The enumeration domain is the same boundary-oriented suite (clamped to ≤ 4096, ±1 probes around every λ boundary min(m, n)=256j and the L2-capacity widths), the coverage ratio uses one reconstruction world in numerator and denominator, and the numerator is now the best of on-the-fly, hybrid, and (correctly feed-charged) materialised routes. The operative bound on the multi-cluster shapes is the on-the-fly re-split load: the switched coverage minimum is the deconstruction-λ floor, 0.50 of the top (473-TFLOPS) roof, cluster-aligned, on Rubin (at 46
large squares such as m=n=666), and 1.00 on GB300—whose 135-TFLOP fp8 roof sits at or below its own re-split floor, so the roof binds first and no dip appears.1 Two conditions on the floor should be stated, both Part-3 measurables: it is a large-k value (the Garner residual drags it to ≈0.35 at k=256, recovering by k≳1024), and it assumes the 16-CTA cluster opt-in (an 8-CTA-portable cluster is a 4×2 grid, which by Eq. (19) lowers the reach to 170.67 and the floor to ≈157 TFLOPS; a 32-CTA cluster is 8×4, reach 341.33, floor ≈314—see Table 6). Crucially, because the planes are never L2-resident, the minimum is no longer conditioned on BL2 : it is a compute (deconstruction) bound, not a bandwidth one, and is therefore robust to the BL2 assumption the previous round leaned on. The per-record engine and enumerator (params.py: rate_rec/env_rec/s0_rec, coverage_min) apply these predicates to every trace record and emit per-record predictions (traces/sched_*.csv), from which Table 11 is computed. Persistent sparse operators. Mode P applied to a sparse value array (operand-stationary O1) changes the dominant value-plus-index stream from ≈12 B/nnz (8+4) to ≈34 B/nnz (30+4) for set S—×2.83—capping delivered throughput near 0.35 of the raw-value memory roof. Operand-stationary conversion is therefore a bandwidth-for-instructions trade, not a free amortisation: strictly worse than fused conversion on GB300 (where the fused SIMT residual fits the keep-up budget), comparable on Rubin (0.35 against the fused 0.36–0.39), and preferable mainly when the plane array is L2-resident or SIMT issue is contended.
7
Validation Plan and Falsifiable Predictions
Everything above is a model projection; its constants come from the NVIDIA memo’s counted kernel [2] (stage provenance in Appendix B) and from official specifications [18, 17]. In the two-way spirit of the TME model [11], we state what measurement would confirm or refute, and note that each outcome identifies its cause: Prediction 1 (Linearity). Below the knee, the delivered throughput of emulated DGEMM on the CRT/fp8 path is linear in n̄—not flat—with slope Pint /(cq r) ≈ 0.39 TFLOPS per unit n̄ at the reference constants. Prediction 2 (Knee position). The signature is storage-mode-dependent, and the sweep must fix and report the mode. On a convert-once sweep—single-cluster tiling, or a materialised kernel with L2-resident planes—the transition to the ceiling occurs at n̄ ≈ n∗ (cq ) of Eq. (12) and the fitted knee measures the delivered cq . On the deployable on-the-fly large-square path the curve instead knees at the cluster reach (n̄ ≈ 256) into the deconstruction-λ floor (≈0.50 of the 473TFLOPS roof, §5.2) and never reaches n∗ ; there cq is read from the rising-branch slope Pint /(cq r) of Prediction 1, not the knee. Prediction 3 (Substrate structure). Under the fourth-term reading, a convert-once size sweep shows two knees at the positions each path’s constants predict—n̄ ≈ 730 on the B200 CRT/int8 path, n̄ ≈ 1,211 on the Rubin CRT/fp8 path—whereas the deployable on-the-fly large-square sweep shows, on each path, the same rising slope up to a plateau at the deconstruction floor near the cluster reach; either signature is constant-predicted. An error-free-slicing int8 kernel on B200 (deconstruction by byte extraction; crossover below n̄ ≈ 25) shows neither at practical sizes and serves as the control arm. Flat sub-ceiling margins on all three paths, in either storage mode, refute the branch hypothesis and indict generic efficiency artifacts instead. 1
The strict enumerated minimum under rigid integer-⌈·⌉ re-splitting is 0.36 (Rubin) / 0.80 (GB300), at (513, 513)—a +1 ragged boundary where a one-element sliver is charged a full extra re-split (area-weighted λ=2.004 but ⌈513/256⌉=3). A ragged-edge-aware schedule (predicated partial clusters, standard in production GEMM) is modeled to hold 0.54 there—a Part-3 validation target, since we enumerate it but do not construct it as a minimum. We therefore headline the cluster-aligned 0.50 plateau and report 0.36 as the only rigid-schedule worst case the enumeration establishes.
47
Prediction 4 (Ladder shift, hierarchical). An Ozaki 2.5 prototype is validated in stages, each falsifiable before the next is attempted: (i) exact residue equivalence of the two-limb reduction against a trusted integer modulo, exhaustively sampled over signed 64-bit inputs per modulus; (ii) isolated conversion throughput at the logical paired K=8, N =r (or single K=8, N =2r, separate limb outputs) shape; (iii) concurrent conversion/fp8 throughput (the Pareto frontier, not two isolated peaks); (iv) end-to-end knee movement toward n̄ ≈ 480–730 on Rubin-class parts; (v) application-level benefit on traced call geometries. A negative control—injecting the same volume of dummy SIMT/int8 work—separates deconstruction service demand from generic small-GEMM inefficiency. The decisive experiment is a single DGEMM size sweep per path (square and tall-skinny shapes to separate n̄ from matrix volume), over the three arms of Prediction 3—Rubin CRT/fp8, B200 CRT/int8, and a slicing int8 control—each sweep fixing and reporting its storage mode (singlecluster/convert-once vs. multi-cluster on-the-fly), since the two expose different but constantpredicted signatures (the n∗ knee vs. the reach-256 deconstruction floor of §5.2); it decides the branch hypothesis, fits each path’s cq , and prices the recovery in one run. The conversionkernel derby of Part 1—including the DP4A and tensor-migrated variants against the memo’s eff complete the ∼7-instruction criterion [2]—and the integer-pipe census that fixes the true Pint set. These are the opening entries of the follow-up implementation and measurement work (Part 3), and are being automated on RIKEN’s Rikyu GB200 NVL4 system so that the sweep re-runs continuously as kernels evolve (the GB200 testbed calibrates the model; it does not by itself validate GB300’s reduced integer-tensor balance or Rubin’s concurrency, which require those parts). Pass/fail criteria. Each stage carries an explicit criterion: (correctness) bit-exact agreement with an arbitrary-precision reference over exhaustive 16-bit and stratified random signed 64bit inputs per modulus, including even-modulus boundaries (−m/2), ±0, subnormals, and the API’s infinity/NaN contract; (accounting) SASS instruction counts per output residue, achieved dp4a issue rate, register count, and occupancy—this directly adjudicates the L1 direct count of §5 and its ≈7.5 target; (reduction shape) useful and issued operations at the exact pairedlimb shapes, isolated and concurrent, fixing ηred ; (concurrency) the paired-rate measurements: isolated fp8, int8-MMA, and dp4a rates; then fp8+int8 and fp8+dp4a at matched occupancy, each pair reported against its serialised, ideal-overlap, and measured predictions—fixing θ and ρsimt,fp8 ; (knee) a piecewise fit with confidence intervals against the size-independent-efficiency alternative on held-out shapes, not visual inspection; (numerics) the source accuracy suite rerun per candidate set with adversarial scaling/conditioning; (applications) the to-be-released traces replayed through measured kernels, GEMM-portion and whole-run Amdahl numbers reported separately. Claim status. For the record, the paper’s claims separate as follows. Algebraically proved: the harmonic-mean crossover under the reduced model; the two-limb byte identity and its 219 accumulation bound; coprimality and exact CRT products of the released candidate sets. Script-checked (scripts in the artifact: lemma_search.py, verify_candidates.py, stail_layout.py): the exhaustive supply bounds; per-set coprimality and exact CRT products; canonical centred digit maps and E4M3 representability of all A/E/D planes and the six square S moduli; the carry-corrected layout of the six nonsquare S-tail moduli, exhaustive over all m2 ordered centred pairs per modulus, together with the canonical (321) and redundant (577) envelopes that bracket it; the signed two-limb identity (exhaustive 16-bit plus stratified 64-bit). Ledger-projected (uncompiled instruction counts): the SIMT residuals (the uniform ≈6.3 and the per-set 5.8/5.3/5.2 for S/E/D, 5.27 for A) and the dp4a route counts (10.3/7.3)—instruction-ledger projections, neither compiled nor measured. Memo-counted: the cq =16 shipping-path SASS estimate (Appendix B). Script-generated (params.py): every 48
number in Tables 3–11, as deterministic evaluation of the stated constants—generation is not independent verification of those constants. Measured (fixed software stack): the LD_PRELOAD call-shape traces, and nothing else. Modeled (conditional): all route knees, envelopes, and composite speedups. Hypothesized: the deconstruction-limited reading of the Rubin specification. To be measured: the L1 reuse schedule, ηred , the concurrency parameter θ, storage-mode boundaries, end-to-end knees, per-set numerical behaviour, and application replay. The artifact (params.py, lemma_search.py, verify_candidates.py, the traces, the per-record predictions traces/sched_*.csv, and the table manifests), with a manifest mapping every table and figure cell to its generating command, is supplied with this submission as a reproduction artifact—an archival Zenodo DOI and commit hash will be added in the first arXiv revision and cited in the camera-ready. In the revised artifact the per-record CSVs (traces/sched_*.csv) emit the selected route, storage mode, reconstruction host, Q0 , Qplane , Qred , and the selected rate as audit fields, and a single-command driver regenerates every table, figure, and report from the raw traces with recorded tool versions. Trace quantiles and composites are FLOP-weighted. params.py is the single parameter source.
8
Conclusion
Ozaki II made fp64 dense matrix multiplication a schedule of fp8 tensor-core operations; Rubin made the result a product specification. This paper has argued that part of the distance between that preliminary specification (∼ 200 TFLOPS) and the raw arithmetic roof (≈ 473 TFLOPS) may have a specific, measurable, and largely removable cause: the deconstruction of streamed fp64 operands into CRT residues, a cost the int8 substrate could largely avoid—by errorfree slicing, or at least Karatsuba-free—and that the fp8 substrate, mandatory once int8 tensor throughput is deprecated, pays in full per modulus on every element. The size crossover n∗ = cq rPfp8 /(αPint ) makes the hypothesis quantitative: at the memo-derived constants the published figure is consistent with the deconstruction-limited branch at model-inferred sizes n̄ ≈ 500–900, and a single size sweep—with a slicing int8 kernel as control—decides the question while measuring cq . Honesty about the ceiling must close the paper as it opened it. On today’s silicon the deliverable for large dense DGEMM —the regime HPL runs in—is not the roof but the deconstruction-λ floor, ≈235 TFLOPS on Rubin (0.50 of the common 473-TFLOPS roof, or equivalently 0.54 of set E’s own 438; ≈8× native, ≈1.9 EFLOPS fp64 sustained on a 10,000-GPU cluster) on the codesigned set E at NB ≳ 1024, or ≈182 TFLOPS (0.38, ≈1.4–1.5 EFLOPS) on the published set S, which is the number backed by the round-to-nearest theorem (Table 8). The rising envelopes of Fig. 3 are achievable only up to one cluster’s reach; larger outputs plateau there. That floor is R Pint /(cq r), liftable to the roof only by the co-design of §5.2—minimally an in-flight-convert copy-engine datapath (Plan A), with a higher-bandwidth L2 under softwareblocked materialisation as the datapath-independent Plan B. The equally important positive is that this floor is a large-square phenomenon: the tall/skinny and small-batch matrices that dominate real solvers—block-Krylov, batched GEMMs, panel factorisations—re-split their big operand only once (λA =1; size-weighted λeff ≈1.25–1.5, bounded, not the ∝edge of a large square), so they stay near the crossover (a mild 0.89–0.94×, not the 0.50 floor) and Ozaki 2.5 is already useful on Rubin today, worth ≈1.6–1.9× (≈2× on the block-Krylov rows) over simple deconstruction with no hardware change (Table 11), and GB300 is roof-bound outright. Only the headline large-square-DGEMM number is half the roof until the hardware changes, and we state it that way—clipping-limited, not fundamental—rather than leaving the roof as the last impression. The constructive contributions are the Ozaki 2.5 method and the modulus/encoding codesign it led to. The method adds no new mathematics: convert once, reduce exactly on the integer tensor side (two limbs, because the actual moduli exceed a byte), keep only the irreducible 49
nonlinearities on SIMT, and hide the rest behind the MMA stream. Under stated assumptions the reduced model projects the Rubin knee moving from ≈ 1,211 to ≈ 480–730 and a conditional 1.8–2.2× at the announced operating region (364–438 TFLOPS at n̄=512, route-dependent—a convert-once envelope, delivered to within a mild clip (λeff ≈1.25–1.5) for tall/skinny but clipped to the ≈235-TF floor for a large square); the measured result, whatever it is, will replace these projections. The codesign study of §5.1 then shows the modulus set itself is a first-class, runtime-switchable performance parameter: an all-byte system legalises the one-pass reduction and is projected faster below n̄ ≈ 540 despite a 20% lower roof, and the hybrid concedes only 7.5% asymptotically while leading the ≈410–620 band—so a library can dispatch by shape and lose almost nothing anywhere. None of this is Rubin-specific: instantiated at Blackwell Ultra (GB300) rates the same machinery gives 135/109/125-TFLOPS roofs but a different winner—the ∼166-TOPS residual integer-tensor rate throttles the tensor-migrated reductions, compressing the small-size band, which the all-byte set leads (its dp4a realisation within ≈7%, needing no integer tensor at all), with the two-limb route taking the roof from n̄ ≈ 290—so even the dispatch itself is platform-dependent, which is precisely what a parameterised model is for. And because GB300’s native fp64 pipe delivers only ≈1.4 TFLOPS, every modeled route exceeds the native reference throughout the evaluated range n̄ ≥ 32: Ozaki 2.5 is a proposition for the Blackwell generation already shipping, not one that waits for Rubin. The measured call-shape traces of §6.1 ground this where it matters: block-Krylov and QR-panel work is block-width-bounded in real runs, and multifrontal fronts sit in the crossover band itself. For workloads whose geometry keeps one output dimension at a block width (tall-skinny and block-vector GEMMs), these gains, if confirmed, apply throughout that band; for what software cannot reach, Part 1’s hardware options remain the outlook—Option-C-class conversion would retire the term at every size—but on current GPUs the switched software dispatch already does well across realistic matrix sizes. Whether the branch hypothesis survives measurement or falls to it, the outcome identifies the responsible term—which is exactly what a performance model is for.
Acknowledgments The author is indebted to Harun Bayraktar, John Gunnels, and Peter Caday of NVIDIA, whose technical note [2] (cited with permission) identified the deconstruction term on which this paper rests, and to Dan Ernst and Matthew Martineau for discussions that helped shape the experimental programme; the follow-up implementation and measurement work that will test the predictions made here is in preparation (Part 3). The author also thanks the RIKEN R-CCS teams standing up the Rikyu GB200 NVL4 measurement harness. This work was undertaken as part of the FugakuNEXT project and related R-CCS initiatives on AI for Science. Disclosure of AI-assisted writing. This manuscript was prepared with assistance from large language models: Anthropic’s Claude (Opus 4.8 and Fable 5, with the final pre-submission pass of Draft 27 by Fable 5.1) served as authoring assistants for drafting, the derivation checks, figure generation, and LATEX mechanics, and OpenAI’s GPT-5.6 (Codex) provided critical technical review across drafts, under the author’s direction. All scientific arguments, performance projections, and conclusions were directed, reviewed, and validated by the author, who takes full responsibility for the content, including any errors of fact or judgment.
A
Proof Obligations per Candidate Modulus Set
Changing the modulus set leaves the Ozaki-II reconstruction framework intact but does not inherit the source paper’s numerical argument automatically. Table 13 fixes the numericalcontract variants by name; every accuracy statement in this paper is scoped to one of them.
50
Table 13: Named numerical-contract variants and their status. variant
conversion rule
status
Ozaki-II-RN
round-to-nearest
Ozaki-II-TZ
implemented truncation, plus the stated endpoint rule
Codesigned A/E/D
per-set digit and endpoint rules
theoretical: the cited error theorem [24] holds under its own hypotheses proposed/unvalidated: no automatic inheritance from RN; a proof or a measured error distribution is required obligation: performance case made here; accuracy is an explicit per-set obligation, items (i)–(vi) below
For each released candidate (A, E, D) the following obligations are stated here and scriptchecked (verify_candidates.py, in the artifact with its output manifest); SASS-level and accuracy-suite confirmation belongs to §7. The checker reports pass/skip/fail per obligation and never returns an all-pass verdict while any obligation is skipped. (Candidate D, the one set not detailed in §5.1, is the 7-bit system {128, 127, 125, 123, 121, 119, 113, 109, 107, 103, 101, 97, 89, 83, 79, 73, 71} (r=17, α′ =52, one-pass reduction with signed-byte constants); the list ships in verify_candidates.py and the artifact manifest.) (i) Range: pairwise coprimality and the exact CRT product (log2 P = 117.8 for A, 113.3 for E, 113.5 for D; the published set’s 111.8; log2 P —total CRT product bits—is the single range metric used throughout this paper, the usable symmetric range being P/2, exactly one bit less)—each meets or exceeds the source requirement. (ii) Digit map: residues carried centred in [−m/2, m/2), even-modulus endpoint included at −m/2 (for m=256, the value −128); under the canonical tie rule d0 ∈ [−8, 7], d1 = (x−d0 )/16 the decomposition is unique on that interval (without a tie rule balanced digits are redundant at ties, e.g., 8 = 0 · 16 + 8 = 1 · 16 − 8). (iii) Representability: every stored digit and every Karatsuba sum actually used is e4m3-exact (checked exhaustively per modulus). (iv) Accumulation: all reduction dot products bounded by 8 · 255 · 255 < 219 , exact in int32; the fp32 MMA accumulation bound and k-blocking rule of the source apply unchanged per set. (v) Signed equivalence: the limb formation with the (−264 ) mod mi correction equals a trusted signed integer modulo on exhaustive 16-bit and stratified 64-bit samples. (vi) Contract: scaling, truncation, ESC/ADP escalation, and reconstruction conditions of the source analysis are reverified per set as part of the validation suite—until then the accuracy contract for A/E/D is a stated obligation, not a result.
B
Provenance of the Memo-Counted Constants
The model’s anchor constants—cq = 16 and the Pint ≈ 75 TOPS normalisation—derive from a private NVIDIA technical note [2], cited with permission. To make the anchor auditable we reproduce the non-confidential stage structure; the note’s authors have been asked to confirm this summary for publication. The confirmation request also asks which emulation realisation the counted SASS belongs to—the released Scheme-I slicing path or a CRT (Scheme-II-class) path— since the stage structure is CRT-flavoured while the released cuBLAS path is publicly described as Scheme-I [10, 16]; §7 carries the reading as conditional pending that answer. Stages, per streamed element and per modulus on the shipping conversion path (SASS-counted estimates, not throughput measurements; architecture and toolchain preliminary): scale and truncate-tointeger, ∼4 per element (amortised 4/r per modulus); per-modulus byte-plane reduction, ∼8–10 (constant division compiling to multiply-high sequences); final reduction to the residue interval, ∼2; Karatsuba split into three fp8 pieces, ∼3; totalling cq ≈ 16 under the stated amortisation. One dp4a counts as one issued SIMT instruction throughout (ISIMT = 75 · 1012 instructions/s in the memo’s normalisation); int8 tensor capacity is counted separately in operations/s (PI8TC , two operations per MAC). Nothing else from the note is used. 51
Reconciliation with the spectral companion. The FFT companion [12] normalises the same pipe at 41.7 T inst/s, the published B300 fp32 TFLOP/s divided by two on the assumption that fp32 and integer issue share lanes on Blackwell. This paper keeps the memo’s 75 as its reference constant because it is the normalisation the memo’s counts were made against; the 1.8× gap between the two is exactly what the Part-3 integer-pipe census measures. Every ISIMT bounded quantity here—the SIMT keep-up budget (6.25 → ≈3.5 instructions per modulus on GB300, 2.3 → ≈1.3 on Rubin at 41.7), the SIMT-residual knees, and the Garner keep-up thresholds—tightens monotonically under the spec-derived value, so the 41.7 case is a uniform downside to the software-only picture and a uniform strengthening of the co-design case. The tensor-side quantities (473 TFLOPS; the 0.497 ratio of the deconstruction-λ term to the roof in Eq. (15); the int8-tensor knees) do not depend on it, but the floor itself does, because it is a minimum over terms: re-running the engine at 41.7 T inst/s, the SIMT-residual term becomes binding on Rubin and the large-square floor falls from 235 to ≈148 TFLOPS (0.31 of 473; routes E and S at 148 and 147, A at 128), and on GB300 the dp4a route L1d that supplies the 135 TFLOPS floor drops to 84, leaving route E at its own 125 TFLOPS roof as the floor. Route E leads on both parts under either normalisation; what the census settles is whether the floor is one half of the roof or one third of it, which makes Pint the single largest sensitivity this paper carries. The 75-normalised figures are the ones this paper quotes; the 41.7 case is the stated sensitivity, and both floor values are engine-checked (paper_numbers.py, check_pint_sensitivity).
References [1] Patrick R. Amestoy, Iain S. Duff, Jean-Yves L’Excellent, and Jacko Koster. A fully asynchronous multifrontal solver using distributed dynamic scheduling. SIAM Journal on Matrix Analysis and Applications, 23(1):15–41, 2001. [2] Harun Bayraktar, John Gunnels, and Peter Caday. The missing fourth term for the emulation TME model: The residue deconstruction cost that temu omits, in the nomenclature of arXiv:2606.06510. Technical report, NVIDIA Corporation, July 2026. Technical note, v2, July 10, 2026; shared with the author in the course of technical review. Cited with permission. [3] Centre for High Performance Computing (CHPC), South Africa. HPL-CUDA: running and tuning GPU-accelerated HPL. CHPC ACE Lab wiki, https://wiki.chpc.ac.za/acelab: hpl_cuda, 2024. Recommends NB = 892–1024 for GPU runs, against 192–256 for x86-only runs. [4] Harvey L. Garner. The residue number system. IRE Transactions on Electronic Computers, EC-8(2):140–147, 1959. [5] Kempner Institute, Harvard University. NVIDIA HPC-Benchmarks: HPL run recipes and results. https://github.com/KempnerInstitute/nvidia-hpc-benchmarks, 2025. Reported single- and multi-GPU HPL runs use NB = 1024. [6] Donald E. Knuth. The Art of Computer Programming, Volume 2: Seminumerical Algorithms. Addison-Wesley, 3rd edition, 1997. [7] Andrew V. Knyazev. Toward the optimal preconditioned eigensolver: Locally optimal block preconditioned conjugate gradient method. SIAM Journal on Scientific Computing, 23(2):517–541, 2001.
52
[8] Xiaoye S. Li and James W. Demmel. SuperLU_DIST: A scalable distributed-memory sparse direct solver for unsymmetric linear systems. ACM Transactions on Mathematical Software, 29(2):110–140, 2003. [9] Glenn K. Lockwood. NVIDIA Rubin: Architecture notes and performance specifications. Glenn’s Digital Garden, 2026. https://www.glennklockwood.com/garden/processors/ R200, accessed May 2026. [10] D. Lu, A. Maeder, M. Luisier, and A. N. Ziogas. EmuGEMM: Fused tensor core kernels for precision emulation in matrix multiplication. arXiv:2606.25453, 2026. [11] Satoshi Matsuoka. FP8 is all you need (Part 1): Debunking hardware FP64 as the HPC holy grail — a tensor–memory equilibrium model and implementation strategy for Ozaki Scheme II on memory-bound workloads in the post-FP64 era. arXiv:2606.06510 [cs.DC], 2026. Revision 37 of September 3, 2026. [12] Satoshi Matsuoka. FP8 is all you need (Part 2): Full-FP64 3-d FFT on FP8-generation tensor cores — the integer-epilogue wall and the minimal hardware that would remove it. arXiv:2606.23698 [cs.MS], 2026. Revision 10 of September 3, 2026. [13] Daichi Mukunoki. DGEMM without FP64 arithmetic: Using FP64 emulation and FP8 tensor cores with Ozaki scheme, 2025. [14] NVIDIA Corporation. Unlocking tensor core performance with floating-point emulation in cuBLAS. NVIDIA Developer Blog, 2025. [15] NVIDIA Corporation. CUDA C++ Programming Guide, 2026. Accessed July 2026. On thread-block-cluster launch (cudaLaunchAttributeClusterDimension with cudaLaunchKernelEx): “the corresponding dimensions of the grid (x, y, and z) must be divisible by the respective dimensions of the specified cluster dimension.” A cluster of 8 thread blocks is the portable maximum; larger clusters are architecture-specific and require the non-portable-size opt-in. [16] NVIDIA Corporation. cuEST: Emulation and mixed precision. https://docs.nvidia. com/cuda/cuest/emulation.html, 2026. Accessed July 2026; exposes slice-count (SchemeI) and modulus-count (Scheme-II-class) controls, selected by compute capability. [17] NVIDIA Corporation. Inside the NVIDIA Vera Rubin platform: Six new chips, one AI supercomputer. NVIDIA Developer Blog, 2026. Lists “Emulated DGEMM” as an official column in Rubin specifications. [18] NVIDIA Corporation. NVIDIA DGX Rubin NVL8: Supercharged infrastructure for agentic AI (product specification). https://www.nvidia.com/en-us/data-center/ dgx-rubin-nvl8/, June 2026. Specification page published June 30, 2026 (product announced at CES 2026); lists 140 PFLOPS dense FP8/FP6 training per eight-GPU system, i.e. 17.5 PFLOPS dense FP8 per Rubin GPU. Values marked preliminary and subject to change. [19] NVIDIA Corporation. NVIDIA GB300 NVL72: Specifications. https://www.nvidia. com/en-us/data-center/gb300-nvl72/, 2026. Accessed 21 July 2026; rack totals with sparsity, unpacked to dense per-GPU rates in the text. [20] NVIDIA Corporation. NVIDIA HGX platform: Specifications. https://www.nvidia. com/en-us/data-center/hgx/, 2026. Accessed 21 July 2026; lists per-Rubin-GPU dense FP8 17.5 PFLOPS, dense INT8 250 TOPS, FP64 33 TFLOPS, and 200 TFLOPS FP64 DGEMM via tensor-core emulation. 53
[21] NVIDIA Corporation. Parallel Thread Execution ISA, 2026. Accessed July 2026 (current version 9.3). The cluster-scope primitives used here— cp.async.bulk.shared::cluster.shared::cta.mbarrier::complete_tx::bytes, mapa, and .cluster-scope mbarrier arrive/wait—were all introduced in PTX ISA 8.0 and require target sm_90 or later. [22] Katsuhisa Ozaki, Takeshi Ogita, Shin’ichi Oishi, and Siegfried M. Rump. Error-free transformations of matrix multiplication by using fast routines of matrix multiplication and its applications. Numerical Algorithms, 59(1):95–118, 2012. [23] Katsuhisa Ozaki, Yuki Uchino, and Toshiyuki Imamura. Ozaki Scheme II: A GEMMoriented emulation of floating-point matrix multiplication using an integer modular technique. arXiv:2504.08009, 2025. [24] Katsuhisa Ozaki, Yuki Uchino, and Toshiyuki Imamura. Error analysis of matrix multiplication emulation using Ozaki-II scheme. arXiv:2602.02549, 2026. [25] A. Schwarz, A. Anders, C. Brower, H. Bayraktar, J. Gunnels, K. Clark, R. G. Xu, S. Rodriguez, S. Cayrols, P. Tabaszewski, et al. Guaranteed DGEMM accuracy while using reduced precision tensor cores through extensions of the Ozaki scheme. In Proceedings of the Supercomputing Asia and International Conference on High Performance Computing in Asia Pacific Region. ACM, 2025. Also available as arXiv:2511.13778. [26] Yuki Uchino. GEMMul8: GEMM emulation using INT8/FP8 matrix engines based on the Ozaki Scheme II. GitHub repository, RIKEN-RCCS, 2025. [27] Yuki Uchino, Katsuhisa Ozaki, and Toshiyuki Imamura. Double-precision matrix multiplication emulation via Ozaki-II scheme with FP8 quantization. arXiv:2603.10634, 2026.
54