Concurrent Use of FP32 and FP64 Computing Units on GPUs and Mixed-Precision Placement in WENO Reconstruction

arXiv:2611.09651v1 [cs.DC] 6 Nov 2026

Concurrent Use of FP32 and FP64 Computing Units on GPUs and Mixed-Precision Placement in WENO Reconstruction Jun Hatayama

Daichi Mukunoki

Yusuke Onishi

Graduate School of Environmental, Life, Natural Science and Technology, Okayama University Okayama, Japan jhatayama@s.okayama-u.ac.jp

Information Technology Center, Nagoya University Nagoya, Japan mukunoki@cc.nagoya-u.ac.jp

Graduate School of Environmental, Life, Natural Science and Technology, Okayama University Okayama, Japan p62d9qiy@s.okayama-u.ac.jp

Atsushi Sekimoto

Toshinori Kouchi

Faculty of Environmental, Life, Natural Science and Technology, Okayama University Okayama, Japan asekimoto@okayama-u.ac.jp

Faculty of Environmental, Life, Natural Science and Technology, Okayama University Okayama, Japan kouchi@okayama-u.ac.jp

Abstract

Keywords

Scientific computing, including computational fluid dynamics (CFD), often requires FP64 arithmetic, although FP64 performance is severely limited on consumer- and workstation-class GPUs, and the shift of GPU resources toward low-precision units for AI workloads may extend this limitation to data-center GPUs. Double-float (DF) emulation, which represents a value using two FP32 values with about 48 bits of significand, outperforms native FP64 where FP64 units are limited. On the other hand, mixed-precision strategies offload to FP32 the operations that do not determine accuracy. In both approaches, the native FP64 units remain idle during FP32 instruction execution because they are physically separate from the FP32 units. We investigated whether this idle resource can be recovered on A100, GH200, GB200, and RTX PRO 6000 Blackwell Max-Q. A microbenchmark shows that co-resident FP32 and nativeFP64 kernels exceed the throughput of the faster kernel alone on A100 but not on GH200 and GB200; on RTX 6000 Blackwell, they do so against FFMA but, against DF emulation, only when the issue density of the DF kernel is reduced. We then tested whether this gain carries over to computing the WENO-Z reconstruction, a standard high-order shock-capturing scheme in CFD, with the nonlinear weights in FP32 and the candidate polynomials in native FP64 (W32P64). The weights are candidates for lower precision, although the precision they need depends on the order, problem, and implementation. We observe substantial overlap on A100, whereas GH200, GB200, and RTX 6000 Blackwell show little overlap. Importantly, the fastest placement does not always preserve the design accuracy of WENO-Z: for a one-dimensional entropy wave problem, W32P64 fails to achieve the design order of accuracy for WENO7 and WENO9. In contrast, DF emulation of the entire reconstruction comes close to it on the grids tested for WENO5, WENO7, and WENO9 while delivering a 2.51–3.94× speedup over FP64 on RTX 6000 Blackwell. For WENO5 on the two-dimensional Euler vortex problem, at the time step used, the W32P64 error at 𝑛 = 256 is 1.26 times that of FP64 and 1/104 of that of the all-FP32 reconstruction.

GPU, WENO-Z, mixed precision, double-float, instruction issue bandwidth, concurrent execution

1

Introduction

Driven by the growing demand from deep learning and AI workloads, GPU vendors have devoted an increasing share of compute resources to low-precision matrix units, while FP64 capability has been steadily deprioritized, generation after generation. On consumer- and workstation-class NVIDIA GPUs, FP64 throughput is only 1/64 of FP32 throughput. However, scientific computing, including computational fluid dynamics (CFD), often requires FP64 arithmetic to capture high-fidelity physics. Two approaches address this mismatch. One is mixed-precision computing, which reduces precision only in the operations that do not govern the accuracy of the solution [3, 8, 9, 14]. The other is high-precision emulation, which represents a high-precision value with multiple low-precision numbers. This paper focuses on doubleword arithmetic [4, 6], which represents a value as the unevaluated sum of two FP32 numbers, giving about 48 bits of significand. On GPUs with limited FP64 units, DF emulation can deliver several times the effective throughput of native FP64. Under either approach, the same thing happens: while an FP32 instruction stream runs, the physically separate native FP64 units sit idle. This paper asks whether that idle resource can be recovered. The ultimate form of recovery is intra-kernel: placing an FP32 instruction and an FP64 instruction side by side in one kernel and letting them overlap in time. In a real compute kernel, however, whether overlap occurs depends on complex interactions among memory access, register pressure, and instruction dependencies. Consequently, when overlap fails to appear, it is difficult to determine whether hardware limitations or software bottlenecks are responsible. We therefore begin with an inter-kernel microbenchmark that runs an FP32-only kernel and a native-FP64-only kernel in separate CUDA streams on the same SM and checks whether the FP32 and FP64 pipes make progress concurrently. Warps from different kernels neither depend on each other nor share registers, so this experiment isolates the hardware-side condition: whether FP32 1

Jun Hatayama, Daichi Mukunoki, Yusuke Onishi, Atsushi Sekimoto, and Toshinori Kouchi

Í Eq. (2), 𝜏/(𝛽𝑘 + 𝜖) and 1/ 𝑚 𝛼𝑚 , the weight computation; the evaluation of the candidate polynomials 𝑝𝑘 , the polynomial computation; Í and the final 𝑘 𝜔𝑘 𝑝𝑘 , the blending. The baseline mixed-precision placement, W32P64, computes the weights in FP32 and keeps the polynomial computation and the blending in FP64.

and FP64 instructions can be issued concurrently on one SM. Our working hypothesis is that, on a GPU where inter-kernel overlap gives a throughput gain, independent FP32 and FP64 instructions within one kernel may overlap as well. As a CFD application, we choose the WENO-Z reconstruction [1], an improved variant of the weighted essentially non-oscillatory (WENO) scheme [7, 10], which is a standard high-order shockcapturing scheme in CFD and is widely used for supersonic and hypersonic flows. The reconstruction builds candidate polynomials from local stencils and blends them with nonlinear weights based on smoothness indicators. The weights act as ratios that select among stencils and are therefore candidates for lower precision than the polynomial values, although the precision they need depends on the order, problem, and implementation. The reconstruction thus splits naturally into two computations that can be given different precisions. Field et al. [5] exploited this split for black-hole simulations, computing the weights in FP64 and the polynomials in FP128. We apply the same split to recover idle hardware instead, computing the weights in FP32 and the polynomials in native FP64. The WENO-Z kernel adds what the microbenchmark removed: its two computations share registers, read the same stencil values, and are ordered by the compiler, so any overlap must survive register pressure, occupancy and instruction order. This paper treats mixed precision not only as a reduction in arithmetic precision but as a redistribution of work across distinct GPU execution resources, whose benefit depends jointly on the microarchitecture and on whether the scheme keeps its design order of accuracy.

2.2

