H2SGEMM: FP32 GEMM Emulation on FP16 Hardware

• H2SGEMM is a high-performance algorithm that decomposes each FP32 operand into two FP16 components to emulate FP32 GEMM on FP16 hardware.
• It employs a two-component decomposition and three-term reconstruction with scaling and term-wise accumulation to recover up to 22 bits of mantissa accuracy.
• Optimized for Huawei Ascend 910A, H2SGEMM leverages cache-aware blocking and double-buffered pipelining to achieve a throughput of 65.3 TFLOP/s.

H2SGEMM is a high-performance algorithm for emulating FP32 general matrix-matrix multiplication (GEMM) using only FP16 computation units on AI accelerators that provide high-throughput FP16 matrix engines but do not provide equally capable native FP32 matrix units. The representative platform is the Huawei Ascend 910A NPU, whose matrix compute capability is concentrated in FP16 “cube” units. H2SGEMM addresses the mismatch between FP16 throughput and FP32 accuracy requirements by decomposing each FP32 operand into two FP16 values, scaling residuals to avoid underflow, reconstructing the product through three FP16 GEMMs, and combining the result with cache-aware blocking and double-buffered pipelining. The method is reported to preserve up to 22 bits of mantissa accuracy and to achieve 65.3 TFLOP/s, or 77% of the theoretical FP32-equivalent peak, on Ascend 910A (Xue et al., 31 Jul 2025).

1. Computational setting and design objective

H2SGEMM targets workloads that still require the numerical quality of FP32 GEMM on hardware whose high-performance matrix path is FP16. The central operation is

$C = AB,$

with FP32 input matrices $A$ and $B$, but the computation is carried out using FP16 matrix engines plus FP32-side orchestration and accumulation. The challenge is that an FP16 multiply has only about 11 bits of significand precision including the hidden bit, whereas FP32 has 24 bits. A naive cast of FP32 inputs to FP16 before GEMM collapses precision dramatically and can also expose the computation to FP16 underflow and overflow behavior (Xue et al., 31 Jul 2025).

The target hardware in the reported implementation is the Ascend 910A NPU. Its matrix compute capability is concentrated in FP16 cube units, and the system-level contribution of H2SGEMM is to make those FP16 units serve as an FP32-equivalent GEMM path. In the terminology used for the method, “FP32-equivalent” does not mean that the hardware executes native IEEE FP32 multiply-adds throughout. It means that the algorithm recovers FP32-level numerical accuracy well enough to match native SGEMM in most cases and, under some low-exponent conditions, exceed its stability (Xue et al., 31 Jul 2025).

2. Two-component FP16 decomposition and three-term reconstruction

The numerical core of H2SGEMM is a two-component decomposition of each FP32 scalar into a high FP16 part and a residual FP16 part. The method represents each FP32 scalar $V$ as

$V \approx V_{\text{high}} + V_{\text{low}}.$

The decomposition is based on round-to-nearest-even conversion from FP32 to FP16. The high part keeps the leading FP16-significant bits, and the low part captures the residual left after subtracting the rounded high part (Xue et al., 31 Jul 2025).

At matrix level, the decomposition is written operationally as

$$A_{\text{half}} = \text{to\_half}(A_{\text{single}}), \qquad R_{A,\text{half}} = \text{to\_half}\big((A_{\text{single}} - \text{to\_single}(A_{\text{half}})) \times sf\big),$$

and analogously for $B$, with $sf = 2^{sb}$ a tunable scaling factor. This gives the approximations

$$A_{\text{single}} \approx A_{\text{half}} + \frac{R_{A,\text{half}}}{sf}, \qquad B_{\text{single}} \approx B_{\text{half}} + \frac{R_{B,\text{half}}}{sf}.$$

Substituting these decompositions into the GEMM yields the three-term reconstruction

$$C \approx A_{\text{half}} B_{\text{half}} + \frac{R_{A,\text{half}} B_{\text{half}}}{sf} + \frac{A_{\text{half}} R_{B,\text{half}}}{sf}.$$

The second-order term

$A$0

is neglected because it is second-order small. This omission is what keeps the compute cost to roughly three FP16 GEMMs. Algorithmically, H2SGEMM is therefore a 2-split operand decomposition leading to a 3-term GEMM reconstruction (Xue et al., 31 Jul 2025).

3. Scaling strategy, underflow control, and mantissa recovery

