Accelerating dense physics solves via FP16 Tensor Core accumulation

Dense linear algebra operations represent the core computational bottleneck across numerical physics workloads, including explicit linear solves, Gaussian Process (GP) regression, reduced-order modeling, and Partial Differential Equation (PDE) Green’s function evaluations. On consumer NVIDIA Ampere GPUs (such as the RTX 3080), these dense matrix multiplications default to half-rate execution on Tensor CoresDedicated matrix-multiply units on NVIDIA GPUs that multiply small fp16 tiles and add up the products, far faster than the general-purpose CUDA cores for dense matmul. when configured with FP32 accumulators.

By enabling FP16 accumulator execution via PyTorch’s allow_fp16_accumulation flag, matrix multiplication throughput increases by 1.91× for matrix inverse applications and 1.84× for PDE Green’s function solves, while maintaining relative errors near $\sim 10^{-3}$ against double-precision (FP64) reference baselines.

Mathematical Equivalence Across Dense Physical Systems

Despite structural differences in underlying physics equations, four distinct scientific computing tasks reduce to identical linear algebra operations:

1. Dense Linear Solves: Operator factorization is performed once, followed by multiple triangular solves across right-hand side (RHS) vector blocks.
2. Gaussian Process / Kernel Ridge Regression: Dense covariance matrices are inverted once, followed by matrix-vector evaluations over evaluation queries.
3. Reduced-Order Models: Physical systems are projected onto reduced basis subspaces, followed by repeated matrix-vector applications across parameter sweeps.
4. PDE Green’s Functions: Discrete differential operators are inverted once to construct Green’s function representations, which are subsequently applied across spatial source distributions.

In each workload, operator construction occurs once, while operator application (a dense GEMMGeneral matrix-matrix multiply, the dense C = A·B kernel that sits at the center of all four workloads.) repeats across time steps, parameter samples, or batch inputs.

Computational Amortization: Formation vs. Application

Profiling an $8192 \times 8192$ matrix inversion demonstrates the performance disparity between operator formation and operator application:

• Operator Formation (torch.linalg.inv): 202.0 ms (executed once in FP32).
• Operator Application ($A^{-1} B$ for $R=2048$): 4.46 ms in FP32-accumulator mode vs. 2.34 ms in FP16-accumulator mode.

Because operator formation occurs once at initialization, optimizing application latency dominates total wall-clock time as the number of RHS evaluations grows.

When $R \ge 100$ applications are performed per formation step, execution time becomes entirely application-dominated, and the speedup ratio approaches $4.46 / 2.34 = 1.91\times$.

NVIDIA Ampere Tensor Core Accumulator Architecture

NVIDIA Tensor Cores process matrix multiply-accumulate (MMAMatrix multiply-accumulate: the tensor-core primitive that multiplies two input tiles and adds the product into a running accumulator. A GEMM is built from many of these.) primitives by multiplying $16 \times 16$ FP16 input tiles and accumulating products into a running sum:

On consumer NVIDIA Ampere hardware (microarchitecture sm_86, e.g., RTX 3080), Tensor Cores support two accumulation modes:

• FP32 Accumulation (fp32-accumulateSumming the running dot product in single precision. Numerically safe on any input, but on sm_86 GeForce it issues at half the rate of fp16 accumulation.): Accumulates intermediate products in FP32 precision. On consumer Ampere, this instruction path executes at half hardware throughput.
• FP16 Accumulation (fp16-accumulateSumming the running dot product in half precision. Same fp16 inputs, same cores, but on sm_86 it issues at the full MMA rate.): Accumulates intermediate products in FP16 precision, issuing at full hardware clock rates.

By default, cuBLASNVIDIA's dense linear-algebra library, the GEMM backend PyTorch calls. It defaults to fp32 accumulation because that is correct on any input. selects FP32 accumulation to prevent potential numerical overflow across arbitrary workloads.

Precision Trade-Offs and Error Floor Bounds

Switching from FP32 to FP16 accumulation introduces small numerical truncation errors during intermediate dot-product summation. However, physical surrogate models, kernel regressions, and discrete PDE solvers already contain inherent discretization noise and input parameter uncertainty.

For well-conditioned physical operators ($\kappa(A) \le 200$), the error floor introduced by FP16 accumulation remains bounded near $\sim 10^{-3}$ relative error. If input data uncertainties exceed $10^{-3}$, using higher-precision accumulators consumes compute throughput without improving model accuracy.

Benchmark Results and Experimental Design

FP16 accumulation is enabled via PyTorch runtime options:

import torch

# Enable full-rate FP16 accumulation on Tensor Cores (PyTorch >= 2.7)
torch.backends.cuda.matmul.allow_fp16_accumulation = True

# Execute operator application GEMM
X = A_inv @ B

To isolate accumulator throughput gains from architectural changes, benchmarks were evaluated across three configurations on an NVIDIA RTX 3080 GPU (10 GB):

