FP64 Performance and the Limits of Ozaki Emulation on HPC Systems

FP64 Performance and the Limits of Ozaki Emulation on HPC Systems Nicholas Malaya

Katarzyna (Kasia) Świrydowicz

arXiv:2611.13332v1 [cs.DC] 9 Nov 2026

Advanced Micro Devices Austin, TX kasia.swirydowicz@amd.com

Noah Wolfe

Paul Bauman

Advanced Micro Devices, Inc. Advanced Micro Devices, Inc. Advanced Micro Devices Austin, TX Austin, TX Austin, TX Paul.Bauman@amd.com Noah.Wolfe@amd.com Nicholas.Malaya@amd.com

Pieter Ghysels

Jonathan Belof

Advanced Micro Devices, Inc. Austin, TX Pieter.Ghysels@amd.com

Advanced Micro Devices, Inc. Austin, TX Jonathan.Belof@amd.com

Abstract—A growing industry narrative holds that native double precision can be deprioritized in High Performance Computing (HPC) accelerators, on the hypothesis that low-precision tensor cores can emulate FP64 on demand, and that most HPC workloads are bandwidth-bound. This paper evaluates some of those claims and explores the importance of FP64 in scientific computing. We show roofline profiling data of representative HPC workloads, across a variety of common computational motifs, to evaluate the arithmetic-intensity and accuracy requirements of these workloads. We also consider FP64 emulation via the Ozaki method to ascertain some of the possibilities and limitations of these approaches. These results indicate that Ozaki emulation should not be evaluated as a general replacement for FP64, but as an optimization for the subset of workloads whose dominant computation is sufficiently large, dense, and General Matrix Multiply (GEMM)shaped. Our analysis indicates this is less than 10% of the HPC ecosystem. Furthermore, even in this favorable regime, standalone GEMM speedup is not enough to accelerate (or avoid performance regression) for the application: numerical behavior, workspace requirements, matrix shape, and integration overhead determine whether the benefit survives at application scale. This paper characterizes those boundaries empirically. Furthermore, we show that dense GEMM accounts for only a small fraction of FP64 computation across HPC workloads, leaving most application-level FP64 work beyond the reach of GEMM-based emulation. Index Terms—floating-point arithmetic, double precision, mixed-precision computing, matrix multiplication, numerical accuracy, performance analysis, GEMM, Ozaki scheme I, Ozaki scheme II, emulation

These claims conflate several distinct questions. Whether a kernel is compute or bandwidth-bound depends on the balance of the target architecture, not only on the application itself. Moreover, FP64 computation includes vector operations, reductions, stencils, sparse linear algebra, and applicationspecific kernels, most of which cannot be expressed as dense matrix multiplication. Even when a workload is GEMMdominated, the speed of an emulated GEMM does not by itself determine application performance; numerical behavior, workspace requirements, matrix shape, and integration overhead must also be considered. This paper examines these limitations at both the workload and kernel levels. We first use roofline measurements from 47 configurations spanning eleven HPC applications and benchmarks to characterize the importance and computational form of FP64 work (Section III). Next, we briefly explain the ideas behind two Ozaki schemes (Section IV). We then evaluate Ozaki-style GEMM emulation across three current GPU architectures, considering accuracy (Section V), performance and memory requirements (Section VI), reproducibility (Section VII), and sensitivity to matrix scaling and shape. Together, these results identify where low-precision GEMM emulation is useful and where substantial native FP64 capability remains necessary. Real scientific applications present substantially broader requirements than dense matrix multiplication alone, including programmability, performance portability, communication, and application-specific optimization [1]. Even High Performance Linpack (HPL) contains critical components beyond doubleprecision GEMM (DGEMM) whose implementation can materially affect performance at scale [2].

I. I NTRODUCTION Double-precision floating-point arithmetic has been fundamental to scientific simulation, numerical solvers, and engineering analysis. At the same time, Graphics Processing Unit (GPU)s designed primarily for artificial intelligence increasingly devote more resources to low-precision tensor operations and memory capacity while providing comparatively less native FP64 throughput. This trend is often justified by two claims: that most HPC workloads are limited by memory bandwidth, and that low-precision tensor cores can emulate FP64 computation when greater accuracy is required.

II. P RIOR W ORK The original Ozaki scheme introduced an error-free transformation for expressing high-precision matrix multiplication as a collection of lower-precision matrix products [3]. Mukunoki et al. [4] subsequently adapted the scheme to FP16 Tensor Cores, demonstrating accurate and reproducible DGEMM on GPUs.