On GPUs with very low native FP64 throughput, W32P64 leaves the polynomial computation and the blending on the bottleneck FP64 pipe. We therefore use double-float (DF) emulation, which achieves FP64-equivalent precision using FP32 operations. DF emulation represents a value as the unevaluated sum of two FP32 numbers, 𝑎 = 𝑎 hi + 𝑎 lo . From the Dekker–Knuth error-free transformations [4, 6], a DF addition takes 11 FP32 instructions, and a DF multiplication takes 7; thus one DF multiply–add (two FP64-equivalent flops) takes 18 FP32 instructions. Because the compiler’s FMA contraction breaks the error-free transformations, we write addition, subtraction, and multiplication with rounding-mode-specified intrinsics, which the compiler does not contract. The resulting significand of approximately 48 bits is shorter than the 53 bits of FP64, and the exponent range remains that of FP32. Hereafter, “FP64equivalent” refers to this precision and operation count; it does not guarantee the exponent range, subnormal handling, or exception behavior of IEEE 754 binary64.

2.3

We use the WENO-Z scheme [1, 2] and call its fifth-, seventh-, and ninth-order variants WENO5, WENO7, and WENO9, respectively. The reconstructed positive/negative flux components at the interface 𝑖 + 1/2 are the weighted sum of candidate polynomials 𝑝𝑘 and nonlinear weights 𝜔𝑘 , 𝑟 −1 ∑︁

𝜔𝑘± 𝑝𝑘±,

Instruction issue and FP32/FP64 execution lanes

An SM consists of four sub-partitions, and the warp scheduler of each sub-partition issues one warp instruction per cycle [11–13]. The number of FP32 and FP64 lanes per sub-partition is 16 and 8 for A100, 32 and 16 for GH200 (Hopper) and GB200 (data-center Blackwell), and 32 and 0.5 for RTX PRO 6000 Blackwell Max-Q (GB202; hereafter RTX 6000 Blackwell), respectively. With a warp width of 32, an FP32 warp instruction occupies the 16-lane pipe of A100 for two cycles and an FP64 warp instruction occupies the 8-lane pipe for four cycles. An FP32 instruction stream on A100 therefore leaves every other issue slot free, and an FP64 stream needs even fewer slots. In terms of issue bandwidth, the two can share the scheduler. On GH200 and GB200, an FP32 stream can issue every cycle, leaving, in principle, no slot for FP64 instructions. On RTX 6000 Blackwell, FP32 can also issue every cycle, but the FP64 pipe accepts only one instruction every 64 cycles, so the two can share if the FP32 stream leaves as few as one issue slot in 64 unused. This is why we expect A100 to behave differently from GH200 and GB200, despite their identical nominal FP64:FP32 ratio of 1:2. We treat this expectation as a hypothesis: whether an instruction of the other pipe can actually use a free slot depends on dispatch ports, operand collectors, register banks, and the fairness of warp selection, which are not fully documented in publicly available NVIDIA materials. Lane counts therefore provide only a necessary condition, and Section 5.4 reports the corresponding warp-scheduler counters.

2 Background and Design Model 2.1 WENO-Z reconstruction and its decomposition

± 𝑞𝑖+1/2 =

Double-float emulation

(1)

𝑘=0

where 𝑟 = 3, 4, 5 corresponds to WENO5, WENO7, and WENO9, and the superscripts − and + denote the left- and right-biased reconstructions, respectively. The candidate polynomials, 𝑝𝑘± , and the smoothness indicators, 𝛽𝑘 , take the standard form [7]: 𝛽𝑘 is a sum of squared derivatives of the 𝑘-th candidate polynomial and is a quadratic form in differences between stencil values. WENO-Z additionally uses a global smoothness indicator 𝜏, (  2) 𝜏 𝛼𝑘 𝛼 𝑘 = 𝑑𝑘 1 + , 𝜔𝑘 = Í , (2) 𝛽𝑘 + 𝜖 𝑚 𝛼𝑚 where 𝑑𝑘 are the linear weights, 𝜖 is a small positive constant, and 𝜏 = |𝛽 0 − 𝛽 2 | for WENO5 [1] and the absolute value of a linear combination of the 𝛽𝑘 for WENO7 and WENO9 [2]. Hereafter we call the evaluation of 𝛽𝑘 , 𝜏, 𝛼𝑘 , and 𝜔𝑘 , including the two divisions in

2.4

Inter-kernel and intra-kernel overlap

FP32 and FP64 instructions can overlap in the following ways. Interkernel overlap: two kernels in separate CUDA streams co-reside on 2

Concurrent Use of FP32 and FP64 Computing Units on GPUs and Mixed-Precision Placement in WENO Reconstruction

one SM, one issuing FP32 and the other FP64 instructions. Because their warps neither depend on each other nor share registers, the requirement reduces to the hardware-side condition: both kernels hold SM resources simultaneously, and the warp scheduler must keep selecting warps from both kernels. A stream that issues on nearly every cycle can leave no slots for the other, so the issue density of each stream, which we set through the length of its dependency chains, is part of this condition even though it is not a hardware property. Intra-kernel overlap: FP32 and FP64 instructions are issued within one kernel, and the overlap comes either from independent instructions within a thread (instruction-level parallelism, ILP) or from other warps (thread-level parallelism, TLP). To yield ILP, the instructions of different types must be interleaved in the machine code of a thread (Section 3.3). On the other hand, TLP requires that the different types of instructions be issued from different warps. Because both paths require free issue slots, we do not expect intra-kernel overlap on a GPU where inter-kernel overlap gives no throughput gain, and software-side conditions can prevent it.

overlap. Note that a DF multiply–add consists of 3 FFMA, 1 FMUL, and 14 FADD/FSUB instructions.

3.3

3 Evaluation method 3.1 Evaluation environment Table 1 lists the four GPUs. A100, GH200, and GB200 share the nominal FP64:FP32 ratio of 1:2. RTX 6000 Blackwell, with a measured ratio of 1:67.3 (nominal 1:64), is bound by native FP64, and therefore serves to evaluate DF emulation. On GH200 and GB200 we also evaluated asynchronous stencil-tile transfer with the Tensor Memory Accelerator (TMA) (Section 5.2).

3.2

WENO-Z kernel and static SASS analysis

For CFD applications, we implement the WENO-Z reconstruction of Eqs. (1)–(2) as a kernel. In FP32, the two divisions in the weight computation use the approximate reciprocal, refined by Í one Newton step for the weight normalization, 1/ 𝑚 𝛼𝑚 . FP64 keeps IEEE division: the FP32 approximate reciprocal (MUFU.RCP) is one instruction close to FP32 precision, whereas the FP64 one (MUFU.RCP64H) is only a seed that needs DFMA Newton steps in any case. -use_fast_math makes the same split. Every placement on every GPU is measured with one hand-written CUDA C kernel, which the solvers of Section 6 also call, and each speedup is against its native FP64 in the same run. We examine the streaming assembler (SASS), the native machine code in CUDA binaries, to determine the instruction order and mix without executing it (static analysis). Because a warp issues the instructions of a thread in program order, a thread can supply overlap only where FP32 and FP64 instructions alternate; if they were grouped, overlap could come only from warps in different phases. In the sm_80 (A100) SASS of the natural version 𝑠 nat of the kernel measured as 𝑠 pure (Section 3.4), the instruction class switches between FP32 and FP64 34 and 93 times at WENO5 and WENO9, respectively, compared with 3 times in its forced-serialization version. In the sm_120 (RTX 6000 Blackwell) SASS, it switches 12 and 30 times. The FP64 instruction counts of the sm_80, sm_90 (GH200), sm_100 (GB200) and sm_120 binaries agree exactly (61, 127, and 187 for WENO5, WENO7, and WENO9), the FP32 counts agree to within 2.3%, and no register spills occur in the measured kernels; the out-of-line routine of the FP64 IEEE division spills 4 bytes on sm_90 and sm_120. Timing differences between generations are therefore not explained by a large difference in instruction counts.

