EmuGEMM: Fused Tensor Core Kernels for Precision Emulation in Matrix Multiplication

Abstract: Modern GPUs devote an increasing silicon budget to low-precision matrix-multiplication units, widening the precision-throughput gap for scientific computing workloads. Ozaki Schemes I and II offer an alternative by reconstructing high-precision general matrix multiplication (GEMM) from low-precision operations, yet existing implementations leave substantial performance untapped. In particular, intermediate results are repeatedly materialized in global memory, making data movement the dominant bottleneck. We present EmuGEMM, fused integer Tensor Core kernels for NVIDIA Hopper and Blackwell GPUs that eliminate redundant memory round-trips in both Ozaki schemes. Using Scheme I, EmuGEMM sustains up to 1,639 Top/s on Hopper (83% of INT8 peak) and 3,654 Top/s on Blackwell (81%). For large matrices, EmuGEMM surpasses cuBLAS TF32 throughput by up to 1.4x on Hopper and 1.7x on Blackwell, at comparable accuracy. Using Scheme II, EmuGEMM extends to complex arithmetic and outperforms cuBLAS ZGEMM by up to 2.3x on Hopper and 5.5x on Blackwell.

• The paper demonstrates fused kernel designs that minimize memory round-trips, enabling high-precision matrix multiplication on low-precision hardware.
• It introduces two variants, EmuGEMM-I and EmuGEMM-II, which use interleaved data layouts and CRT-based reconstruction to deliver 1.2× to 5.5× throughput improvements over cuBLAS.
• The methodology provides continuous precision-throughput trade-offs, offering a practical solution for HPC applications facing precision constraints.

EmuGEMM: Fused Tensor Core Kernels for Precision Emulation in Matrix Multiplication

Motivation and Architectural Context

Recent GPU architectures, notably NVIDIA Hopper and Blackwell, increasingly optimize for low-precision matrix-multiplication units (MMUs) such as INT8 Tensor Cores (TCs), achieving peak throughputs up to 112$\times$ higher than FP64 on Blackwell relative to previous architectures (Figure 1). The widening precision-throughput gap presents significant challenges for HPC workloads requiring high precision (FP32/FP64), including computational chemistry, numerical linear algebra, and climate modeling, which cannot compromise solution accuracy by simply reducing precision. Emulation schemes, specifically Ozaki Scheme I (mantissa splitting) and Scheme II (modular/CRT-based integer decomposition), are critical for bridging this gap by reconstructing high-precision GEMM from low-precision operations.

Figure 1: GPU data formats and TC characterization, indicating the INT8 throughput advantage and architectural interface evolution on Hopper and Blackwell.

Kernel Fusion and Data Layout Innovations

Existing implementations of Ozaki emulation schemes are profoundly bottlenecked by repeated memory round-trips: intermediate products are frequently materialized and reloaded into GMEM, resulting in substantial operand and accumulator traffic. EmuGEMM addresses this by leveraging kernel fusion and a hardware-aware interleaved data layout, thereby reducing slice loads from $\mathcal{O}(p^2)$ to $\mathcal{O}(p)$ in Scheme I and eliminating INT32 round-trips entirely in Scheme II.

EmuGEMM-I (Scheme I) utilizes an interleaved layout aligning INT8 slices with the TC tile granularity, permitting all $p(p+1)/2$ slice-pair products to be efficiently computed within a single persistent kernel. The triangular MMA scheduling and register-resident accumulator design ensure no intermediate traffic leaves the chip, enabling a shift-reduce epilogue for final reconstruction.

Figure 2: End-to-end EmuGEMM-I pipeline, detailing Ozaki decomposition, operand interleaving, and on-chip triangular MMA execution with shift-reduce epilogue.

EmuGEMM-II (Scheme II) directly fuses modular reduction into the GEMM kernel. Each INT32 product is reduced to INT8 residue on-chip, avoiding the $4{\times}$ output amplification of naive implementations. The CRT-based 3M kernel for complex arithmetic computes all three real products within a fused pipeline, retaining intermediate INT8 results and forming the complex output via error-free integer operations.

Figure 3: EmuGEMM-II kernel structure for modular reduction and fused 3M complex GEMM.

Resource Analysis and Microarchitectural Considerations