Residual scaling is required because the low part can be too small to be representable in FP16, especially when the original FP32 input has a small exponent. Without scaling, the residual can underflow or become subnormal, and the decomposition cannot preserve the intended extra mantissa bits. H2SGEMM addresses this by multiplying the residual by $A$1 before FP16 storage and compensating for that scaling in the reconstructed GEMM (Xue et al., 31 Jul 2025).

The method extracts two explicit engineering rules from its error analysis.

Under RN rounding, when the absolute value of an FP32 number is below $A$2, the low part must be scaled to preserve 22 bits. If the value is below $A$3, decomposition fails without sufficient scaling.

Under RN rounding, if large FP32 values are present, the scaling factor should not exceed $A$4 to avoid overflow in the lower component.

The admissible range for the scaling exponent $A$5 is given as

$A$6

When the data range is unknown, the recommended conservative choice is

$A$7

The stated numerical goal is preservation of up to 22 bits of mantissa accuracy. The argument is that the FP16 high part captures roughly 11 significant bits including the hidden bit, while the scaled residual preserves much of the remaining FP32 mantissa. The paper also states an important boundary condition: if the FP32 number is below the FP16 subnormal threshold $A$8, then both high and low parts would need scaling, and that case is beyond the current scope (Xue et al., 31 Jul 2025).

4. Accumulation order and numerical stability

H2SGEMM distinguishes two accumulation orders for the expanded product. In element-wise evaluation, all three contributions are formed and summed in the ordinary GEMM reduction flow. In term-wise evaluation, the algorithm first computes the three component matrices independently,

$A$9

and then sums these matrices afterward (Xue et al., 31 Jul 2025).

The term-wise accumulation scheme is a central numerical feature. The stated rationale is that $B$0 and $B$1 are much smaller than $B$2 and $B$3, so mixing all contributions too early can round away the corrections. By computing the correction GEMMs separately, H2SGEMM preserves their significance longer. The paper describes the benefit as effectively summing smaller-magnitude terms first, followed by larger terms, reducing rounding loss compared with adding tiny values into already-large accumulators (Xue et al., 31 Jul 2025).

This matters especially in low-exponent regimes, where underflow and significance loss are more severe. In that regime, the term-wise scheme is reported to improve numerical stability over conventional FP32 GEMM. The paper states that, with $B$4, term-wise H2SGEMM can surpass SGEMM for small exponents and that as $B$5 increases, term-wise H2SGEMM scales better than SGEMM and other H2SGEMM variants in relative error (Xue et al., 31 Jul 2025).

5. Hardware mapping, blocking, and pipelined execution on Ascend 910A

The reported implementation maps the three GEMM terms onto the Ascend 910A NPU, which contains 32 AI cores, high-throughput FP16 cube units, and a memory hierarchy involving L1, L0A, L0B, L0C, and UB. The hardware provides about 256 TFLOP/s FP16 peak, and because H2SGEMM uses three major FP16 GEMMs per logical FP32 GEMM, the theoretical FP32-equivalent peak is defined as

$B$6

Decomposition and re-scaling are done by vector operations and format conversion, while the heavy matrix products use the cube unit. Intermediate and final accumulation is done in FP32 (Xue et al., 31 Jul 2025).

The performance-critical optimization is cache-aware blocking. The stated bottleneck is not L1-to-L0 bandwidth, which is high, but main-memory to L1 traffic. The blocking strategy therefore aims to maximize L1 residency and reuse of $B$7 blocks while reserving L1 space for double buffering of $B$8 blocks. The hardware constraints on block sizes are

$B$9

$V$0

$V$1

$V$2

The factor of 6 reflects the storage pressure from multiple terms and buffers in the expanded computation (Xue et al., 31 Jul 2025).

The second major optimization is a double-buffered pipeline. In the single-buffered case, compute must wait for global memory $V$3 UB, conversion, UB $V$4 L1, L1 $V$5 L0, and cube compute. In the double-buffered case, one buffer is used by the cube for the current tile while another buffer prefetches and converts the next tile. The paper describes four overlapped stages for $V$6: global memory to UB, conversion from FP32 to FP16, UB to L1 transfer, and movement of $V$7 and $V$8 into L0 for cube execution (Xue et al., 31 Jul 2025).

The best measured blocking configuration is

$V$9

6. Accuracy evaluation, achieved throughput, and operational limits

Accuracy is evaluated against FP64 ground truth using

$V \approx V_{\text{high}} + V_{\text{low}}.$0

The test methodology uses uniform random input matrices with controlled exponent ranges, including symmetric sampling from

$V \approx V_{\text{high}} + V_{\text{low}}.$1

and nonnegative sampling from

$V \approx V_{\text{high}} + V_{\text{low}}.$2