Inter-kernel microbenchmarks

Three microbenchmarks share a common harness. Each kernel launches one block per SM. The two kernels are launched simultaneously from separate CUDA streams and timed with CUDA events, with the iteration count calibrated so that each kernel runs about 300 ms standalone. Each block records its SM ID and its start and end times using clock64(). We define the retention of each kernel as its throughput during the interval in which the execution windows of the two kernels overlap on the same SM, divided by its standalone throughput. The source code and compiler options are identical on the four GPUs except for the -arch option; the CUDA versions differ and are listed in Table 1. Pure kernels (Section 4.1): an FFMA kernel with eight independent accumulators per thread runs against a DFMA kernel of the same structure (25–28 registers per thread, no shared memory). We sweep the FP32-side configuration and the FP64-side warp count. DF emulation against native FP64 (Section 4.2): the FP32 side is DF emulation, counted in FP64-equivalent flops (one DF operation as one flop), against a native FP64 kernel fixed at 1024 threads/SM (32 warps/SM), the configuration of its standalone peak. Two configurations set the issue density of the DF side, both at 256 threads/SM: eight accumulators per thread, which keeps eight independent chains, and two accumulators per thread, for which data dependencies leave gaps in the issue stream. Instruction-type probe (Section 4.3): this measurement repeats DFMA against a single instruction type such as FFMA, FADD, FMUL, or DF emulation, to determine which instruction types prevent

3.4

Measuring overlap in the WENO-Z kernel

In an ordinary kernel, memory access and arithmetic proceed concurrently, and the run time alone cannot separate their contributions. We therefore use a compute-repeat configuration: the stencil input is read from global memory once, the reconstruction is repeated 𝑁 times on the same input, and the result is written back once. The slope of a least-squares linear fit of run time against 𝑁 gives the arithmetic time per reconstruction, excluding memory access. Let 𝑠 32 and 𝑠 64 denote the slope of the weight computation alone and the polynomial computation alone, respectively, and 𝑠 pure the slope of the kernel that performs both computations. We define the overlap indicator as 𝑠 32 + 𝑠 64 − 𝑠 pure 𝐸= . (3) min(𝑠 32, 𝑠 64 ) The numerator is the time saving relative to serial execution, 𝑠 32 +𝑠 64 , and the denominator is its value when the shorter stream is fully hidden; 𝐸 = 0 corresponds to fully serial execution, and 𝐸 < 0 to a kernel slower than serial execution. Because the partial kernels may differ from the combined kernel in register usage and compiler optimization, 𝐸 is an indicator of overlap, not a measurement of the overlap time. When one stream is much longer than the other, the denominator is small, and a timing change of a few percent of 3

Jun Hatayama, Daichi Mukunoki, Yusuke Onishi, Atsushi Sekimoto, and Toshinori Kouchi

Table 1: GPUs used and their main characteristics. Measured peaks are the standalone values (FP32 / FP64) from the pure kernels of Section 4.1.

Architecture / compute capability SMs FP32:FP64 lanes per SM Measured peak [TFLOPS] (FP32 / FP64) Measured FP64:FP32 CUDA / driver

A100 80GB PCIe GH200 120GB

GB200

RTX PRO 6000 Blackwell Max-Q

Ampere / 8.0 108 64:32 19.10 / 9.424 1:2.0 12.5 / 570.211.01

Blackwell / 10.0 152 128:64 65.40 / 38.696 1:1.69 13.3 / 610.57.04

Blackwell GB202 / 12.0 188 128:2 99.80 / 1.483 1:67.3 13.2 / 595.58.03

Hopper / 9.0 132 128:64 52.63 / 31.679 1:1.66 13.0 / 595.58.03

the longer stream moves 𝐸 by a large amount; we therefore report 𝑠 64 /𝑠 32 alongside 𝐸. A positive 𝐸 does not distinguish intra-thread ILP from TLP supplied by other warps. To separate them, we compare against a forced-serialization kernel that removes only the intra-thread overlap while matching the instruction count, bitwise results, and occupancy. In the serialized version (slope 𝑠 ser ), the input of the polynomial computation depends on the output of the weight computation; the independent version (slope 𝑠 ind ) adds the same dependent instructions but applies them to an input that is already final. The added instructions appear in both versions, so their cost cancels in the difference 𝑠 ser − 𝑠 ind , the only quantity we use.

of native FP64, whereas on A100, GH200, and GB200, native FP64 achieves about nine times the throughput of DF emulation. On RTX 6000 Blackwell, when native FP64 runs concurrently with the ordinary DF kernel with eight accumulators per thread, the DF emulation side retains 99.2%, while the FP64 side almost stops during the overlap. We call this state starvation: DFMA instructions obtain almost no issue opportunities, and their effective throughput falls to nearly zero. Reducing the accumulator count from eight to two, so that data dependencies leave natural gaps in the issue stream of the DF kernel, allows DF emulation and FP64 to retain 87.2% and 81.6%, respectively, reaching 5.40 FP64-equivalent TFLOPS, 6.7% above DF emulation alone. On A100, both configurations retain 87–88% on the DF side and 92–94% on the FP64 side, and the eight-accumulator configuration performs best at 9.795 TFLOPS, 3.8% above the 9.441 TFLOPS of native FP64 alone. Concurrency itself thus works well. However, DF emulation alone reaches only about 1/9 of native FP64 throughput, so adding independent DF work yields only a small gain. GH200, despite an FP64:FP32 ratio close to that of A100, drops to 47–63% retention on both sides, and in every configuration the combined throughput falls at least 40.2% below that of native FP64 alone. GB200 falls further, to 30.3% and 52.4% of native FP64 alone with eight and two accumulators per thread. Without interference, the gain over the faster standalone kernel is bounded by the throughput ratio of the slower kernel to the faster one: about 29% on RTX 6000 Blackwell (6.7% measured), 11.3% on A100 (3.8% measured), 10.8% on GH200, and 10.7% on GB200 (negative measured gain on both). Gains must therefore be judged against the faster standalone kernel, not against co-residency or wall-time reduction. We refer to a combined throughput exceeding that of either standalone kernel as an inter-kernel throughput gain, distinct from temporal overlap.

4

Microbenchmark: concurrent FP32/FP64 use across kernels 4.1 Pure FFMA/DFMA kernels We first measure pure compute kernels, which most directly test the expectation of Section 2.3. Across the sweep of the FP32-side configuration and the FP64-side warp count, blocks of both kernels co-resided on every SM of all four GPUs, and their execution windows overlapped completely. Over 128–1024 FP32 threads and 4–16 FP64 warps per SM, on A100, the FP32 side retains 89–92% and the FP64 side 77–89%, saturating at 12 warps/SM. In the best configuration (256 FP32 threads, 384 FP64 threads per SM), the two kernels sustain 90.6% and 89.1%, respectively, and the wall time is 44.0% shorter than when the two kernels run back-to-back. On RTX 6000 Blackwell, FP64 retention likewise recovers as the FP64 warp count increases, from 40.2% at 4 warps to 94.8% at 16, and the best configuration shortens the wall time by 47.1%. GH200 behaves differently: over the same range, FP64 retention remains at only 15–65%, while the FP32 side keeps 79–97%, and even the best configuration (128 FP32 threads, 128 FP64 threads per SM) reduces the wall time by only 27.9%. GB200 behaves like GH200: over the same range, FP64 retention stays at 18–37% while the FP32 side keeps 89–92%, and the best configuration, the same as GH200’s, reduces the wall time by 15.8%. This pattern agrees with the lane-count expectation of Section 2.3.

