EmuGEMM-I: Fused Ozaki Scheme I on NVIDIA GPUs

• EmuGEMM-I is a fused Tensor Core implementation of Ozaki Scheme I that emulates high-precision GEMM by decomposing inputs into INT8 slices and reconstructing outputs on chip.
• It integrates all slice-pair products and the shift-reduce reconstruction into a single persistent kernel, dramatically cutting global memory traffic.
• Optimized for NVIDIA Hopper and Blackwell GPUs, it achieves up to 89% of INT8 peak throughput and significant speedups over native FP32/TF32 operations.

Searching arXiv for the relevant papers on EmuGEMM and related Ozaki-scheme GEMM emulation. arXiv search query: "EmuGEMM fused tensor core kernels precision emulation matrix multiplication" EmuGEMM-I is a fused Tensor Core implementation of Ozaki Scheme I for general matrix multiplication on NVIDIA Hopper and Blackwell GPUs. It reconstructs higher-precision GEMM from exact INT8 slice products with INT32 accumulation, but differs from earlier Ozaki implementations in a decisive systems-level respect: all required slice-pair products for an output tile, the triangular INT32 accumulators, and the final shift-reduce reconstruction are kept on chip inside a single persistent kernel, rather than being repeatedly materialized in global memory. The result is a Scheme-I realization that moves the bottleneck away from intermediate data movement and back toward Tensor Core execution (Lu et al., 24 Jun 2026).

1. Definition, scope, and position in precision emulation

EmuGEMM-I is the Scheme-I component of the broader EmuGEMM framework. Its purpose is to recover higher-precision GEMM from much faster low-precision INT8 Tensor Core operations while avoiding the memory bottlenecks that made earlier Ozaki implementations underperform their theoretical potential. In the formulation used here, the input matrices are first decomposed into multiple INT8 mantissa “slices”; EmuGEMM-I then computes all required slice products in one fused kernel, reconstructs the final tile on chip, and writes only the final FP32 or FP64 output tile to global memory (Lu et al., 24 Jun 2026).

The method is motivated by the widening precision–throughput gap on modern GPUs. Hopper and Blackwell devote substantially more hardware to low-precision matrix engines than to FP64 or FP32 arithmetic, so direct high-precision GEMM increasingly leaves available compute underused. EmuGEMM-I addresses that asymmetry with a structured emulation strategy rather than a new arithmetic format. Its defining contribution is therefore implementation-centric: it preserves the arithmetic structure of Ozaki Scheme I while changing the dataflow so that the expensive intermediate INT32 state never leaves the chip during the fused compute phase.

The intended operating region is the low-to-mid precision range. The evaluation positions EmuGEMM-I SGEMM at slice counts $p=1$ to $p=4$, and EmuGEMM-I DGEMM at $p=5$ to $p=8$. This already indicates a central tradeoff: Scheme I scales precision by increasing the number of slices, but its compute cost grows quadratically with $p$. A plausible implication is that EmuGEMM-I is strongest where a modest slice count is sufficient and where eliminating global-memory round trips materially changes the roofline regime.

2. Numerical foundation: Ozaki Scheme I and shift-reduce reconstruction

EmuGEMM-I implements Ozaki Scheme I exactly as a mantissa-splitting emulation method. Each operand is decomposed into $p$ INT8 slices, each slice representing a $\beta$-bit mantissa segment with $\beta \le 8$, and the matrices are scaled so that rows of $A$ and columns of $B$ have aligned exponents. The paper writes the decomposition as

$p=4$0

Here $p=4$1 and $p=4$2 are the INT8 slice matrices, and $p=4$3 and $p=4$4 are power-of-two scaling vectors. The decomposition is described as “error-free up to a residual that diminishes with increasing $p=4$5,” and each additional slice contributes about $p=4$6 bits of precision (Lu et al., 24 Jun 2026).

Because slice $p=4$7 carries weight $p=4$8 and slice $p=4$9 carries weight $p=5$0, the slice product $p=5$1 has weight $p=5$2. Scheme I groups all terms with equal total shift $p=5$3 into triangular accumulators

$p=5$4

The final result is then reconstructed as

$p=5$5

The paper calls this weighted summation a shift-reduce. Since the weights are exact powers of two, the reconstruction introduces no extra rounding beyond the decomposition residual.

The compute phase therefore requires

$p=5$6

independent INT8 GEMMs in total. That quadratic count is the central algorithmic liability of Scheme I, but it is also what EmuGEMM-I exploits structurally: instead of launching those GEMMs separately, it executes the triangular schedule inside one persistent kernel.

The implementation also adopts the improved Scheme-I slicing approach used by cuBLAS: a signed leading slice (SINT8) with unsigned remaining slices (UINT8), gaining one bit of precision per slice over the original signed-everywhere formulation. This detail matters because EmuGEMM-I is not proposing a new numerical decomposition; it is realizing an already favorable decomposition more efficiently.

3. Fused kernel structure and interleaved data layout

