Many scientific simulations, from molecular chemistry to materials science, spend the majority of their execution time computing high-precision dense general matrix multiplications (GEMM). These operations are called DGEMM for real 64-bit floating-point numbers and ZGEMM for complex 64-bit floating-point numbers. cuBLAS is widely used by HPC and AI applications to accelerate GEMM and other Basic Linear Algebra Subprograms (BLAS) on NVIDIA GPUs. Rapid advances in lower-precision hardware have motivated exploration of a new class of algorithms that deliver full double-precision results on the faster integer Tensor Cores.
With CUDA Toolkit 13.4.2, cuBLAS no longer has to choose between accuracy and speed. cuBLAS breaks each double-precision number into integer pieces, multiplies them on the GPU’s Tensor Cores, and reassembles the exact result. The answers keep full double-precision accuracy, but the work runs on faster hardware.
Two new algorithms make this more efficient: Ozaki-II needs fewer integer multiplications than the earlier Ozaki-I approach, and the 2M method cuts a complex matrix multiplication from four real multiplications to two. This builds on the Automatic Dynamic Precision (ADP) framework introduced in Unlocking Tensor Core Performance with Floating Point Emulation in cuBLAS, which delivers higher double-precision (FP64) matrix-multiplication performance by mapping FP64 arithmetic onto integer Tensor Cores while preserving native FP64 accuracy.
cuBLAS delivers automatic acceleration of DGEMM and ZGEMM with no application code changes. Because speed matters only when accuracy is not compromised, we validate both on real-world applications: coupled cluster singles and doubles (CCSD) in quantum chemistry, and BerkeleyGW in materials science. This post explains how the algorithms work, shows benchmark and application results, and covers how to get started.
Performance highlights #
cuBLAS in CUDA Toolkit 13.4.2 delivers the following performance gains:
- cuBLAS DGEMM with Ozaki-II runs 5.5x and 3.9x faster than native FP64 on NVIDIA VR NVL72 and NVIDIA GB200, respectively. It’s also 2x faster than Ozaki-I on NVIDIA GB200.
- CCSD case study using cuBLAS DGEMM with Ozaki-II achieves 2x and 1.4x end-to-end speedup over native FP64 on NVIDIA VR NVL72 and NVIDIA GB200 in FP64. In domain-specific mixed-precision mode, end-to-end speedup over native FP64 is 5x and 3.5x on NVIDIA VR NVL72 and NVIDIA GB200 respectively.
- cuBLAS ZGEMM with 2M algorithm achieves a 5x speedup over native Complex FP64 (CFP64) on NVIDIA GB200.
- BerkeleyGW case study using cuBLAS ZGEMM with 2M algorithm achieves 6.6x and 3.9x speedup over native ZGEMM on NVIDIA VR NVL72 and NVIDIA GB200 respectively.
- cuBLAS DGEMM with Ozaki-II achieves up to 175, 147, and 37 TFLOPS on NVIDIA VR NVL72, NVIDIA GB200, NVIDIA RTX PRO 6000 respectively.
From Ozaki-I to Ozaki-II
Ozaki-I and Ozaki-II schemes share a common first step, but differ in how the big integer multiplication is performed. Both the schemes align the exponents of the input matrices and convert them to a shared-exponent, variable-length integer representation so that products accumulate exactly on integer Tensor Cores. cuBLAS uses heuristics to automatically choose between Ozaki-I and Ozaki-II for best performance. The integer width is chosen by ADP at runtime or can also be set explicitly by the user.
Ozaki-I breaks each big integer into 8-bit slices and multiplies slices of A against slices of B in a triangular pattern, so the number of Tensor Core GEMMs grows with the square of the slice count.
Ozaki-II uses a faster big-integer multiplication algorithm based on the Chinese Remainder Theorem (CRT). It reduces each quantized operand modulo a set of small pairwise-coprime integers, performs exactly one integer GEMM per modulus, and reconstructs the exact product from the residues. With Ozaki-II, the GEMM count is linear in the number of moduli, whereas Ozaki-I grows with the square of the slice count, so the advantage widens as the accuracy target rises. At the accuracy required to match FP64, Ozaki-I needs roughly twice as many Tensor Core GEMMs as Ozaki-II.
Because Ozaki-II materializes one residue pair per modulus, a naive implementation would carry a much larger workspace than Ozaki-I. However, cuBLAS already applies the memory-parsimonious tiling introduced for Ozaki-I in CTK 13.3 to Ozaki-II scheme as well, discarding each residue pair immediately after its GEMM. This results in a more manageable workspace size.
cuBLAS ADP framework is designed to provide accuracy guarantees comparable to, or better than, standard floating-point GEMM algorithms. More information on the accuracy of fixed-point methods can be found in the references. If at run-time, the number of mantissa bits required to preserve FP64-equivalent accuracy exceeds the library default or a user-supplied maximum, cuBLAS falls-back to a code path that utilizes the native FP64 MMA instructions. So, applications can safely use cuBLAS in their production deployments without worrying about numerical accuracy. Figure 1, below, shows the speedup of cuBLAS DGEMM using Ozaki-II over Ozaki-I on GB200.
DGEMM performance #
NVIDIA Rubin
Rubin introduces a new numeric format, TI16 (12-bit integer in sign-magnitude format with range of [-2047, +2047]), that is aptly suited for the acceleration of Ozaki scheme computations. cuBLAS 13.4 uses TI16 type for acceleration of Ozaki-I and Ozaki-II computations through the tcgen05.mma instruction on NVIDIA Rubin GPUs (compute capability 10.7). Figure 2, below, shows the speedup of cuBLAS DGEMM using Ozaki-II over native FP64 on VR NVL72.
DGEMM application performance – CCSD in quantum chemistry
To demonstrate the potential application performance of Ozaki-II DGEMM, we consider the use case of coupled cluster singles and doubles (CCSD). CCSD is a foundational post-DFT method in electronic structure that together with the (T) perturbative correction (not considered here) forms the “gold standard” of quantum chemistry.
Moreover, the response of CCSD to various perturbations yields the EOM-CCSD treatment of electronic excited states and the CCSD lambda treatment of property computations. CCSD is noteworthy for its very high cost of dense FP64 tensor contractions, which the Ozaki-II scheme in cuBLAS can significantly accelerate.
A small CCSD demonstration was implemented using cuBLAS in CTK 13.4.2 for the main tensor contractions, with the coupled cluster problem setup accelerated by cuEST v0.2.2. A standard FP64 single-GPU T1-transformed restricted CCSD approach was constructed in the density fitting representation. Figure 3, below, shows performance data for frozen-core CCSD/cc-pVTZ for five mid-sized members of the S22 database — these systems contain 24–28 atoms and roughly 590–630 atomic orbitals.
Ozaki-II DGEMM with ESC guaranteeing full FP64 precision causes performance to rise monotonically with generation, yielding a total speedup of 4.7x of Rubin over A100. Another compelling possibility is to use an application-specific mixed precision mode with a fixed Ozaki-II mantissa bit budget, here selected to be 32 mantissa bits. This also yields monotonic performance gains, with total speedup now at 11.5x native FP64 A100. With ESC at full FP64 quality, results are indistinguishable from native FP64 precision – errors in the CCSD correlation energy are <10<sup>-14</sup> E<sub>h</sub>.
Ozaki-II emulation provides the opportunity to use lower fixed precision with 32 mantissa bits while still providing sufficient accuracy for the vast majority of practical applications; errors in the energy are <10<sup>-10</sup> E<sub>h</sub> and the amplitude equations can be converged to a residual of 10<sup>-9</sup> or better.
On Rubin for uracil dimer frozen-core CCSD/cc-pVTZ, fixed precision with 32 mantissa bits achieves an end-to-end throughput of 173 TFLOPS; the throughput of individual GEMMs is as high as 368 TFLOPS with an average GEMM throughput of 200 TFLOPS. Similar throughputs are obtained during the solution of the CCSD lambda response equations needed for molecular property calculations (not shown in plot).
Complex GEMMs: from 4M to 3M to 2M
A complex GEMM is conventionally computed as four real GEMMs (4M) of the same size:
A = Ar + i(Ai); B = Br + i(Bi)
C = A x B = (Ar·Br − Ai·Bi) + i(Ar·Bi + Ai·Br)
The well-known 3M algorithm reduces this to three real GEMMs plus quadratic-time pre- and post-processing. In floating-point arithmetic, 3M is used cautiously because forming the sums Ar + Ai and Br + Bi, and subtracting the intermediate products during reconstruction can degrade numerical stability.
With the Ozaki-II scheme, the additions, subtractions, and products are all exact integer operations, so the choice of complex multiplication method no longer affects the accuracy of the result. Multiplying two complex numbers cannot be done in fewer than three real multiplications, but matrix products in this exact integer setting are not bound by that result: the 2M algorithm computes a complex product with only two real GEMMs instead of four. The end result is that a ZGEMM now costs roughly twice a real GEMM of the same size instead of four times. The full derivation is available in the 2M paper.
Starting with CUDA Toolkit 13.4, cuBLAS uses the 2M algorithm to accelerate ZGEMM and Figure-4 shows the speedup achieved by the 2M + Ozaki-II scheme over the traditional 4M method on GB200.
ZGEMM application performance – BerkeleyGW
The BerkeleyGW code is used by researchers to study physical properties of materials that emerge as a result of how electrons change energy states. It is a massively parallel code that has been used at full scale on leadership class supercomputers. Using GPUs with BerkeleyGW can lead to an 86x performance speedup over the CPU-only implementation and can be even further accelerated with the latest cuBLAS enhancements to ZGEMM.
By default, cuBLAS uses its ADP framework to determine the parameters that will provide accuracy guarantees comparable to, or better than, standard complex floating-point (CFP64) GEMM algorithm in the CHISUM routine of the BerkeleyGW Epsilon module. This is done automatically for users without any application code changes. We evaluated accuracy and performance in the following configurations: Native ZGEMM, Ozaki-II 2M with ADP and Ozaki-II 2M with 53 bits, which all produced the same Epsilon matrix elements for up to 10 decimal places as specified in the NERSC-10 harness.
On GB200, moving from native ZGEMM to Ozaki-II 2M with 53 bits delivers a 3.9x end-to-end speedup. Similarly, we observe a 6.6x speedup on VR NVL72. Overall, we observe a 5.8x speedup going from native ZGEMM on Blackwell to ZGEMM using Ozaki-II 2M with 53 bits on VR NVL72. Figure-5 shows the end-to-end performance speedup running the Epsilon module of BerkeleyGW with a reference problem size (Si998) on 32 Blackwell and Rubin GPUs (2 Grace CPUs + 4 B200 GPUs per node and 2 Vera CPUs + 4 R100 GPUs per node).
DGEMM performance on GB200, RTX PRO 6000
Figure 6, below, shows the speedup of DGEMM with Ozaki-II over native FP64 on GB200. At the smallest shapes, ADP has determined that the Ozaki-II path offers no advantage and dispatched to native FP64, but the performance gains grows steadily with larger M, N and K values reaching 3.9x at the largest shapes measured.
The upper-bound of the performance gain depends on how much integer Tensor Core throughput a GPU has relative to its native FP64 units, and that ratio varies widely across the product line with NVIDIA RTX PRO 6000 Blackwell Server Edition at the far end of that range. Figure 7, below, shows the speedup of DGEMM with Ozaki-II over native FP64 on RTX PRO 6000 Blackwell SE. As can be observed from the figure, Ozaki-II benefits the full spectrum of matrix sizes tested, with even the smallest gaining 10x over native FP64 and the largest reaching 23x.
Getting started #
- Download and install the latest NVIDIA CUDA Toolkit
- Read the floating point emulation section of the cuBLAS documentation for full API details, including cublasSetEmulationStrategy() and manual mantissa bit control.
- Explore the cuBLAS Emulation samples on GitHub
- Evaluate the gains on your own workloads and tell us what you find, we want feedback on both the API design and the precision/performance trade-off. You can reach us on the NVIDIA Accelerated Computing Developer Forum .
Acknowledgments
Jiqun Tu contributed to this post.