Kernel fusion necessitates increased on-chip storage for operand buffering and accumulators. EmuGEMM-I scales per-CTA resource consumption with $p$, requiring efficient management of register (RF/TMEM) and SMEM budgets, especially on Hopper versus Blackwell (where TMEM alleviates RF pressure). Maximizing TC throughput, the effective pipeline depth is quadratic in $p$, ensuring saturation even as per-slice resource budgets shrink.

Figure 4: Mapping of EmuGEMM variants on throughput--intensity plane, demonstrating compute-bound operation and superior utilization relative to naive emulation.

Numerical Performance and Precision Gains

EmuGEMM achieves up to 1,639~Top/s on Hopper (83% INT8 peak) and 3,654~Top/s on Blackwell (81%), outperforming cuBLAS native INT8 by 1.2$\times$ and cuBLAS Scheme I emulation by 1.4$\times$. For large matrices, EmuGEMM-I surpasses native cuBLAS TF32 throughput by up to 1.4$\times$ on Hopper and 1.7$\mathcal{O}(p^2)$0 on Blackwell at comparable accuracy. EmuGEMM-II, utilizing Scheme II modular reduction and CRT reconstruction, outperforms cuBLAS ZGEMM by 2.3$\mathcal{O}(p^2)$1 (Hopper) and 5.5$\mathcal{O}(p^2)$2 (Blackwell).

Figure 5: Effective throughput and precision of real and complex GEMM across architectures, schemes, and slice/modulus counts. EmuGEMM consistently dominates cuBLAS native and emulated baselines.

EmuGEMM-II’s fused CRT approach enables high-precision regimes (e.g., FP64) to not only outperform cuBLAS emulation but native FP64 units on newer architectures dominated by low-precision TCs. The precision-throughput-memory trade-off is characterized by linear workspace scaling in $\mathcal{O}(p^2)$3; Scheme II further incurs output-side INT8 residues for CRT, increasing total workspace but dominating throughput at high precision.

Figure 6: Emulated DGEMM and ZGEMM via 3M, comparing EmuGEMM-II with GEMMul8 variants and cuBLAS references, highlighting throughput and precision superiority.

Figure 7: Throughput–precision–memory trade-off visualization at scale, emphasizing EmuGEMM’s continuous trade-off and workspace scaling by slice/modulus counts.

Theoretical and Practical Implications

The EmuGEMM library advances the state of the art by systematically eliminating the fundamental data movement bottlenecks in floating-point emulation schemes on GPU MMUs. This approach ensures that scientific workloads can fully leverage the underlying hardware throughput. As found empirically, architectures with an increasingly larger INT8-to-FP64 ratio favor emulation over native FP64, and this hardware trajectory is projected to continue. The fused kernel concept, including interleaved data layouts and persistent triangular or modular scheduling, is extensible to newer GPU generations, alternative low-precision floating-point units (FP8/BF16), and potential adaptations on CPUs and specialized accelerators.

From a practical standpoint, the EmuGEMM methodology enables continuous precision-throughput trade-offs spanning TF32 to FP64, allowing applications to select optimal operating points instead of the discrete precision levels offered by native hardware routines. This is particularly critical for memory-constrained workloads favoring Scheme I at moderate precision and for high-precision workloads where Scheme II is preferable despite increased workspace demands.

Limitations and Future Directions

Reported benchmarks focus on square matrices and random datasets with high conditioning; while prior theoretical work validates the Ozaki schemes’ numerical properties, comprehensive studies on irregular problem shapes and real application datasets are required. Currently, EmuGEMM supports only Hopper and Blackwell; porting to other architectures remains a priority for the research community. Automatic dynamic precision (ADP) features from cuBLAS are absent but orthogonal, and could be integrated to further improve runtime adaptability.

Future developments in hardware (e.g., increased TMEM capacity or new low-precision units) and software (e.g., reduction operations for FP8-based emulation) offer promising avenues, with EmuGEMM’s data-movement optimization strategies readily applicable.

Conclusion

EmuGEMM establishes a new baseline for performance and accuracy in GPU-based floating-point emulated matrix multiplication. By fusing kernels and optimizing data movement for Ozaki Schemes I and II, EmuGEMM not only exceeds the throughput of all existing emulation libraries but outperforms native GEMM routines at comparable precision on INT8-dominated microarchitectures. This positions EmuGEMM as a foundational tool for HPC applications facing growing precision-throughput disparity, and suggests that high-fidelity scientific computing will remain viable on low-precision GPU hardware.