The principal obstacle in prior Scheme-I implementations was not the INT8 arithmetic itself but the repeated materialization of intermediate INT32 outputs. In a naive realization, each of the $p=5$7 slice-pair products is launched as a separate GEMM kernel, those kernels write INT32 outputs to global memory, and a later reconstruction kernel reads them back. The paper models the global-memory traffic of that approach as

$p=5$8

where $p=5$9 is the output element size: 4 for FP32 and 8 for FP64. The dominant inefficiency is the $p=8$0 INT32 traffic term. In the fused design, the traffic becomes

$p=8$1

This removes the entire INT32 intermediate term. The paper states that the arithmetic intensity rises by a factor of $p=8$2; for FP64 emulation with $p=8$3, that is a $p=8$4 increase (Lu et al., 24 Jun 2026).

To make that fusion practical, EmuGEMM-I uses an interleaved layout along the contraction dimension $p=8$5. For operand $p=8$6, the decomposition kernel writes

$p=8$7

Thus, for $p=8$8, the stored order is

$p=8$9

with each block $p$0 columns wide. $p$1 is interleaved analogously along rows. This layout ensures that all $p$2 slices for one $p$3-chunk can be fetched by one TMA descriptor per operand, that every slice begins at an MMA-aligned address in shared memory, and that per-slice offsets are compile-time constants.

Within the kernel, the dataflow is persistent and tiled. For each output tile, EmuGEMM-I initializes $p$4 on-chip accumulators, iterates over the $p$5-dimension in interleaved blocks, issues the triangular set of MMAs for slice pairs $p$6, and keeps the $p$7 INT32 accumulators resident on chip for the full $p$8-loop. Only after the loop ends does the epilogue convert the accumulators to floating point, perform the shift-reduce, apply the row and column scales, and store the final FP32 or FP64 tile.

A crucial boundary is explicit: the initial decomposition / packing / interleaving preprocessing is not fused. EmuGEMM-I is therefore not a one-kernel path from raw floating-point inputs to final output. What is fused is the expensive middle of the pipeline: all low-precision block multiplications, all intermediate INT32 accumulation, and the reconstruction epilogue.

4. Hopper and Blackwell mapping, on-chip storage, and saturation model

EmuGEMM-I is specialized for NVIDIA Hopper and Blackwell integer Tensor Cores. On Hopper, it uses wgmma.mma_async, which is warpgroup-level, with $p$9, $p$0, and configurable $p$1. On Blackwell, it uses tcgen05.mma, with $p$2, $p$3, and configurable $p$4. In both cases, operands are staged in shared memory, but the accumulator placement differs sharply: Hopper keeps accumulators in the register file, whereas Blackwell keeps them in TMEM, a dedicated tensor memory separate from the register file (Lu et al., 24 Jun 2026).

That architectural distinction directly affects scalability with slice count $p$5. On Hopper, both thread-local state and triangular accumulators compete for the same register-file capacity, so increasing $p$6 intensifies RF pressure. On Blackwell, accumulators move to TMEM, leaving the register file primarily for addresses, epilogue state, and temporaries. This allows larger effective on-chip working sets at the same slice count and helps explain why the low-slice performance advantage is larger on B200 than on GH200.

The paper makes the on-chip scaling law explicit. Relative to a standard INT8 GEMM, the fused operand buffer scales as

$p$7

and the fused accumulator footprint scales as

$p$8

Given maximum accumulator capacity $p$9, the largest $\beta$0-stacking factor is

$\beta$1

Ordinarily each $\beta$2-step issues $\beta$3 MMA instructions. Because Scheme I computes $\beta$4 products per step, EmuGEMM-I has effective pipeline depth

$\beta$5

The stated consequence is that even if increasing $\beta$6 reduces $\beta$7, the triangular factor $\beta$8 grows quadratically, so $\beta$9 remains above the hardware saturation threshold of 16 MMA instructions per $\beta \le 8$0-step. This is the core saturation argument for fused Scheme I.

The kernel is also overlapped. Because it is persistent, TMA can prefetch the next output tile’s operands while the current tile’s shift-reduce epilogue executes. On Hopper the accumulators remain register-resident through the $\beta \le 8$1-loop; on Blackwell they reside in TMEM and are drained to registers during epilogue in pipelined fashion.

5. Performance envelope and accuracy interpretation

The headline throughput results are stated at two levels. In the end-to-end evaluation, EmuGEMM-I sustains up to 1,639 Top/s on Hopper (83% of INT8 peak) and 3,654 Top/s on Blackwell (81%). In the kernel-only roofline study on GH200, excluding slicing and reconstruction preprocessing, it reaches 1,752 Top/s, or 89% of INT8 peak; this exceeds cuBLAS native INT8 GEMM at 1,379 Top/s (70%) by 27% (Lu et al., 24 Jun 2026).

For large matrices, EmuGEMM-I surpasses cuBLAS TF32 throughput while maintaining similar effective precision. At $\beta \le 8$2 and $\beta \le 8$3, GH200 reaches 516 Tflop/s versus 365 Tflop/s for cuBLAS TF32, and B200 reaches 1,203 Tflop/s versus 710 Tflop/s, corresponding to $\beta \le 8$4 and $\beta \le 8$5 speedups, respectively. The paper describes this operating point as having “similar precision to TF32.”