Ootomo et al. [5] extended this approach to INT8 matrixmultiplication units, while Ozaki et al. [6] introduced the Chinese Remainder Theorem (CRT)-based Ozaki-II scheme. More recent implementations have explored FP8 Tensor Cores and software-emulated FP64 arithmetic [7], dynamic selection of the required precision [8], and kernel fusion to reduce intermediate memory traffic [9]. Complementary work has examined the numerical and architectural limitations of these approaches. Dongarra et al. [10] survey the broader shift toward reduced- and mixed-precision scientific computing. Abdelfattah et al. [11] analyze integerbased matrix multiplication and show that poorly scaled inputs can reduce accuracy or require additional products, while Uchino et al. [12] provide deterministic error bounds for Ozaki-II. Most recently, a series of papers from Matsuoka [13] advances the broader claim that FP8-based emulation can replace native FP64 across HPC kernels. Our work complements these studies by evaluating both sides of that claim: the accuracy and performance boundaries of standalone Ozaki GEMM, and the limited fraction of application-level FP64 computation that is sufficiently large, dense, and GEMMshaped. We additionally quantify the substantial workspace overhead of Ozaki emulation, which can increase the total High Bandwidth Memory (HBM) footprint to approximately 3× that of native FP64 GEMM, limiting the problem sizes that can reside on the accelerator. III. T HE P REVALENCE OF M ATRIX M ULTIPLICATION IN THE C OMPUTATIONAL S CIENCES To place GEMM emulation in the broader context of HPC workloads, Figure 1 shows the measured FP64 throughput and arithmetic intensity of 47 workloads spanning eleven proxyapplication families. Our roofline dataset contains 47 workload configurations drawn from 11 applications and benchmarks, comprising 131 plotted kernel instances. Table I summarizes the scientific or computational domain represented by each application. Most workloads lie below the HBM ridge point and are therefore primarily bandwidth-sensitive. However, a distinct compute-intensive subset remains sensitive to native FP64 throughput. Moreover, roofline position alone does not establish applicability to Ozaki-style emulation: only workloads whose FP64 computation is dominated by sufficiently large dense GEMMs can directly benefit. Thus, memory bandwidth, native FP64 throughput, and GEMM emulation address complementary rather than interchangeable parts of the HPC workload portfolio. This imbalance may become substantially more pronounced in future accelerators designed primarily for Artificial Intelligence (AI). NVIDIA’s Vera Rubin GPU provides 33 TFLOP/s of native FP64 throughput1 and 22 TB/s of HBM4 bandwidth, while AMD’s Instinct MI455X provides only 5 TFLOP/s of 1 NVIDIA separately reports 200 TFLOP/s of FP64 DGEMM throughput for Rubin using Tensor Core-based emulation. This rate is specific to GEMM and does not represent the native FP64 ceiling for general vector, reduction, sparse, or application-specific kernels.

TABLE I: Applications and benchmarks represented in the MI350X roofline dataset. The dataset comprises 47 application configurations and 131 plotted kernel instances. Application

Scientific or computational domain

XSBench TeaLeaf

Monte Carlo neutron transport and reactor physics Heat conduction using implicit stencil-based solvers Compressible Eulerian hydrodynamics High-order finite-element operator evaluation and PDE discretization High-order finite-element Lagrangian hydrodynamics Lagrangian hydrodynamics on unstructured meshes Classical molecular dynamics and materials modeling Deterministic discrete-ordinates particle transport Dense linear algebra and system-level FP64 benchmarking Sparse iterative linear algebra and system benchmarking Device-memory bandwidth benchmarking

CloverLeaf CEED Bake-Off Laghos PENNANT LAMMPS Kripke rocHPL rocHPCG BabelStream