4.2

4.3

Instruction-type probe

The essential difference between the DF emulation kernel and the pure FFMA kernel is the instruction mix. We therefore isolated its effect by changing only the instruction type of the competing kernel (Table 3). On RTX 6000 Blackwell, GH200, and GB200, saturated singletype FADD or FMUL streams block DFMA almost completely (FP64 retention of 0.1–0.4%). On A100, the same probes keep FP64 retention at 87–93% for every instruction type, and no starvation is observed; with the competing kernel fixed, varying the FP64 warp count from 8 to 32 changes the retention by only 1–2 points.

DF emulation and native FP64

Table 2 gives the standalone throughput of each kernel and the results of concurrent execution in two streams. On RTX 6000 Blackwell, the faster DF configuration achieves 3.4 times the throughput 4

Concurrent Use of FP32 and FP64 Computing Units on GPUs and Mixed-Precision Placement in WENO Reconstruction

Table 2: Concurrent execution of DF emulation and native FP64 in two streams, in FP64-equivalent TFLOPS. The standalone columns give each kernel run alone in the same configuration. The DF side and FP64 side columns are the throughput during the overlap interval, with the retention relative to standalone in parentheses. The last column is the combined throughput divided by the largest standalone value of the GPU; a configuration below 1 is faster without concurrency. Standalone [TFLOPS] GPU

Configuration

A100

8 accumulators 1.071 2 accumulators 1.071

9.441 0.943 (88.0%) 9.441 0.933 (87.1%)

8.852 (93.8%) 8.686 (92.0%)

9.795 9.619

1.038 1.019

GH200

8 accumulators 2 accumulators

31.43 31.43

17.07 (54.3%) 14.66 (46.6%)

18.80 16.80

0.598 0.535

GB200

8 accumulators 4.130 2 accumulators 4.105

38.695 3.344 (81.0%) 8.392 (21.7%) 38.695 2.645 (64.4%) 17.619 (45.5%)

11.736 20.264

0.303 0.524

∼0 (starved) 1.215 (81.6%)

5.03 5.40

0.994 1.067

RTX 6000 8 accumulators 2 accumulators

DF

Concurrent [TFLOPS]

FP64

3.03 3.38

5.06 4.80

1.49 1.49

DF side

1.74 (57.3%) 2.14 (63.3%)

5.02 (99.2%) 4.18 (87.2%)

Table 3: Retention of the native-FP64 kernel (32 warps/SM) during the overlap interval with a competing kernel (256 threads/SM). Competing stream

A100 GH200 GB200 RTX 6000

Continuous FFMA Continuous FADD Continuous FMUL Alternating FMUL/FADD DF (ordinary) DF (all-FMA) Low-density FADD (1 acc.)

89.0% 90.7% 89.4% 88.3% 93.0% 89.4% 87.2%

20.3% 0.4% 0.1% 0.3% 54.1% 51.0% 67.9%

18.8% 0.3% 0.3% 0.3% 41.4% 40.9% 71.3%

5

FP64 side Combined vs. faster standalone

Mixed-precision reconstruction inside the WENO-Z kernel

This section implements the reconstruction of Section 2.1 as a single kernel and measures three precision placements, FP64, W32P64, and DF emulation. We measure DF emulation only on RTX 6000 Blackwell, because on A100, GH200, and GB200 native FP64 achieves about nine times the throughput of DF emulation (Section 4.2). Table 4 highlights two key points. First, W32P64 is faster than FP64 on all GPUs, with speedups of 2.73–3.04×, 1.94–2.78×, and 1.93–2.88× on A100, GH200, and GB200, respectively. With IEEE division in the FP32 weights, they drop to 1.82–2.25×, 1.35–1.59×, and 1.29–1.62×. On RTX 6000 Blackwell, W32P64 gives the largest speedup, 4.47–5.92× (4.54–5.89× with IEEE division), but, as Section 6 shows, it loses the design order of accuracy for WENO7 and WENO9; the fastest placement that preserves the design order is DF emulation, with a speedup of 2.51–3.94×. These speedups also include the changes in arithmetic and division cost and in register usage, so they do not isolate the contribution of overlap. Second, the temporal overlap of the arithmetic depends strongly on the GPU generation, and this tendency is consistent with the results of the microbenchmarks in Section 4.

93.2% 0.4% 0.4% 0.4% 0.6% 75.4%† 100.1%

On RTX 6000 Blackwell, FADD blocks DFMA but FFMA does not (0.4% against 93.2%). Because 14 of the 18 instructions of a DF multiply–add are FADD/FSUB, this is consistent with the starvation of native FP64 by DF emulation, but not by FFMA, on that GPU. We also rewrote every addition and multiplication of the DF kernel as an FMA, which preserves rounding and instruction count. On RTX 6000 Blackwell, this rewrite raised FP64 retention to 75.4% (marked † in Table 3). The same transformation embedded in a different program, however, starved the FP64 kernel completely again, and on A100 and GH200 it lowered the retention, while on GB200 it changed it by only 0.5 points. The rewrite is therefore not a portable remedy. In summary, an inter-kernel throughput gain appears on A100 but not on GH200 and GB200. On RTX 6000 Blackwell, it appears against FFMA but, against DF emulation, only when the issue density of the DF kernel is reduced. By the logic of Section 2.4, intrakernel overlap in the WENO-Z kernel can therefore be expected on A100, not on GH200 or GB200, and on RTX 6000 Blackwell only conditionally.

5.1

A100: overlap and its source

For WENO9, 𝑠 32 = 115.84, 𝑠 64 = 60.59, and 𝑠 pure = 125.49 𝜇s/iteration, giving 𝐸 = 0.84, and 𝐸 = 0.89 for WENO5, at 𝑠 64 /𝑠 32 = 0.52. These values and the Nsight Compute counters are consistent with most of the shorter FP64 polynomial computation overlapping the FP32 weight computation: in the mixed kernel, FP64 pipe utilization falls from 98–99% standalone to 48–49%, while the FP32 FMA pipe utilization falls only from 98% to 91–93%. As expected from Section 2.3, the FP32 stream on A100 leaves issue slots free (41–46%, Table 6), which is the precondition for the overlap. Whether subsequent instructions from the same thread or other warps fill these slots cannot be determined from this measurement alone; we therefore ran the forced-serialization control in Section 3.4 (Table 5). 5

Jun Hatayama, Daichi Mukunoki, Yusuke Onishi, Atsushi Sekimoto, and Toshinori Kouchi