In the low-precision or SGEMM-oriented regime, the gains over prior Scheme-I implementations are also explicit. At $\beta \le 8$6, $\beta \le 8$7, on GH200, EmuGEMM-I = 138 Tflop/s, cuBLAS Scheme-I emulation = 71 Tflop/s, and cuBLAS native FP32 = 50 Tflop/s. This corresponds to $\beta \le 8$8 over cuBLAS Scheme-I emulation and $\beta \le 8$9 over native FP32. For complex GEMM with Scheme I, at $A$0, $A$1, on GH200, EmuGEMM-I CGEMM = 129 Tflop/s and cuBLAS native FP32 CGEMM = 54 Tflop/s, for a $A$2 speedup.

The high-slice regime is less favorable. At $A$3, $A$4, on GH200, EmuGEMM-I = 44 Tflop/s, cuBLAS Scheme-I emulation = 32 Tflop/s, and cuBLAS native FP64 = 58 Tflop/s. EmuGEMM-I still exceeds prior Scheme-I emulation, but it does not beat native FP64 on Hopper at the highest slice counts. That is the paper’s main empirical reason for treating Scheme II as necessary for the highest-precision region.

Accuracy is measured as effective bits of precision, defined as the absolute value of the base-2 exponent of the relative error with respect to a reference product. “Comparable accuracy” therefore does not mean bitwise equivalence to TF32, FP32, or FP64. It means similar effective precision on the tested datasets. The paper also stresses that Scheme I itself is not exact overall: the exact phase is the INT8$A$5INT8-to-INT32 arithmetic, while the decomposition remains “error-free up to a residual that diminishes with increasing $A$6.”

The evaluation scope is correspondingly specific. It covers square matrices with sizes 2048–16384. For small problems, preprocessing overhead dominates, and Ozaki-scheme implementations underperform high-precision Tensor Cores. EmuGEMM-I improves substantially over earlier Ozaki implementations in that regime, but it is still not presented as a small-matrix method.

6. Relation to Scheme II, strengths, and explicit limitations

EmuGEMM-I is best understood as one half of a two-scheme design space. The paper contrasts Scheme I and Scheme II as follows (Lu et al., 24 Jun 2026):

Aspect	EmuGEMM-I / Scheme I	Scheme II
Decomposition	Floating-point mantissa splitting	Modular reduction
INT8 GEMMs	$A$7	$A$8
Scaling	Quadratic	Linear
Reconstruction	Shift-reduce	CRT
Precision gain	$A$9 bits	$B$0 bits

This distinction places EmuGEMM-I within a broader Ozaki lineage. Earlier work on Ozaki Scheme II introduced the modular-reduction and Chinese Remainder Theorem formulation precisely to avoid the triangular cross-product growth of Scheme I, replacing it with one GEMM per modulus and a CRT reconstruction stage (Ozaki et al., 10 Apr 2025). Later large-scale INT8-matrix-engine studies report that Scheme II on GH200 can achieve about 1.4× speedup over native DGEMM and 3.0× over native SGEMM for sufficiently large problems, together with substantial power-efficiency gains, but those results concern the CRT-based line rather than EmuGEMM-I itself (Uchino et al., 6 Aug 2025).

Within that division of labor, EmuGEMM-I’s strengths are clear. It is highly effective in the low-to-mid precision range, where a relatively small $B$1 already delivers useful precision and where fusing the triangular schedule yields very high Tensor Core utilization. Its core systems insight is that Scheme I was not primarily limited by INT8 arithmetic throughput; it was limited by intermediate INT32 traffic. By removing that traffic, EmuGEMM-I converts a memory-bound emulation method into a compute-efficient one.

Its limitations are equally explicit. Scheme I requires $B$2 GEMMs, so compute cost grows quadratically with slice count. Fusion also multiplies operand-buffer and accumulator requirements by $B$3, which can force smaller tiles or shallower pipelines, especially on Hopper where accumulators and thread-local state share the register file. The implementation does not include cuBLAS-style automatic dynamic precision (ADP), so slice count must effectively be chosen manually. The current implementation is Hopper- and Blackwell-only, the evaluation does not report non-square or irregular cases, and the authors do not provide a new full numerical analysis of Scheme I, instead relying on prior Ozaki literature for the arithmetic structure. The optimization target is data movement rather than the mathematics of the decomposition itself.

EmuGEMM-I is therefore best characterized as a hardware-aware realization of a pre-existing emulation method. It does not alter the numerical identity of Scheme I; it reorganizes the kernel so that the triangular INT32 state stays on chip, the interleaved slices arrive in MMA-aligned form, and the final shift-reduce occurs before write-back. In that form, Scheme I becomes competitive with, and in specific large-matrix regimes faster than, native cuBLAS TF32 and earlier Scheme-I emulation, while remaining less attractive than Scheme II near the full-FP64 end of the precision spectrum.