FP64 throughput with 23.3 TB/s of HBM4 bandwidth [14], [15]. These devices therefore offer considerably less native FP64 compute throughput per unit of memory bandwidth than the MI355X, despite their large increases in memory and lowprecision tensor performance. The corresponding HBM ridge points, computed as the ratio of peak FP64 throughput to peak memory bandwidth, are approximately 1.50 FLOP/byte for Rubin and 0.21 FLOP/byte for MI455X. These values are substantially below the MI355X ridge point shown in Figure 1. Under a first-order roofline projection, many kernels that are bandwidth-bound on MI355X would consequently become compute-bound on Rubin and MI455X. Approximately 50 of the 131 plotted kernel instances would cross the HBM ridge on Rubin and roughly 100 on MI455X. Thus, the observation that current HPC workloads are frequently bandwidth-bound does not imply that native FP64 throughput can be reduced without consequence: changing the machine balance changes the bottleneck and would render many existing HPC codes compute-bound. For FP64 kernels that cannot be expressed as dense GEMM, low-precision tensor-core emulation does not remove this limitation. IV. M ATRIX M ULTIPLICATION AND THE O ZAKI S CHEME A. Ozaki-I and Ozaki-II emulation GEMM is a fundamental Basic Linear Algebra Subprograms (BLAS) operation: given matrices A and B with matching dimensions (and an optional D), compute C = AB + D. Ozaki scheme I [3] and Ozaki Scheme II [6] are used for emulating double precision GEMM using lower precision computing (FP8 or INT8). In Ozaki-I scheme, the matrices PsAare iteratively scaled PsB and split into slices to form A = i=1 Ãi and B = i=1 B̃i , where Ãi and B̃i are scaled slices. Slices without their scaling factors fit into prescribed precision (usually INT8 or FP8,

Fig. 1: Roofline analysis of 47 workloads spanning eleven HPC proxy-application families on AMD Instinct MI350X/MI355X. Each marker represents a measured workload or kernel configuration. The roofline ceilings reflect the attainable throughput imposed by the different levels of the memory hierarchy and by peak FP64 compute throughput. Most workloads are bandwidthsensitive, while a smaller but important subset approaches the compute ceiling.

but FP16 can be used). They are multiplied through, scaled back and added together to form C. This algorithm requires a maximum of sA · sB low-precision GEMMs (some GEMMs, which would not contribute to the result can be dropped). One needs 7-8 such slices to preserve double precision accuracy if using INT8 slices for matrices with large exponent range, see [5]. As mentioned before, Ozaki-II scheme is based on an entirely different principle: it uses CRT. Rows of A and columns of B are scaled by a diagonal matrix with powers of 2 on the diagonal so that the elements are integers; several methods exist for choosing the scaling factors [6], [16]. After the scaling step, one uses a set of s coprime numbers called moduli (typically close to the upper limit of INT8, i.e., 255, 253, etc) and computes s slices based on the modulo arithmetic. The ”matching” slices (originating from the division by the same number) are multiplied in the modulo arithmetic. The result is reconstructed using Garner’s algorithms (CRT reconstruction) and scaled back to FP64. Ozaki-II typically requires less GEMMs than Ozaki-I scheme (s for s moduli). It can also be implemented using FP8 for the GEMM step [17]. To preserve FP64 accuracy, one needs 15-16 moduli (depending on the exponent range). As can be inferred from the above, Ozaki-II scheme requires less GEMMs and is expected to be faster (especially if high

number of slices is needed), and thus it immediately gained popularity in HPC community, which lead to highly optimized open source implementation. Thus we limit our studies on Ozaki-I to its accuracy. B. Terminology By FP64 GEMM or native GEMM we understand GEMM computed using BLAS dgemm function on CPU or GPU with input and output matrices stored in double precision accuracy. We use [format] Ozaki-N k, where [format] is either FP8 or INT8, N is either I or II and k = 2, . . . 16 is the number of slices or moduli, to denote specific Ozaki scheme variants. V. ACCURACY OF STANDALONE EMULATED GEMM A. Test setup In the set of tests in the current and the next sections we use matrices with tightly controlled condition number. Let A be an n × m real-valued matrix. We define the condition number as σmax (A) κ (A) = , (1) σmin (A) where σmin and σmax are the smallest and largest singular values of A, respectively. The condition number is an internal property of the matrix; it does not change with the algorithm or floating point precision in which the operations on the

matrix are carried out. Condition number can be considered a sensitivity measure: the higher it is, the more the matrices are sensitive to small perturbation(s) in the computations. In this work, we only present results for square matrices, in both accuracy and performance experiments. We studied nonsquare matrices, however, the performance results depended only on the internal dimension (the dimension along which inner product is performed): this dimension needed to be sufficiently large to see speedups. This aligns with the square matrix results. In the first set of tests we generate (pseudo-random) matrices and enforce the condition number, see Algorithm 1. In Algorithm 1, qr command performs QR decomposition using dgeqrf followed by dorgqr to ensure orthonormal columns (which guarantees that both QU and QV have condition numbers equal 1). The linspace command in Line 3 generates an array of n evenly spaced real numbers starting at 100 and ending at 10k . Random matrices are generated and stored in FP64; the method given in Algorithm 1 is executed entirely on the GPU by using appropriate library (ROCm or CUDA) functions for QR decomposition and a custom kernel to generate S. A in Step 6 is computed by two standard FP64 dgemms.

Fig. 2: Accuracy of various Ozaki scheme II computations using two (distinct) matrices A and B with prescribed condition number (D = 0). Matrix size 4096 × 4096.

Algorithm 1 Algorithm used to generate test matrices Require: Positive integers m and n; condition number 10k 1: U = rand(m, n) ▷ mean 0.0, std 1.0 2: V = rand(n, n) ▷ mean 0.0, std 1.0

3: S = diag linspace 100 , 10k , n ▷ dim n × n 4: [QU , RU ] = qr(U ) 5: [QV , RV ] = qr(V ) 6: A = QU SQT V 7: return A B. Basic tests To measure accuracy, we generate a series of matrices with condition numbers ranging from 1e0 to 1e17 and multiply them with native FP64 GEMM, INT8 Ozaki-II 12 GEMM, INT8 Ozaki-II 16 GEMM, FP8 Ozaki-II 10 GEMM, and FP8 Ozaki-II 12 GEMM. We use RIKEN’s GEMMul8 code [18], commit 9c1dba457fd5. We measure the accuracy with the formula ∥Cozaki − CFP64 ∥F εF = , (2) ∥CFP64 ∥F where ∥ · ∥F denotes Frobenius norm. The results are shown in Figure 2. The GEMMul8 code allows for two pathways: accurate (shown in figure) and fast (not shown). The difference lies in how the scaling factors in Ozaki-II are computed; the performance difference is generally within 5% especially in high moduli count. One can observe that as long as enough moduli is used, GEMM accuracy stays consistent and in an acceptable range, even if multiplying two extremely ill-conditioned matrices. Using more moduli does not improve the accuracy any further. In further experiments, we used various matrix shapes and sizes,

