cuBLAS FP64 Algorithms: Trading Extra Products for Faster Scientific GEMM

NVIDIA’s November 13 cuBLAS post describes new FP64 algorithm paths in CUDA Toolkit 13.4.2. ^[1] The key idea is algorithmic: represent high-precision values as combinations of lower-precision pieces, perform more products efficiently, then reconstruct a result within an error strategy. This can make some dense matrix operations faster, but it is not a blanket replacement for conventional double-precision arithmetic.

Ozaki-II and 2M

Ozaki Scheme II decomposes inputs into integer pieces and uses modular or residue-style arithmetic to recover products. Its usefulness depends on matrix dimensions, data ranges, and the required error bound. The vendor’s Ozaki-II 2M approach targets complex GEMM and uses two real GEMMs after decomposing inputs into integer pieces; those intermediate integer operations are exact. This differs from standard 4M and classical 3M complex multiplication. The scalar three-multiplication lower-bound argument does not describe this exact integer-matrix algorithm. cuBLAS selects algorithmic paths dynamically under its precision controls.

NVIDIA reports speedups up to 5.5x for DGEMM on a Vera Rubin NVL72 configuration and up to 3.9x on GB200, with up to 5x for complex ZGEMM in selected cases. These are vendor measurements for chosen sizes and hardware. They do not predict an application’s speedup without its matrix shapes, data, compiler/runtime versions, and error tolerance.

Validate the numerical contract

For a scientific application, compare results against a trusted FP64 reference using the application’s own residual, conservation, or convergence criteria. Sweep matrix dimensions and value distributions, including ill-conditioned cases. Record absolute and relative error, repeatability, and any downstream change in solver iterations. A fast GEMM that shifts a simulation outside its validated tolerance is not an optimization.

Also profile the whole application. GEMM may not dominate runtime; memory traffic, communication, sparse kernels, and synchronization can limit end-to-end gains. Keep a conventional algorithm available as a reference and document which precision mode was used for each published result.

The development is notable because it puts more adaptive numerical methods behind a familiar library interface. The engineering responsibility moves to making the accuracy/performance tradeoff explicit and testable.

Error bounds belong beside speed numbers

Algorithmic precision is not adequately described by the label FP64 alone. Decomposition and reconstruction introduce a numerical pathway different from a conventional direct double-precision product. The relevant question is whether its error stays within the application’s validated tolerance for the matrices it actually uses. The Ozaki Scheme II research literature motivates the approach, while cuBLAS packages an implementation and selection policy for specific NVIDIA platforms. A library benchmark does not substitute for checking a simulation’s numerical contract.

Build a test matrix suite from production behavior: square and rectangular shapes, sizes around algorithm-selection thresholds, varied magnitudes, cancellation-heavy values, and matrices from real simulation states. Compare output to both a trusted implementation and the application-level invariants. Track relative error where the reference is nonzero and absolute error near zero; one metric can be misleading for mixed-scale values. For iterative solvers, compare convergence history and final residual, not just one matrix output.

The 2M complex method reduces the work to two real GEMMs through its integer decomposition and exact intermediate operations. Whether it helps depends on matrix shape, conversion/reconstruction work, memory movement, and the library path selected. Profile kernels and the whole application. If GEMM is 15 percent of elapsed time, even a large GEMM speedup has a bounded overall effect.

Record the Toolkit version and algorithm/precision controls so that a later library update can be checked for changed selection behavior. NVIDIA’s headline figures are useful leads for where to test first; they should remain attributed to the benchmark configuration that produced them.

For related systems context, see GPU scheduling and contention; resource contention can affect measured throughput, although it does not alter the numerical method’s accuracy contract.

Keep algorithm selection observable

If a library selects an algorithm automatically, performance changes can be caused by a shape threshold or changed runtime version rather than an application code edit. Include representative matrix sizes in performance monitoring and capture the selected mode when the library exposes it. Compare against a known conventional path so a speedup can be attributed to the new algorithm rather than different input generation, affinity, or system load.

For reproducibility, record input-generation seeds, scaling, matrix dimensions, transposition, leading dimensions, workspace, stream behavior, GPU model, and Toolkit build. For mixed scientific workloads, report both microbenchmark performance and solver-level outcomes. Our deterministic training analysis covers a related tension: performance optimizations need an explicit reproducibility and correctness contract.