Table 4: Run time per precision placement on the four GPUs, speedup over FP64, and the overlap indicator. 𝐸 is measured with the compute-repeat configuration of Section 3.4, at WENO5 and WENO9 only. Each parenthesis is that row’s FP64 cell over the cell’s 𝑡 min . †: 𝑠 64 /𝑠 32 = 13.8 is extremely asymmetric and 𝐸 is not evidence of overlap. DF emulation on A100, GH200, and GB200 is out of scope. 𝑡 min [𝜇s] (speedup over FP64 in parentheses) GPU

WENO order

FP64

W32P64

DF

𝐸

A100 A100 A100

5 7 9

606.2 739.3 1078.3

199.7 (3.04) 270.3 (2.73) 392.2 (2.75)

– – –

0.89 – 0.84

GH200 GH200 GH200

5 7 9

264.2 269.7 431.3

95.1 (2.78) 138.9 (1.94) 202.8 (2.13)

– 0.12 – – – −0.02

GB200 GB200 GB200

5 7 9

220.3 233.1 370.4

76.4 (2.88) 120.5 (1.93) 174.2 (2.13)

– 0.03 – – – −0.13

RTX 6000 RTX 6000 RTX 6000

5 7 9

3303.1 557.5 (5.92) 4563.6 1021.3 (4.47) 6773.9 1415.6 (4.79)

838.0 (3.94) 1.00† 1629.3 (2.80) – 2694.5 (2.51) −0.45

Table 5: Forced-serialization control on A100, 𝑠 in 𝜇s/iteration, with 𝑠 null = 𝑠 32 + 𝑠 64 the zero-overlap reference. 𝐸 nat comes from the natural version 𝑠 nat at fixed occupancy and differs from 𝐸 in Table 4 because 𝑠 nat is an occupancy-pinned copy of the kernel (0.85 against 0.84 at WENO9). “dependent” joins the two streams through the blending once per reconstruction, so 𝑠 null and 𝐸 nat are undefined for it. 𝑠 ser and 𝑠 ind match in static and executed instruction count, occupancy, and busy time of both pipes. WENO order Placement WENO5 WENO5 WENO9 WENO9

𝑠 nat

𝑠 ser

𝑠 ind

independent 42.19 44.16 44.05 60.71 0.89 dependent 46.41 47.83 47.71 – – independent 124.80 124.23 125.82 176.48 0.85 dependent 129.81 131.95 132.58 – –

The natural version is 29–31% faster than the zero-overlap reference, 𝑠 null = 𝑠 32 + 𝑠 64 , which is consistent with temporal overlap. Forced serialization, however, recovers none of that gap: the fraction returned toward 𝑠 null , (𝑠 ser − 𝑠 ind )/(𝑠 null − 𝑠 ind ), is only +0.01 for WENO5 and −0.03 for WENO9, and the serialization changes the run time by +0.24% for WENO5 and −1.27% for WENO9. This suggests that TLP sustains the overlap on A100, but does not rule out intra-thread ILP. We took this measurement at 9.8 (WENO5) and 5.9 (WENO9) resident warps per scheduler.

5.2

𝑠 null 𝐸 nat 𝑠 ser /𝑠 ind − 1 [%] Outcome +0.24 +0.24 −1.27 −0.48

serial slower within noise serial faster serial faster

stream leaves GH200 almost no free issue slots, and FP64 instructions therefore extend the critical path instead of filling gaps. The FP64 retention remained low in Section 4.1, even as more FP64 warps were added, supporting the same interpretation. GB200, which has the lane counts of GH200 (Section 2.3), gives the same result: 𝐸 = 0.03 for WENO5 and −0.13 for WENO9, at 𝑠 64 /𝑠 32 = 0.43–0.44. For WENO9, 𝑠 pure = 52.02 𝜇s/iteration is 3.9% above 𝑠 32 + 𝑠 64 = 50.08 𝜇s/iteration, and the FP64 pipe utilization of the mixed kernel is 28.0%, against 28.8% for fully serial execution. The loss of overlap on GH200 is therefore not specific to Hopper: it occurs on a later generation with the same lane counts, as does the low FP64 retention of Section 4.1. TMA-based asynchronous stencil-tile transfer, in contrast, was effective on GH200: it was 4.0%, 1.5%, and 3.3% faster than the nontiled version for WENO5, WENO7, and WENO9, and 10.1%, 5.0%, and 4.3% faster than filling shared memory manually in each thread; on GB200 the corresponding gains are 3.8%, 4.7%, and 5.2%, and 10.2%, 10.7%, and 12.6%.

GH200 and GB200: loss of overlap

The same compute-repeat configuration gives the opposite result on GH200. 𝐸 is 0.12 for WENO5 and −0.02 for WENO9, at 𝑠 64 /𝑠 32 = 0.44: the FP32 weight computation and the FP64 polynomial computation execute almost serially. For WENO9, 𝑠 pure = 59.87 𝜇s/iteration against 𝑠 32 + 𝑠 64 = 59.58 𝜇s/iteration, which agree within 0.5%, and Nsight Compute reports an FP64 pipe utilization of 29%, against 29.1% expected for fully serial execution. Increasing the FP64-side work did not restore simultaneous issue. These observations are consistent with the expectation of Section 2.3: an FP32 instruction 6

Concurrent Use of FP32 and FP64 Computing Units on GPUs and Mixed-Precision Placement in WENO Reconstruction

Table 6: Warp-scheduler counters for the weight computation alone (FP32), the polynomial computation alone (FP64), and both in one kernel (mixed), in the compute-repeat configuration. Issue-slot utilization is the fraction of cycles on which a sub-partition’s scheduler issued. In the mixed kernel, “eligible” is warps eligible per scheduler per active cycle and “not sel.” warps eligible but not selected per issue-active cycle. issue-slot utilization [%]

mixed kernel

GPU

WENO

FP32

FP64

mixed

eligible

not sel.

A100 A100

5 9

58.8 54.4

28.7 26.3

67.7 62.1

4.56 2.75

5.71 3.41

GH200 GH200

5 9

88.8 87.7

72.4 69.9

84.9 81.3

6.98 4.22

7.18 4.17

GB200 GB200

5 9

88.1 86.2

71.6 68.8

82.9 79.3

6.82 4.68

7.18 4.87

RTX 6000 RTX 6000

5 9

91.6 88.1

6.8 6.8

11.9 11.2

0.25 0.23

1.11 1.06

5.3

mixed kernel. Thus, a shortage of eligible warps cannot explain the difference between these GPUs. These counters are consistent with the reading of Section 2.3 but do not establish it. RTX 6000 Blackwell shows their limit: its schedulers issue on only 11–12% of cycles, so issue bandwidth does not limit that GPU, and its mixed kernel has almost no eligible warps (0.23–0.25 per scheduler). These counters do not show what decides whether one pipe can use a slot the other leaves, and we do not claim to have identified it.

6

The speedups of Section 5 do not decide whether a placement is usable. Each placement was therefore embedded in an Euler solver to test whether it preserves the design order of the scheme, the order of accuracy it is constructed to attain on smooth solutions: 5, 7, and 9 for WENO5, WENO7, and WENO9. Against FP64 we compare DF emulation (labeled DF in Table 4 and in the figures), W32P64, and all-FP32. We tested two smooth fields, an entropy wave and an Euler vortex, because the computation of the smoothness indicator is sensitive to reduced precision, and our objective is to assess this sensitivity. Accordingly, we make no claims about shock-capturing performance. Both fields are governed by the two-dimensional compressible Euler equations in conservative form,