Fig. 3: Accuracy of Ozaki I scheme computations using two (distinct) matrices A and B with prescribed condition number (D = 0). Matrix size 4096 × 4096. The 6 and 8 slice results differ significantly from the IEEE machine accuracy.

but we omit the plots, as they look analogous. In this simple test, ill-conditioning does not break Ozaki scheme GEMM the results are independent on the condition number and only depend on the number of moduli. We do an analogous test using Ozaki I scheme. In this case we use ozablas code [19], commit 043719e. Note that ozablas implements Ozaki I with INT8 backend and no FP8 backend is available. In case of Ozaki I, less slices are needed; using more slices than shown in Figure 3 does not improve the accuracy any further. The conclusions are identical: the accuracy depends on the number of slices and not on the condition number, These tests do not prove that Ozaki scheme I and Ozaki

Fig. 4: Results of the scaling test using GEMMul8 code.

scheme II themselves are accurate – if the matrix condition number is as large as 1e17, the FP64 results cannot be trusted – they might be very different from the exact (”pen and paper”) analytic result. The tests we performed simply show that native and emulated GEMM produce very similar outputs for similar inputs.

Fig. 5: Pairwise Frobenius-norm differences under the productpreserving scaling As = AD and Bs = D−1 B. Native FP64 remains invariant, whereas Ozaki-II becomes increasingly sensitive as the operand exponent range m grows. Thus, algebraically equivalent matrix products can produce materially different Ozaki II results solely because of operand scaling. TABLE II: System configurations.

C. Test that breaks the accuracy

Name

GPU

CPU

Software Stack

To further stress-test emulation we perform another experiment, based on [20]. We generate random matrices A and B with mean value of 0 and standard deviation of 1, each matrix of size 1024 × 1024. Next, we scale the matrices: As = AD and Bs = D−1 B, where D is a diagonal matrix with entries 2t with t ∈ [−m, m] is an integer chosen randomly and m = 0, 2, . . . , 54. In exact arithmetic, AB = As Bs . We compute six outputs, as follows:

MI355X B200 GB300

AMD MI355X NVIDIA B200 NVIDIA GB300

AMD EPYC 9575F Intel Xeon 6960P Arm Neoverse-V2

ROCm 7.2.1, GCC 13.3.0 CUDA 13.3.73, GCC 13.3.0 CUDA 13.1, GCC 13.3.0

C1 = AB

in FP64,

(3)

C2 = As Bs

in FP64,

(4)

C3 = As Bs

using INT8 Ozaki II-16,

(5)

C4 = AB

using INT8 Ozaki II-16,

(6)

C5 = As Bs

using FP8 Ozaki II-12,

(7)

C6 = AB

using FP8 Ozaki II-12.

(8)

We compare the norms of Ci s on Figure 4. It is clear from the graph that both C3 and C5 start differing from the correct result (C1 ) at m ≈ 34. The norm deviates more severely if INT8 arithmetic is used (line with red triangles). In Figure 5 we compare the resulting matrices pair-wise. If C3 or C5 is used, the difference eventually becomes very large; even modest m decreases accuracy. At the same time, there is no difference between C1 and C2 (flat line at y = 0) and C4 and C6 are good approximations of C1 . This data shows that Ozaki scheme II is sensitive to scaling even if we use a large number of slices. VI. P ERFORMANCE OF STANDALONE EMULATED GEMM In our next set of tests, we evaluate the performance. We use three computational systems, see Table II for details.

According to the theoretical specifications, on Nvidia B200 FP8 and INT8 GEMM achieves up to 4.5 POPS (Peta OPS) and FP64 GEMM achieves 40 TFLOPs, which gives a speedup of 112.5× per single GEMM. Hence, we can expect speedup up to 7.04× for INT8 Ozaki-II 16 and up to 9.38× for FP8 Ozaki-II 12. On B300, we get only 1.3 TFLOPS for FP64 GEMM according to the technical specs [21], 10 PFLOPS per FP8 and 330 TOPs for INT8. The ratio per single GEMM is 254× for INT8 and 7692× for FP8, which gives up to 16× speedup for INT8 Ozaki-II 16 and up to 641× for FP8 Ozaki-II 12. On MI355X, the theoretical peak FP64 matrix performance is 78.6 TFLOPS, the peak INT8 is 5 POPs. and peak FP8 is 5 PFlops. This leads to 63.6 ratio, which caps the speedup at 3.98 for INT8 Ozaki-II 16. These numbers are not achievable in practice (Ozaki II does more than just GEMMs: the matrices need to be preprocessed and the results need to be assembled, and there is memory traffic and cost of kernel launches) but they constitute the theoretical upper limit. Results are shown in Tables III, IV and V. The workspace size was measured with GEMMul8’s own function gemmul8::workSize. Note that for the 40960 × 40960 matrix we got an out-of-memory error on B200. Here is what we can infer from the data: 1. INT8 Ozaki II scheme is faster than native GEMM only for large enough matrices (4096 × 4096 and larger on AMD MI355X and Nvidia B200). At this size, one needs