with particular emphasis on non-positive offset exponents (Xue et al., 31 Jul 2025).

The reported numerical findings are explicit. Vanilla HGEMM has relative error around $V \approx V_{\text{high}} + V_{\text{low}}.$3, which is stated to be too large for precision-sensitive applications. H2SGEMM without scaling improves accuracy but does not consistently match SGEMM near small exponents. Scaling with $V \approx V_{\text{high}} + V_{\text{low}}.$4 is insufficient. Scaling with $V \approx V_{\text{high}} + V_{\text{low}}.$5 reduces error by 1–2 orders of magnitude in the low-exponent regime and is identified as the decisive factor for precision recovery. With $V \approx V_{\text{high}} + V_{\text{low}}.$6, element-wise H2SGEMM becomes comparable to SGEMM, while term-wise H2SGEMM can surpass SGEMM for small exponents (Xue et al., 31 Jul 2025).

The headline performance result is 65.3 TFLOP/s on Ascend 910A for FP32-equivalent GEMM, corresponding to 77% of the theoretical FP32-equivalent peak. The double-buffered pipeline increases throughput from 41.7 TFLOP/s to 65.3 TFLOP/s, a gain of about 57%. The method is also reported to sustain around 60 TFLOP/s over a wide range of matrix sizes, and in some scaling experiments to slightly outperform the vendor CANN SGEMM measured on a newer Ascend 910B3 with native FP32 support (Xue et al., 31 Jul 2025).

The method also has explicit limits. Very small FP32 inputs below the FP16 subnormal threshold $V \approx V_{\text{high}} + V_{\text{low}}.$7 are beyond the current scope because both high and low parts would need scaling. The neglected term $V \approx V_{\text{high}} + V_{\text{low}}.$8 remains an approximation, not an exact reconstruction. The cache-aware blocking and buffering are tightly tuned to Ascend’s L1/L0/UB hierarchy and cube constraints. The claim is recovery of native FP32 accuracy in practice and up to 22 bits of mantissa accuracy, not guaranteed bit-for-bit equality in all cases (Xue et al., 31 Jul 2025).

7. Nomenclature and relation to adjacent GEMM literatures

In the literature, H2SGEMM refers to FP32 GEMM emulation on Ascend NPUs via two FP16 components per operand, three FP16 GEMM terms, precision-recovery scaling, and hardware-aware execution (Xue et al., 31 Jul 2025). This should be distinguished from $V \approx V_{\text{high}} + V_{\text{low}}.$9-matrix–dense-matrix multiplication. Later work on block Krylov solvers presents efficient $$A_{\text{half}} = \text{to\_half}(A_{\text{single}}), \qquad R_{A,\text{half}} = \text{to\_half}\big((A_{\text{single}} - \text{to\_single}(A_{\text{half}})) \times sf\big),$$0-matrix-vector and $$A_{\text{half}} = \text{to\_half}(A_{\text{single}}), \qquad R_{A,\text{half}} = \text{to\_half}\big((A_{\text{single}} - \text{to\_single}(A_{\text{half}})) \times sf\big),$$1-matrix-matrix multiplication that exploit BLAS level-3 routines such as GEMM, but that work does not use the name H2SGEMM explicitly and addresses hierarchical matrix formats rather than FP32 emulation on FP16 cube hardware (Christophersen, 21 Sep 2025).

It should also be distinguished from SGEMM performance-portability research on GPUs. Multi-versioning approaches for SGEMM select a small set of precompiled kernel variants to cover multiple devices and matrix shapes, addressing autotuning overfitting and runtime compilation cost. That line of work concerns dense SGEMM kernel selection and dispatch across heterogeneous GPUs rather than numerical reconstruction of FP32 GEMM on FP16-only matrix engines (Hochgraf et al., 21 Jul 2025).

Under this narrower and now specific usage, H2SGEMM is best understood as a software–hardware co-designed FP32 GEMM emulation method for accelerators whose practical dense-matrix throughput is concentrated in FP16 units. Its defining characteristics are the two-component FP16 encoding of each FP32 operand, the three-term reconstruction with a neglected second-order term, the tunable scaling strategy centered on $$A_{\text{half}} = \text{to\_half}(A_{\text{single}}), \qquad R_{A,\text{half}} = \text{to\_half}\big((A_{\text{single}} - \text{to\_single}(A_{\text{half}})) \times sf\big),$$2, the term-wise accumulation scheme, and the Ascend-specific cache and pipeline optimizations that make the reconstruction operationally competitive (Xue et al., 31 Jul 2025).