RTX 6000 Blackwell: moving off the bottleneck pipe

On RTX 6000 Blackwell, the native FP64 throughput is about 1/67.3 of FP32, and the FP64 WENO-Z kernel is bound by the FP64 pipe (FP64 pipe utilization 86.4–86.9%, issue-slot utilization 7.1–8.2%). Applying the DF emulation of Section 2.2 to the whole reconstruction gives speedups of 3.94×, 2.80×, and 2.51× over the FP64 baseline for WENO5, WENO7, and WENO9. A placement that converts only the weight computation to DF emulation and leaves the polynomial computation and the blending in native FP64 achieves a speedup of only 2.39–3.05×. Although most FP64 instructions are in the weight computation, converting it alone leaves the remaining FP64 instructions on the bottleneck pipe, and among the DF placements measured, the speedup is determined by moving these instructions to the FP32 pipe as well. This GPU shows even more clearly that free issue capacity is not a sufficient condition for overlap. For W32P64, 𝐸 = −0.45 for WENO9: no overlap appears, even though the run time of the mixed kernel is 7.8–8.9× the issue-slot demand of the two streams combined, so issue bandwidth does not limit it. At WENO9, adding the FP32 instructions raises the register count from 58 to 72 and lowers occupancy from 65.5% to 57.5%; register usage and the number of eligible warps may limit the overlap. The inter-kernel throughput gain of Section 4.2 at reduced issue density does not carry over to one kernel. For DF emulation, we report only the 2.51–3.94× speedup and make no claim about temporal overlap within it.

5.4

Numerical Accuracy at the Solver Level

𝜕Q 𝜕E 𝜕F + + = 0, 𝜕𝑡 𝜕𝑥 𝜕𝑦

(4)

with Q = (𝜌, 𝜌𝑢, 𝜌𝑣, 𝑒) ⊤ , E = (𝜌𝑢, 𝜌𝑢 2 + 𝑝, 𝜌𝑢𝑣, (𝑒 + 𝑝)𝑢) ⊤ , and F = (𝜌𝑣, 𝜌𝑢𝑣, 𝜌𝑣 2 + 𝑝, (𝑒 + 𝑝)𝑣) ⊤ . Here 𝜌 is the density, (𝑢, 𝑣) the velocity, 𝑝 the pressure, and 𝑒 = 𝑝/(𝛾 − 1) + 𝜌 (𝑢 2 + 𝑣 2 )/2 the total energy per unit volume of an ideal gas with 𝛾 = 1.4; for the one-dimensional field, 𝑣 = 0 and nothing depends on 𝑦. The onedimensional field is computed by a one-dimensional solver and the vortex by a two-dimensional one. The scheme is a nodal finitedifference WENO-Z scheme with global Lax–Friedrichs flux splitting, where the WENO reconstruction is performed independently in each coordinate direction using one-dimensional operators. Both fields are non-dimensionalized with the domain length 𝐿, the background density 𝜌 0 , and the speed of sound 𝑐 0 (pressure by 𝜌 0𝑐 02 , time by 𝐿/𝑐 0 ) and are carried by a uniform free stream at the Mach number 𝑀0 = 0.05. Time integration uses the classical fourthorder Runge–Kutta method with Δ𝑡 = 0.03 ℎ/𝑀0 , where ℎ = 1/𝑛 is the spacing of a grid with 𝑛 points in each direction in which the field varies. This is a convective CFL number 𝑢Δ𝑡/ℎ of 0.03 and an acoustic CFL number (𝑢 + 𝑐)Δ𝑡/ℎ of 0.63, with 𝑐 the speed of sound. Both fields are periodic and are advected over exactly one period, 𝑡 𝑓 = 1/𝑀0 , after which the exact solution is the initial field. Errors are reported as 𝐿2 (𝜌 − 𝜌 exact )/𝜌 0 , the root-mean-square nodal density error.

Warp-scheduler counters

6.1

Table 6 tests the lane-count expectation of Section 2.3 directly. The FP32 stream alone leaves 41–46% of the issue slots free on A100, 11–12% on GH200, and 12–14% on GB200, the order the lane counts suggest. Adding the FP64 stream raises the issue rate on A100 and lowers it on GH200 and GB200. The GH200 and GB200 mixed kernels also have more eligible warps per scheduler than the A100

Design order on a one-dimensional entropy wave

The first test case is a one-dimensional entropy wave. On the periodic domain 𝑥 ∈ [0, 1) the initial condition is   𝜌 (𝑥, 0) = 1 + 𝐴 sin(2𝜋𝑥), 𝐴 = 0.2, (𝑢, 𝑝) = 𝑀0, 𝛾1 . (5) 𝜌0 7

Jun Hatayama, Daichi Mukunoki, Yusuke Onishi, Atsushi Sekimoto, and Toshinori Kouchi

/0

test conditions, W32P64, the placement with the largest speedup on RTX 6000 Blackwell in Section 5, does not preserve the design order at WENO7 and WENO9. Of the three reduced-precision placements tested, the only one that stays close to the design order there on these grids is DF emulation.

FP64, WENO9, n = 64

1.2 1.0 0.8 0.0

0.2

0.4

x/L

0.6

0.8

6.2

1.0

We next test WENO5 using a two-dimensional field, an advected isentropic vortex; WENO7/9 are not evaluated, because the FP32weight placements already lose the design order there in a onedimensional problem. On the periodic domain [0, 1) 2 a vortex of core radius 𝑅𝑐 = 0.05 is superimposed on the uniform flow (𝑢, 𝑣, 𝑝) = (𝑀0, 0, 1/𝛾),    2 𝑅2 𝑇 𝑢𝜃 𝑅 𝑅 = 𝛽 v 𝑀0 exp − 2 , = 1 − 12 (𝛾 − 1)𝑀02 𝛽 v2 exp − 2 , 𝑐0 𝑅𝑐 𝑇0 2𝑅𝑐 𝑅𝑐 (6) where 𝑅 is the distance from the vortex center, 𝑢𝜃 the tangential velocity, 𝑇0 the background temperature, and 𝛽 v the vortex strength, unrelated to the 𝛽𝑘 of Eq. (2). Homentropy gives 𝜌/𝜌 0 = (𝑇 /𝑇0 ) 1/(𝛾 −1) and the equation of state 𝛾𝑝 = (𝜌/𝜌 0 ) (𝑇 /𝑇0 ). The peak density deficit is 𝛿𝜌/𝜌 0 ≃ 𝑀02 𝛽 v2 /2 = 5.0 × 10−3 at the 𝛽 v = 2 used here (Figure 1(b)), and the peak vortex Mach number is about 0.06, so the field remains smooth. The grids are 𝑛 = 64, 128, 256, i.e. 3.2, 6.4, and 12.8 points per 𝑅𝑐 . Figure 3 shows the results. At 𝑛 = 64 the four placements agree within 1.3% at 1.730–1.752 × 10−4 . The local order, the slope of log 𝐿2 against log ℎ between successive grids, is 5.88 for FP64 and 5.52 for DF emulation from 𝑛 = 128 to 256, near the design order 5; phase error remaining at 𝑛 = 128 and the time step limit this estimate to about one digit. At this time step the W32P64 errors are 1.91 and 1.26 times those of FP64 at 𝑛 = 128 and 256. We quote no local order for W32P64, because its 𝑛 = 256 error depends on the time step. The all-FP32 reconstruction stops converging at a local order of 0.29; its 𝑛 = 256 error of 8.8 × 10−5 is 1.8% of the peak density deficit, 131 times that of FP64 and 104 times that of W32P64. At WENO5 and this time step, then, moving only the weights to FP32 costs a factor of 1.26 in the 𝑛 = 256 error, whereas moving the whole reconstruction to FP32 costs a factor of 131. This is consistent with the observation of Field et al. [5] that the weights need less precision than the polynomial values.