∼ 850 MB (for Ozaki II-16) of workspace storage. 2. The workspace size grows with the matrix size; for instance, to store 40960 × 40960 A and B (and resulting C) we need about 37.5 GB. However, the workspace size (for INT8 Ozaki-II 16) increases memory requirement three-fold (to about 120 GB). The increase is only larger for FP8 Ozaki II. 3. The speedup achieved by Ozaki is
larger with larger matrix size. GEMM performs O n3 operations (where n is the matrix size); all the supplementary operations scale with n2 , hence at large values of n GEMMs dominate over everything else, including a constant cost of kernel launches that limits the performance at small matrix sizes. 4. While FP8 Ozaki-II 12 leads to about 5× higher speedup (compared to INT8 Ozaki II-16) on GB300, almost 2× as much additional memory is required. 5. On GB300, INT8 Ozaki II never reaches performance of native FP64 delivered by B200. FP8 Ozaki needs very large (8192 × 8192 or larger) matrix to reach it. 6. The speedups (emulated GEMM vs native GEMM) are very high on GB300 mostly because the native FP64 performance is greatly reduced. The workspace reported in Table III is not a minor scratch allocation: it consumes a substantial fraction of GPU HBM. For the 40960 × 40960 case, storing A, B, and C requires approximately 37.5 GB, while Ozaki-II requires an additional 64.0 GB with 12 moduli or 83.2 GB with 16 moduli. The resulting total footprint is therefore approximately 100–120 GB, or 2.6–3.1× that of native FP64 GEMM.

mean that Ozaki scheme is ”more correct” (closer to the analytical true solution): it simply means that for the same input we can expect exactly the same output (correct or not, this test does not check correctness); this is not the case for FP64 computations for high condition numbers, see Figure 6, subplot (e).

VII. ROBUSTNESS OF STANDALONE EMULATED GEMM

Our results demonstrate that native FP64 remains an essential architectural capability for HPC. The roofline analysis shows that FP64 performance, not only memory bandwidth, determines application performance. Reducing a device’s ratio of FP64 throughput to memory bandwidth lowers its roofline ridge point, causing kernels that are bandwidth-bound on current HPC GPUs to become compute-bound on future AIspecialized accelerators. GEMM emulation cannot address this gap broadly. In our dataset, ten of the eleven applications and benchmarks have no dominant DGEMM kernel and therefore depend primarily on native vector FP64 performance. Even when a workload is amenable to Ozaki-style emulation, its practical benefit is constrained by accuracy,workspace, problem size, and integration overheads. Accelerators such as Vera Rubin and MI455X may be highly effective for AI, but their reduced native FP64 throughput prevents them from serving as general replacements for HPC-oriented GPUs across this workload portfolio. These results motivate continued investment in science-focused accelerators such as MI430X, which preserves substantial hardware-based FP64 throughput while also supporting modern AI workloads. Future work will extend this study from standalone GEMM and selected application integrations to a broader set of production HPC workloads. We will investigate how the observed accuracy, scaling, and workspace sensitivities propagate

A. Is Ozaki reproducible? A common misconception about emulation is that emulated results are not bit-wise reproducible: the results of two emulated GEMM operations might not be the same for the same inputs; this is a valid concern as even small variations in the results can lead to unstable behavior in numerically sensitive applications. The published literature clearly disproves this [4], [12] and both Ozaki schemes are reproducible by design. To confirm reproducibility in practice, we run the following experiment: we generate a set of 4096 × 4096 matrices with prescribed condition numbers. We save the matrices and multiply them using INT8 Ozaki-II 12, INT8 Ozaki-II 16, FP8 Ozaki-II 10, FP8 Ozaki-II 12, and FP64 on three compute systems from Table II. We save the computed Cs (five results per matrix pair) and compare them per method. The final results are shown in Figure 6. In all cases of Ozaki II scheme, the results are identical (note here Frobenius norm is not scaled). The results of native FP64 differ for the native GEMM and the difference grows with the growth of the condition number. We conclude that not only Ozaki scheme II produces bit-wise reproducible results: it is actually more numerically reliable than the native GEMM, see Figure 6. The results are identical, regardless of the platform and number of moduli used. It is worth emphasizing that it does not