1. CUDA Core FP32: Standard FP32 execution on general CUDA cores (TF32 disabled).
2. Tensor Core FP32-Accumulate: FP16 input tensors with FP32 accumulator (cuBLAS baseline).
3. Tensor Core FP16-Accumulate: FP16 input tensors with FP16 accumulator (optimized configuration).

torch.manual_seed(0)

# Configuration 1: Pure FP32 (CUDA Cores)
torch.backends.cuda.matmul.allow_tf32 = False
X_fp32 = A_fp32 @ B_fp32

# Configuration 2: FP16 Inputs, FP32 Accumulator (Tensor Cores)
torch.backends.cuda.matmul.allow_fp16_accumulation = False
X_acc32 = A_fp16 @ B_fp16

# Configuration 3: FP16 Inputs, FP16 Accumulator (Tensor Cores)
torch.backends.cuda.matmul.allow_fp16_accumulation = True
X_acc16 = A_fp16 @ B_fp16

Changing the accumulator precision yields a 1.91× speedup on matrix inversion applications and a 1.84× speedup on PDE Green’s function applications compared to matched FP16-input baselines.

Relative Latency Comparison

Disambiguating Cross-Path vs. Same-Kernel Speedups

Comparing FP16 accumulation directly against CUDA Core FP32 execution yields a nominal speedup of $5.35\times$ ($12.51\text{ ms} \to 2.34\text{ ms}$). However, this metric combines two distinct architectural transitions:

1. Moving computation from general CUDA cores to Tensor Cores ($22.0 \to 61.6\text{ TFLOP/s}$, a $2.8\times$ gain).
2. Frictional accumulation rate doubling within Tensor Cores ($61.6 \to 117.5\text{ TFLOP/s}$, a $1.91\times$ gain).

To isolate the specific impact of the accumulator toggle, reported speedups ($1.91\times$ and $1.84\times$) reflect comparisons against matched Tensor Core FP16-input baselines.

Numerical Accuracy vs FP64 Gold References

To measure numerical precision loss, outputs from each execution mode were compared against double-precision (FP64) reference solutions ($X_{\text{gold}}$) computed via torch.linalg.solve:

Relative error $\|X - X_{\text{gold}}\| / \|X_{\text{gold}}\|$ and residual norms $\|A X - B\| / \|B\|$ confirm that FP16 accumulation maintains numerical errors bounded between $2.3 \times 10^{-3}$ and $3.4 \times 10^{-3}$ for well-conditioned operators.

Validity Boundaries and Condition Number Constraints

The application of FP16 accumulation is bounded by three specific numerical conditions:

1. Condition Number Sensitivity: Test operators must remain well-conditioned. The PDE benchmark evaluates a screened PoissonThe operator I − α∇², a Poisson operator shifted away from singularity by a small α so its inverse stays well-conditioned. operator ($I - \alpha \nabla^2$ with $\alpha = 2 \times 10^{-3}$), yielding a condition numberHow much a matrix amplifies error when inverted. Near 1 is benign; a large value magnifies whatever the accumulator rounds off. $\kappa(A) \approx 146$. Unshifted Poisson operators ($\nabla^2 u = f$) exhibit ill-conditioned spectra ($\kappa(A) \gg 10^4$), causing FP16 truncation errors to amplify significantly.
2. Iterative Time-Stepping Accumulation: Explicit time-integration schemes (e.g., forward Euler PDE integration) compound rounding errors across sequential steps. FP16 accumulation should be restricted to single-pass matrix-vector applications or stationary solves.
3. Residual Refinement Limits: Standard iterative refinement ($X_{k+1} = X_k + A^{-1}(B - A X_k)$) fails to recover lost precision under FP16 accumulation because small residual updates round to zero within FP16 dynamic ranges.

Verification and Reproducibility Protocol

All benchmark measurements adhere to the following experimental constraints:

• Reference Solves: Evaluated against FP64 reference ground truth (torch.linalg.solve).
• Timing Measurement: Captured using asynchronous GPU CUDA events (torch.cuda.Event) across 10 trials (median reported).
• Environment: NVIDIA RTX 3080 GPU (Ampere sm_86, 10 GB), PyTorch 2.12, CUDA 13.0.

Generalization to Dense Matrix Operations Across Domains

The performance characteristics of FP16 accumulator execution extend to any compute-bound dense GEMM workload tolerant of $\sim 10^{-3}$ relative error:

Conclusion

Enabling FP16 accumulation via PyTorch’s allow_fp16_accumulation option recovers 1.84–1.91× throughput on consumer NVIDIA Ampere GPUs for compute-bound dense matrix operations. For well-conditioned operators ($\kappa(A) \le 200$), the resulting numerical error floor remains bounded at $\sim 10^{-3}$, providing significant acceleration for surrogate modeling, uncertainty quantification, and numerical PDE solves.

Glossary