(a) entropy wave 1.00

1.000 0.998

/0

y/L

0.75 0.50 0.25

Grid convergence on Euler vortex advection

0.996

0.00 0.00 0.25 0.50 0.75 1.00 x/L

(b) Euler vortex

Figure 1: The two smooth flow fields on which the placements are compared, as the density normalized by the background density. (a) The entropy wave: the line is Eq. (5), the initial field and, by periodicity, the exact solution after one period; the circles are the computed FP64 solution after one period at WENO9, 𝑛 = 64. (b) The Euler vortex of Eq. (6) after one period, computed in FP64 at WENO5, 𝑛 = 256, 𝛽 v = 2. The velocity and the pressure are uniform, so the density perturbation is advected without change of shape (Figure 1(a)). The grids are 𝑛 = 16, 32, 64; on finer grids the truncation error falls below the noise floor of DF emulation. Figure 2 shows the results, with the fitted order defined as the least-squares slope of log 𝐿2 against log ℎ over 𝑛 = 16, 32, 64. FP64 comes close to the design orders on these grids, with 4.98/6.98/8.96 at WENO5/7/9, and DF emulation follows at 4.98/6.98/8.95. At 𝑛 = 64 the errors of the two are identical at WENO5, agree to five significant digits at WENO7, and differ by a factor of 1.01 at WENO9 (1.690 × 10−11 against 1.669 × 10−11 ). We attribute this difference to the noise floor of the roughly 48-bit DF significand, not to a loss of design order. The FP32-weight placements, in contrast, keep the design order only at WENO5, with a fitted order of 4.91 for both and 𝑛 = 64 errors within a factor of 1.10 of FP64. At WENO7 the fitted orders of W32P64 and all-FP32 fall to 1.78 and 1.64, and at WENO9 to −2.94 and −3.26. At both orders the error increases from 𝑛 = 32 to 64; at WENO9, that of W32P64 grows from 4.97 × 10−7 to 2.52 × 10−4 , seven orders of magnitude above the 1.67 × 10−11 of FP64. This loss of order may be attributed to the smoothness indicators. The 𝛽𝑘 of WENO7/9 are quadratic forms with many cross terms, which can cancel more strongly as refinement makes the field locally smoother. In FP32, 𝛽𝑘 may lose its significant digits while the truncation error is still falling, so 𝜏/(𝛽𝑘 + 𝜖) in Eq. (2) can be dominated by round-off. Consequently, in our implementation and

6.3

The solver time step on the roofline

Figure 4 places the four stages of one Runge–Kutta step of the onedimensional solver on the roofline of each GPU. Each stage moves at least four times the largest L2 cache, and the measured DRAM traffic is within 0.4% of the 408 bytes per point and step that the stages must move. On A100, GH200, and GB200, W32P64 achieves a speedup of 1.61–2.36× over FP64 at WENO5 and WENO9, whereas DF emulation makes it 4.0–12.3× slower. On RTX 6000 Blackwell, DF emulation achieves a speedup of 2.27–2.63× and W32P64 a speedup of 3.71–4.18×, against 4.47–5.92× for the W32P64 kernel alone (Table 4). The WENO step is not necessarily memory-bound: whether it is depends on the scheme order and the placement. We compare the time of each step with the longest of its FP64 FLOP at the FP64 peak, its FP32 FLOP at the FP32 peak, and its DRAM bytes at the peak 8

Concurrent Use of FP32 and FP64 Computing Units on GPUs and Mixed-Precision Placement in WENO Reconstruction

exact

)/

ρ

0

10−3