B. Corner cases GPU-accelerated BLAS libraries such as ROCm rocBLAS [22] perform large GEMMs by splitting the work into tiles on several levels; bottom level tile sizes are aligned with the sizes needed for hardware supported Matrix Fused Multiply-Add (MFMA) operations. These operations act on chunks of data with shapes such as, for example, 16x16x4 (16 × 4 matrix multiplied by 4 × 16 matrix) for MI355X FP64 MFMA and 32x32x32 (32 × 32 matrix multiplied by 32 × 32 matrix) for MI355X INT8 MFMA. When the matrix size is not aligned with the MFMA size, more tiles are launched to cover full matrix multiplication. Hence, we expect a slight performance drop when using misaligned matrix sizes such as 4097 × 4097. We note that Ozaki II performs 12 or 16 GEMMs instead of one, so we expect performance to drop further as the effects of misaligned matrix size accumulate, see Figure 7. We can expect this to be worse for Ozaki I as it executes more GEMMs than Ozaki II. In the figure the performance drops less than 5% for FP64 GEMM and can drop more than 25% for the emulated GEMMs. VIII. C ONCLUSIONS AND FUTURE WORK

(a) FP8 Ozaki II-10

(b) FP8 Ozaki II-12

(c) INT8 Ozaki II-12

(d) INT8 Ozaki II-16

(e) Native FP64

Fig. 6: Comparison of results across three GPUs (AMD MI355X, NVIDIA B200, and NVIDIA GB300).

TABLE III: AMD MI355X: INT8 GEMMul8 Ozaki II performance and workspace for C = AB with A ̸= B (both square, same size) with condition number κ = 103 . Performance is the mean of 50 timed runs after 10 warmup iterations. Native

Ozaki II (12 moduli)

Ozaki II (16 moduli)

Size

TFLOP/s

TFLOP/s

WS (MB)

Speedup

TFLOP/s

WS (MB)

Speedup

384 512 1024 2048 4096 8192 16384 24576 32768 40960

3.94 9.64 38.33 58.53 65.54 68.61 69.13 65.76 66.26 64.36

0.44 1.09 7.79 38.48 86.06 131.85 161.49 173.41 176.10 183.51

40.1 41.8 71.0 188.0 656.0 2560.0 10240.1 23040.1 40960.1 64000.2

0.11× 0.11× 0.20× 0.66× 1.31× 1.92× 2.34× 2.64× 2.66× 2.85×

0.39 0.96 6.66 30.76 66.21 102.55 123.88 129.93 134.36 136.31

42.6 44.8 83.0 236.0 848.0 3328.0 13312.1 29952.1 53248.1 83200.2

0.10× 0.10× 0.17× 0.53× 1.01× 1.49× 1.79× 1.98× 2.03× 2.12×

TABLE IV: Nvidia B200: INT8 and FP8 GEMMul8 Ozaki II performance and workspace for C = AB with A ̸= B (both square, same size) with condition number κ = 103 . Performance is the mean of 50 timed runs after 10 warmup iterations. Native

INT8 Ozaki II (12 moduli)

INT8 Ozaki II (16 moduli)

FP8 Ozaki II (10 moduli)

FP8 Ozaki II (12 moduli)

Size

TFLOP/s

TFLOP/s

WS (MB)

Speedup

TFLOP/s

WS (MB)

Speedup

TFLOP/s

WS (MB)

Speedup

TFLOP/s

WS (MB)

Speedup

384 512 1024 2048 4096 8192 16384 24576 32768

5.65 11.77 27.10 34.20 35.18 35.64 35.94 36.41 36.11

0.61 1.40 9.35 39.57 96.22 147.96 174.70 152.20 161.46

40.1 41.8 71.0 188.0 656.0 2560.0 10240.1 23040.1 40960.1

0.11× 0.12× 0.34× 1.16× 2.74× 4.15× 4.86× 4.18× 4.47×

0.15 1.21 7.74 31.30 73.79 109.42 131.10 119.49 126.60

42.6 44.8 83.0 236.0 848.0 3328.0 13312.1 29952.1 53248.1

0.03× 0.10× 0.29× 0.92× 2.10× 3.07× 3.65× 3.28× 3.51×

0.47 1.11 6.79 25.88 52.03 67.64 75.82 78.62 79.66

48.1 51.5 110.0 344.0 1280.0 5120.0 20480.1 46080.1 81920.1

0.08× 0.09× 0.25× 0.76× 1.48× 1.90× 2.11× 2.16× 2.21×

0.42 1.01 6.07 22.40 44.05 57.05 63.72 65.95 66.77

51.5 55.5 126.0 408.0 1536.0 6144.0 24576.1 55296.1 98304.1

0.07× 0.09× 0.22× 0.66× 1.25× 1.60× 1.77× 1.81× 1.85×

TABLE V: Nvidia GB300: INT8 and FP8 GEMMul8 Ozaki II performance and workspace for C = AB with A ̸= B (both square, same size) with condition number κ = 103 . Performance is the mean of 50 timed runs after 10 warmup iterations. Native

INT8 Ozaki II (12 moduli)

INT8 Ozaki II (16 moduli)

FP8 Ozaki II (10 moduli)

FP8 Ozaki II (12 moduli)

Size

TFLOP/s

TFLOP/s

WS (MB)

Speedup

TFLOP/s

WS (MB)

Speedup

TFLOP/s

WS (MB)

Speedup

TFLOP/s

WS (MB)

Speedup

384 512 1024 2048 4096 8192 16384 24576 32768

0.49 0.88 1.05 1.10 1.10 1.10 1.10 1.10 1.10

0.41 0.94 4.06 7.34 9.83 11.01 11.58 11.77 11.88

40.1 41.8 71.0 188.0 656.0 2560.0 10240.1 23040.1 40960.1

0.85× 1.07× 3.86× 6.69× 8.94× 10.00× 10.50× 10.67× 10.77×

0.35 0.81 3.22 5.76 7.59 8.45 8.87 9.01 9.09

42.6 44.8 83.0 236.0 848.0 3328.0 13312.1 29952.1 53248.1

0.73× 0.92× 3.07× 5.25× 6.90× 7.67× 8.05× 8.17× 8.24×

0.44 0.99 5.17 16.66 35.35 58.73 75.67 84.30 88.08

48.1 51.5 110.0 344.0 1280.0 5120.0 20480.1 46080.1 81920.1

0.90× 1.12× 4.92× 15.19× 32.16× 53.32× 68.65× 76.44× 79.85×

0.39 0.89 4.63 14.69 30.82 50.62 63.75 71.17 74.19

51.5 55.5 126.0 408.0 1536.0 6144.0 24576.1 55296.1 98304.1

0.80× 1.01× 4.41× 13.39× 28.03× 45.95× 57.84× 64.53× 67.26×

Fig. 7: GEMMul8 on AMD MI355X: influence of misaligned matrix sizes on performance

through complete solvers and multi-node runs. We will also explore whether AI-for-science surrogates and other datadriven models can safely exploit reduced- and mixed-precision arithmetic while preserving the accuracy required by scientific workflows. When MI430X becomes available in 2027, we plan to validate our roofline projections and compare native highthroughput FP64 execution, mixed-precision approaches, and Ozaki-style emulation on the same architectures. AI-for-science surrogate models provide another complementary direction, potentially replacing selected high-cost solver components or enabling reduced- and mixed-precision execution while retaining a conventional high-fidelity correction mechanism [23]. IX. R EPRODUCIBILITY All results in this work are meant to be reproduced. We already provided details of hardware and compilers, and the commit hashes of the open source repositories used in this work. The code used in this paper can be found in this repository: https://github.com/kswirydo/ozaki tests/ R EFERENCES [1] N. Malaya, B. Messer, J. Glenski, A. Georgiadou, J. G. Lietz, K. Gottiparthi, M. Day, J. Chen, J. Rood, L. Esclapez, J. W. III, G. R. Jansen, N. Curtis, S. Nichols, J. Kurzak, N. Chalmers, W. Freitag, P. T. Bauman, A. Fanfarillo, R. D. Budiardja, T. Papatheodore, N. Frontiere, D. McDougall, M. Norman, S. Sreepathi, P. C. Roth, D. Bykov, N. Wolfe, P. Mullowney, M. Eisenbach, M. T. H. de Frahan, and W. Joubert, “Experiences readying applications for exascale,” in Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, ser. SC ’23. Association for Computing Machinery, 2023, pp. 1–13. [Online]. Available: https://doi.org/10.1145/3581784.3607065 [2] H. Lu, M. Matheson, N. Chalmers, A. Kashi, N. Malaya, and F. Wang, “Insights from optimizing HPL performance on exascale systems: A comparative analysis of panel factorization,” in Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, ser. SC ’25. Association for Computing Machinery, 2025, pp. 397–410. [Online]. Available: https://doi.org/10.1145/3712285.3759875 [3] K. Ozaki, T. Ogita, S. Oishi, and S. M. Rump, “Error-free transformations of matrix multiplication by using fast routines of matrix multiplication and its applications,” Numerical Algorithms, vol. 59, no. 1, pp. 95–118, 2012. [4] D. Mukunoki, K. Ozaki, T. Ogita, and T. Imamura, “DGEMM using tensor cores, and its accurate and reproducible versions,” in High Performance Computing: 35th International Conference, ISC High Performance 2020, ser. Lecture Notes in Computer Science, vol. 12151. Springer International Publishing, 2020, pp. 230–248. [5] H. Ootomo, K. Ozaki, and R. Yokota, “DGEMM on integer matrix multiplication unit,” The International Journal of High Performance Computing Applications, vol. 38, no. 4, pp. 297–313, 2024. [6] K. Ozaki, Y. Uchino, and T. Imamura, “Ozaki Scheme II: A GEMMoriented emulation of floating-point matrix multiplication using an integer modular technique,” The International Journal of High Performance Computing Applications, 2025. [7] D. Mukunoki, “DGEMM using FP64 arithmetic emulation and FP8 tensor cores with Ozaki scheme,” in Proceedings of the Supercomputing Asia and International Conference on High Performance Computing in Asia Pacific Region Workshops. Association for Computing Machinery, 2026, pp. 303–311. [Online]. Available: https://doi.org/10. 1145/3784828.3785017 [8] A. Schwarz, A. Anders, C. Brower, H. Bayraktar, J. Gunnels, K. Clark, R. G. Xu, S. Rodriguez, S. Cayrols, P. Tabaszewski, and V. Podlozhnyuk, “Guaranteed DGEMM accuracy while using reduced-precision tensor cores through extensions of the Ozaki scheme,” 2025, sCA/HPCAsia 2026. [Online]. Available: https://arxiv.org/abs/2511.13778