L2 (

ρ−ρ

10−6 10

h

5

h

−9

7

h

10−12

1/64

1/32

1/16

h = 1/n

1/64

1/32

1/16

1/64

h = 1/n

(a) WENO5 FP64

(b) WENO7 DF W32P64

9

1/32

1/16

h = 1/n

(c) WENO9 FP32

Figure 2: Grid convergence on the one-dimensional entropy wave at (a) WENO5, (b) WENO7, and (c) WENO9: the normalized density error 𝐿2 (𝜌 − 𝜌 exact )/𝜌 0 after one period against the grid spacing ℎ = 1/𝑛. Dotted lines show the slope of the design order. A100 shows both an inter-kernel throughput gain and intrakernel overlap; GH200 and GB200 show neither. On RTX 6000 Blackwell, the gain appears between kernels against FFMA and, against DF emulation, once the issue density of the DF kernel is reduced, but no overlap appears within a single kernel with FP32 weights and FP64 polynomials (W32P64). On the GPUs evaluated here, intrakernel overlap appeared only with an inter-kernel throughput gain, but the gain does not imply it. On RTX 6000 Blackwell, even with free issue bandwidth, the added FP32 instructions raise register usage and lower occupancy, which may leave too few eligible warps to sustain the overlap. Moreover, the issue-slot counters are consistent with the lane counts on A100, GH200, and GB200 but not on RTX 6000 Blackwell, so whether free slots can actually be used must be measured on each GPU. The main findings are as follows.

10−5

L2 (

ρ−ρ

exact

)/

ρ

0

10−4

10−6 design 1/256

h

5

1/128

1/64

h = 1/n

FP64

DF

W32P64

FP32

(1) In the microbenchmark, A100 sustains about 90% of the standalone throughput on both the FP32 and FP64 sides with pure kernels, reducing wall time by 44.0%, and gains 3.8% over native FP64 alone when running concurrently with a DF emulation kernel. RTX 6000 Blackwell gains 6.7% over DF emulation alone when the issue density of the DF kernel is reduced, whereas on GH200 the combined throughput falls at least 40.2% below native FP64 alone in every configuration, and on GB200 at least 47.6%. (2) Inside the WENO-Z kernel, A100 gives an overlap indicator of 𝐸 = 0.89 and 0.84 for WENO5 and WENO9, respectively, against 𝐸 ≤ 0.12 on GH200 and GB200. A forcedserialization control barely changes the run time, which is consistent with TLP but does not rule out intra-thread ILP. (3) On RTX 6000 Blackwell, the native FP64 pipe is the bottleneck, and applying DF emulation to the whole reconstruction reaches a 2.51–3.94× speedup over native FP64. Converting only the weight computation to DF emulation achieves only 2.39–3.05× speedup; the remaining gain comes from also moving the FP64 instructions left in the polynomial computation onto the FP32 pipe. (4) The fastest placement does not always preserve the design order of accuracy. For the one-dimensional entropy wave, W32P64 retains 4.91 for WENO5 but drops to 1.78 for

Figure 3: Grid convergence on Euler vortex advection at WENO5, 𝛽 v = 2: the normalized density error 𝐿2 (𝜌 − 𝜌 exact )/𝜌 0 after one period against ℎ = 1/𝑛. The dotted line is the slope of the design order ℎ 5 .

bandwidth. Among FP64, W32P64, and DF emulation, only W32P64 at WENO5 on A100 and GH200 runs within twice its DRAM time (1.43–1.52 times). FP64 runs at 1.42–1.87 times its FP64-pipe time on A100 and RTX 6000 Blackwell, and DF emulation at 2.10–4.67 times its FP32-pipe time on every GPU. Register count, and hence occupancy, differs by placement and was not equalized: on RTX 6000 Blackwell, FP64 uses 56–76 registers at 50–66% occupancy and W32P64 37–58 at 66–83%.

7

Conclusion

We investigated whether idle computing units can be utilized by concurrently executing instructions of different precisions on the same SM and, if so, whether this technique can be applied to a practical CFD kernel. 9

Jun Hatayama, Daichi Mukunoki, Yusuke Onishi, Atsushi Sekimoto, and Toshinori Kouchi

performance [GFLOP/s]

FP64 105

DF

W32P64

FP32

104

FP64

FP32

WENO5

WENO9

FP32

FP32

FP64

FP64

FP32

FP64

103 100

102

104

100

102

104

100

102

104

100

102

104

arithmetic intensity [FLOP/B] arithmetic intensity [FLOP/B] arithmetic intensity [FLOP/B] arithmetic intensity [FLOP/B] A100 GH200 GB200 RTX 6000 Blackwell Figure 4: The four stages of one classical Runge–Kutta step of the one-dimensional solver at WENO5 (circles) and WENO9 (squares), 𝑛 = 224 , on the roofline of each GPU; colors distinguish the placements. The intensity is FP64 plus FP32 FLOP (an FMA as two) per DRAM byte, and the ceilings are published peaks (solid: DRAM, FP64; dashed: FP32), the FP64 one of RTX 6000 Blackwell being FP32 over 64.

References

WENO7 and −2.94 for WENO9; for WENO9 with 𝑛 = 64, its error is seven orders of magnitude larger than that of FP64. DF emulation stays close to the design order in all three cases on the grids tested. For the two-dimensional Euler vortex with WENO5, at the time step used, the W32P64 error is 1.26× that of FP64 at 𝑛 = 256, while the all-FP32 placement stops converging (local order 0.29, with an error 104× that of W32P64).

[1] Rafael Borges, Monique Carmona, Bruno Costa, and Wai Sun Don. 2008. An improved weighted essentially non-oscillatory scheme for hyperbolic conservation laws. J. Comput. Phys. 227, 6 (2008), 3191–3211. [2] Marcos Castro, Bruno Costa, and Wai Sun Don. 2011. High order weighted essentially non-oscillatory WENO-Z schemes for hyperbolic conservation laws. J. Comput. Phys. 230, 5 (2011), 1766–1792. [3] Yanxiang Chen, Pablo de Oliveira Castro, Paolo Bientinesi, Niclas Jansson, and Roman Iakymchuk. 2026. Enabling mixed-precision in spectral element codes. Future Generation Computer Systems 174 (2026), 107990. [4] Theodorus Jozef Dekker. 1971. A floating-point technique for extending the available precision. Numer. Math. 18, 3 (1971), 224–242. [5] Scott E Field, Sigal Gottlieb, Zachary J Grant, Leah F Isherwood, and Gaurav Khanna. 2023. A GPU-accelerated mixed-precision WENO method for extremal black hole and gravitational wave physics computations. Communications on Applied Mathematics and Computation 5, 1 (2023), 97–115. [6] Yozo Hida, Xiaoye S Li, and David H Bailey. 2001. Algorithms for quad-double precision floating point arithmetic. In Proceedings 15th IEEE Symposium on Computer Arithmetic. ARITH-15 2001. IEEE, 155–162. [7] Guang-Shan Jiang and Chi-Wang Shu. 1996. Efficient implementation of weighted ENO schemes. J. Comput. Phys. 126, 1 (1996), 202–228. [8] Martin Karp, Ronith Stanly, Timofey Mukha, Luca Galimberti, Siavash Toosi, Hang Song, Lisandro Dalcin, Saleh Rezaeiravesh, Niclas Jansson, Stefano Markidis, et al. 2025. Effects of lower floating-point precision on scale-resolving numerical simulations of turbulence. J. Comput. Phys. (2025), 114600. [9] Aditya Kashi, Hao Lu, Wesley Brewer, David Rogers, Michael Matheson, Mallikarjun Shankar, and Feiyi Wang. 2026. Mixed-precision numerics in scientific applications: survey and perspectives. The Journal of Supercomputing 82, 5 (2026), 287. [10] Xu-Dong Liu, Stanley Osher, and Tony Chan. 1994. Weighted essentially nonoscillatory schemes. J. Comput. Phys. 115, 1 (1994), 200–212. [11] NVIDIA. 2020. NVIDIA A100 Tensor Core GPU Architecture. White Paper. NVIDIA Corporation. https://images.nvidia.com/aem-dam/en-zz/Solutions/data-center/ nvidia-ampere-architecture-whitepaper.pdf [12] NVIDIA. 2023. NVIDIA H100 Tensor Core GPU Architecture. White Paper. NVIDIA Corporation. https://resources.nvidia.com/en-us-hopper-architecture/nvidiah100-tensor-c [13] NVIDIA Corporation. 2025. CUDA C++ Programming Guide. NVIDIA Corporation. https://docs.nvidia.com/cuda/archive/13.0.0/cuda-c-programmingguide/index.html Version 13.0. [14] Bálint Siklósi, Pushpender K Sharma, David J Lusher, István Z Reguly, and Neil D Sandham. 2025. Reduced and mixed precision turbulent flow simulations using explicit finite difference schemes. Future Generation Computer Systems (2025), 108111.

To realize the performance benefits of mixed-precision computation and DF emulation in CFD solvers, effective redistribution of computational resources is essential. Although NVIDIA GPUs generally offer limited FP64 throughput, increasing the instruction issue bandwidth could enable the concurrent execution of instructions with different precisions, thereby utilizing otherwise idle execution units and improving overall performance. Future work will investigate the impact of W32P64 on high-order turbulence statistics and shock-containing flows, as well as extend the proposed resource redistribution strategy to FP16 and integer execution units. Ultimately, our goal is to develop novel CFD algorithms that fully exploit the available computational and memory resources of each GPU architecture.

Acknowledgments This research was supported by the Collaborative Research Project on Computer Science with High-Performance Computing at Nagoya University and the Core Facility of Okayama University (CFPOU Engineering 001), which provided computational resources. In addition, this work received financial support from JST SPRING, Japan (Grant Number JPMJSP2126), and the Public Interest Incorporated Foundation Ohmoto Ikueikai. During this work, the authors used Claude Code, Claude, Codex, and ChatGPT for code development and manuscript writing. 10

Related documents

• Objective Coefficient Rounding and Almost Symmetries in Binary Programs#984101
• Time-to-Digital Conversion in the Lattice CertusPro FPGA for the CAPE Instrument#623773
• The Price of Precision: A Calibrated Cost Model for Mixed-Precision LLM Serving#1326879
• Comparison of SITA-Standard Compared to SITA-Fast Visual Fields.#290886
• Improvement of the energy resolution via an optimized digital signal processing in GERDA Phase I#611970