[9] D. Lu, A. Maeder, M. Luisier, and A. N. Ziogas, “EmuGEMM: Fused tensor core kernels for precision emulation in matrix multiplication,” 2026. [Online]. Available: https://arxiv.org/abs/2606.25453 [10] J. Dongarra, J. Gunnels, H. Bayraktar, A. Haidar, and D. Ernst, “Hardware trends impacting floating-point computations in scientific applications,” 2024. [Online]. Available: https://arxiv.org/abs/2411. 12090 [11] A. Abdelfattah, J. Dongarra, M. Fasi, M. Mikaitis, and F. Tisseur, “Analysis of floating-point matrix multiplication computed via integer arithmetic,” arXiv preprint arXiv:2506.11277, 2025. [Online]. Available: https://arxiv.org/abs/2506.11277 [12] Y. Uchino, K. Ozaki, and T. Imamura, “Error analysis of matrix multiplication emulation using Ozaki-II scheme,” arXiv preprint arXiv:2602.02549, 2026. [Online]. Available: https://arxiv.org/abs/2602. 02549 [13] S. Matsuoka, “FP8 is all you need (part 1): Debunking hardware FP64 as the HPC holy grail,” 2026. [Online]. Available: https: //arxiv.org/abs/2606.06510 [14] NVIDIA Corporation, “NVIDIA Vera Rubin NVL72,” https://www. nvidia.com/en-eu/data-center/vera-rubin-nvl72/, 2026, preliminary specifications; accessed: 2026-08-07. [15] Advanced Micro Devices, Inc., “AMD Instinct MI455X GPU,” https: //www.amd.com/en/products/accelerators/instinct/mi400/mi455x.html, 2026, accessed: 2026-08-07. [16] S. Kawakami and D. Takahashi, “Improved scaling for fast mode of Ozaki scheme II,” arXiv preprint arXiv:2606.29129, 2026. [17] Y. Uchino, K. Ozaki, and T. Imamura, “Double-precision matrix multiplication emulation via Ozaki-II scheme with FP8 quantization,” 2026. [Online]. Available: https://arxiv.org/abs/2603.10634 [18] Y. Uchino, “GEMMul8,” https://github.com/RIKEN-RCCS/GEMMul8, Jul. 2026. [19] I. Nyrkov, “ozablas,” https://github.com/IlyaNyrkov/ozablas, Jul. 2026. [20] J. Demmel, X. Li, J. Langou, W. Pereira, M. Gates, and C. R. Gonzalez, “How to grade the accuracy of an implementation of the BLAS; short update on exception handling,” Paper presented at the BLIS Retreat 2024, Austin, TX, September 2024. [Online]. Available: https://www. cs.utexas.edu/∼flame/BLISRetreat2024/slides/Grading BLAS.pdf [21] NVIDIA Corporation, “NVIDIA Vera Rubin NVL72,” https://www. nvidia.com/en-us/data-center/gb300-nvl72, 2026, official specifications; accessed: 2026-08-23. [22] Advanced Micro Devices, Inc., “rocBLAS: Next-generation BLAS implementation for the ROCm platform,” https://github.com/ROCm/ rocm-libraries, 2026. [23] O. Obiols-Sales, A. Vishnu, N. Malaya, and A. Chandramowlishwaran, “CFDNet: A deep learning-based accelerator for fluid simulations,” in Proceedings of the 34th ACM International Conference on Supercomputing, ser. ICS ’20. Association for Computing Machinery, 2020, pp. 1–12. [Online]. Available: https://doi.org/10.1145/3392717. 3392772

Related documents

• FP8 is All You Need (Part 1): Debunking Hardware FP64 as the HPC Holy Grail (Sep 3rd version)#986869
• Compiler and Hardware Co-Design for Accelerator Architectures#1078170
• Compiler-Assisted Performance Analysis Workflows for HPC Applications / Sebastian Kreutzer ; Gutachter: Christian Bischof, Jesus Labarta Mancho#1145408
• Proxics: an efficient programming model for far memory accelerators#120579
• PATHA: Performance Analysis Tool for HPC Applications#611942