The hot loop of physics is a dense matrix multiply.

Dense linear solves, kernel and covariance methods, reduced-order models, and PDE Green's-function applies share one operation, multiplying a big matrix by a block of many columns. On a consumer RTX 3080 that multiply runs at half the speed the chip can reach, because the standard math library keeps the running sum in the safer of two precisions. Where about 1e-3 relative error is acceptable, one flag doubles the rate.

Written to be read cold. Each term is explained where it appears, with mechanics behind the expandable panels and a glossary at the end. Every number was measured on one RTX 3080 (Ampere, sm_86).

Matrix-inverse apply 1.91× faster

PDE Green's-function apply 1.84× faster

at ~1e-3 relative error · one flag, no custom CUDA

The shared setup

A GEMM (general matrix-matrix multiply) computes C = A·B, where every entry of C is a dot product of a row of A with a column of B. It is the densest, most regular numerical work a computer does, and GPUs are built to run it fast.

The hardware. One consumer NVIDIA RTX 3080 (Ampere, compute capability sm_86). Its tensor cores multiply fp16 inputs while keeping the running sum in either fp16 or fp32. On this chip the fp16-accumulate path issues at full rate and the fp32-accumulate path at half rate. cuBLAS, the standard library under PyTorch, picks fp32 accumulation by default, so the ordinary path runs at half the speed the silicon can reach.

01 where the time goes

Why physics compute is a matrix multiply

Four workloads that look different on paper run the same operation underneath. Each forms an operator or factorization once, then applies it to a block of many right-hand sides, and that block apply is a large dense GEMM.

Dense linear solve

Solving A x = b for a single right-hand side is cheap. Real problems solve A X = B for hundreds or thousands of columns at once, one per source, target, or parameter sample. Factor or invert A once, then every solve is X = A⁻¹B, a dense matrix multiply across the whole block.

Covariance and kernel matrices

Gaussian-process regression, kernel ridge regression, RBF field interpolation, and data assimilation build a dense n×n kernel matrix K, where K_ij measures similarity between points i and j. Fitting and prediction solve against K + λI. The matrix is genuinely dense, and prediction multiplies the solved operator against many query blocks. This is the exact system measured below.

Reduced-order models

A reduced-order model compresses a large simulation onto a small basis of a few hundred modes, then advances the small dense system in time. Each step multiplies the reduced operator by the current state across many parameter settings or ensemble members, again a dense matmul. This one is explained here rather than benchmarked.

PDE Green's-function apply

For a linear PDE with fixed coefficients, the solution operator is a Green's function that returns the field from the sources. Discretized, that operator is a dense matrix M, and solving for many source configurations is X = M F, one column per configuration. The screened-Poisson operator below is measured this way.

Why "many right-hand sides" puts the work on the tensor cores

A matmul's speed limit depends on its arithmetic intensity, the ratio of math done to bytes moved from memory. Applying an operator to a single column reads the whole matrix to do very little math, so it waits on memory. Applying it to a block of R columns reuses each loaded matrix value R times, so the math dwarfs the traffic.

For the inverse apply below (n=8192, R=2048), one call does 2·n²·R ≈ 2.7×10¹¹ floating-point operations against roughly 0.2 GB of memory traffic, about 1,400 operations per byte. That sits far inside the regime where math throughput, not memory bandwidth, is the ceiling. So the operation runs at the tensor cores' arithmetic limit, and the accumulator rate is what sets that limit.

02 the lever

One flag flips the accumulator

A matmul sums thousands of products into a running total called the accumulator. Its precision is a separate choice from the input precision, and on this GPU it sets the speed.

Tensor cores take fp16 inputs and keep the running sum in fp16 or fp32. fp32 holds about seven decimal digits and rounds gently; fp16 holds about three and rounds hard. cuBLAS keeps the sum in fp32 by default so results stay accurate for any input. On sm_86 GeForce silicon the fp32-accumulate path issues at half the rate of the fp16-accumulate path, so the default runs the tensor cores at half speed.

Choosing fp16 accumulation runs the same cores at full rate. In PyTorch it is one line, no custom CUDA and no kernel to build, on version 2.7 or newer:

# fp16 inputs already; flip the accumulator to full rate
torch.backends.cuda.matmul.allow_fp16_accumulation = True

The cost is numerical. The fp16 sum accumulates rounding and lands near 1e-3 relative error on these problems, against about 3.5e-4 for fp32 accumulation and about 1e-6 for a true-fp32 solve. The rest of this page measures the speed gained and the accuracy given up.

Why the fp16 sum lands near 1e-3, and why fp32 accumulation is the safe default

A dot product of length n adds n products one after another. Once the partial sum grows large, a new product smaller than the sum's fp16 step size rounds away entirely, a loss called swamping. fp16 carries about 10 mantissa bits, so after thousands of adds the discarded low bits total a relative error near 1e-3. The measured applies land at 2.3e-3 (PDE, inner dimension 9,216) and 3.4e-3 (solve, inner dimension 8,192); the exact figure depends on the operator's conditioning and the spread of magnitudes in the sum, not on the length alone.

fp32 carries 23 mantissa bits, so the same sums stay near 3.5e-4 to 3.9e-4, and the true-fp32 CUDA-core path reaches about 1e-6. Keeping the accumulator in fp32 is the safe default because it protects any input; flipping it to fp16 is a deliberate trade for full-rate throughput where the extra error is below what the problem already carries.

One lever, three domains

The same fp16-accumulate flag carries across a family of projects. The flag is identical; the workloads differ.

this repo

physics-frontier

Dense physics solves on an RTX 3080. Matrix-inverse apply 1.91×, PDE Green's-function apply 1.84×, each at about 1e-3 relative error.

sibling · same lever

nerf-frontier

The same full-rate fp16-accumulate flag, applied to neural-field rendering. Method and measurements live in its own repo: nerf-frontier.

sibling · same lever

search-frontier

The same flag, applied to vector search. Method and measurements live in its own repo: search-frontier.

03 measured

Two workloads, measured on one RTX 3080

Both apply a dense operator to a block of many right-hand sides. The fair baseline is the fp32-accumulate tensor path, the cuBLAS default; the lever is the same kernel with the fp16 accumulator. True fp32 is the CUDA-core path, shown for context, and comparing the lever against it overstates the win.

The bars below replay the two applies of the inverse workload at their measured wall-clock times. Press replay to watch fp16-accumulate finish first.

apply A⁻¹ to 2048 right-hand sides · fp32-acc vs fp16-acc

—0.00 ms · — TFLOP/s

Matrix inverse · apply A⁻¹ to 2048 right-hand sides

A dense Gaussian-process / RBF kernel matrix A = K + λI, 8,192 centers, condition number 101. The inverse is formed once in fp32 (202 ms), then each apply to a 2,048-column block is a 274.9 GFLOP GEMM. Bars show achieved throughput (longer is faster); the right column is wall-clock per apply.

PDE · screened-Poisson, 4096 sources

The screened-Poisson operator (I − α∇²) on a 96×96 grid, 9,216 unknowns, condition number 146. Form the dense Green's operator M = (I − α∇²)⁻¹ once, then each apply to 4,096 source fields is a 695.8 GFLOP GEMM.

≈1.9×
at one flag · same tensor cores

The fp16 accumulator runs the cores at full rate. 1.91× on the inverse apply (117.5 vs 61.6 TFLOP/s) and 1.84× on the Green's-function apply (120.0 vs 65.2 TFLOP/s), each at about 1e-3 relative error against an fp64 gold solve. The gain is throughput the fp32-accumulate default already left unused on this chip.

Why the apply, not the factorization, is the number that matters

Forming the inverse of the 8,192×8,192 kernel costs 202 ms once, in fp32. Each apply to a 2,048-column block costs 2.34 ms. The factorization is paid a single time; the apply runs again for every new block of right-hand sides, every batch of sources, every ensemble member. Over a sweep of many blocks the apply is where the wall-clock accrues, which is why it is the operation put on the full-rate path.

The true-fp32 CUDA-core path finishes the same apply in 12.51 ms (solve) and 30.09 ms (PDE). Against it the fp16-accumulate path looks 5.35× and 5.19× faster, but that compares tensor cores to CUDA cores. The tensor path itself, fp32-accumulate to fp16-accumulate, is 1.91× and 1.84×.

04 accuracy

When ~1e-3 is fine, and when it is not

The lever trades numerical precision for speed. Whether the trade is acceptable depends on the problem, not the hardware.

When ~1e-3 is acceptable

fp16 accumulation lands at 2.3e-3 (PDE) and 3.4e-3 (solve) relative to an fp64 gold solve; fp32 accumulation holds 3.5e-4 to 3.9e-4. For Gaussian-process surrogates, uncertainty-quantification ensembles, RBF field interpolation, and machine-learning-for-physics, the inputs already carry more than 1e-3 of noise, so the accumulator error sits below the problem's own floor. Both operators here are well-conditioned (condition number 101 and 146) and each apply is independent, so the error does not amplify and nothing compounds across solves.

If you need better than 1e-3

The fp32-accumulate tensor path holds 3.5e-4 to 3.9e-4 and still runs about 2.8× the CUDA-core path (4.46 vs 12.51 ms on the solve, 10.67 vs 30.09 ms on the PDE). It is the accurate middle setting when 1e-3 is too loose but the true-fp32 path is slower than you want.

Where the lever does not apply

Ill-conditioned operators. A well-conditioned operator keeps the fp16 error bounded; a near-singular one amplifies it. The plain Poisson inverse magnifies small perturbations, so the benchmark uses the screened operator (I − α∇²), which is symmetric positive-definite and well-posed. On a stiff or near-singular system, stay on the fp32 path.

Long time-marching. Stepping a PDE forward in fp16 drifts. The smoothest mode of the screened operator barely changes per step, so its rounding error is never damped and accumulates over a long trajectory. The lever is for independent solves, where each apply carries its own error.

Sub-1e-3 dense inverses. The full fp16 path stores the inverse in fp16 and accumulates in fp16, landing at 3.4e-3 relative error and a 7.9e-3 residual. Improving that would need iterative refinement, and for a dense inverse each refinement step is itself a full GEMM, so a correction costs as much as the original solve. Newton–Schulz inversion fares worse in fp16, its early iterations increment below fp16 resolution and stall.

Why the smoothest mode barely damps, and why Newton–Schulz stalls

The screened operator (I − α∇²) has eigenvalues near 1 for the smoothest modes and larger for the roughest. Its inverse multiplies the smoothest mode by nearly 1, so an fp16 time step neither grows nor shrinks it, and that mode's rounding error persists step after step until it dominates a long trajectory. An independent solve never sees this, because it applies the operator once and reports a single result.

Newton–Schulz builds an inverse by repeated matmuls, X ← X(2I − AX). Its early iterations make tiny corrections, and in fp16 those corrections fall below the accumulator's resolution and round to zero, so the iteration never gets moving. Both cases are reasons to keep those specific paths in fp32; the measured lever covers the block solve, where neither failure occurs.

The apply is the recurring cost

Form A⁻¹ once in 202 ms; each 2,048-column block is a 2.34 ms GEMM. The block apply, run for every batch, is where wall-clock accrues.

Full rate is one flag away

fp16 accumulation runs sm_86 tensor cores at 2× the fp32-accumulate default. Measured 1.91× and 1.84×, no custom CUDA.

~1e-3 is the price

fp16-acc lands at 2.3–3.4e-3 vs fp64, below the noise floor of surrogates and UQ, above what stiff solves and time-marching tolerate.

Measure against the right baseline

Versus the CUDA-core path the apply looks 5.2–5.4× faster, but that compares tensor cores to CUDA cores. The lever itself is 1.84–1.91×.

05 glossary

Every term and tool, defined

GEMM: General matrix-matrix multiply, C = A·B. The dense, regular core of most numerical linear algebra and the operation these workloads reduce to.
Dense matrix: A matrix whose entries are mostly non-zero, so every one is stored and multiplied. The opposite of a sparse matrix, which skips zeros.
Tensor core: A GPU unit that does a small matrix-multiply-and-accumulate per instruction, fed fp16 inputs, summing into an fp16 or fp32 accumulator.
Accumulator: The running total a matmul adds its products into. Its precision is chosen separately from the input precision and, on sm_86, sets the throughput.
Accumulate precision: Whether the running sum is kept in fp16 or fp32. fp16-accumulate is full rate on sm_86; fp32-accumulate is half rate but more accurate.
fp16 / fp32 / fp64: 16-, 32-, and 64-bit floating-point. Roughly 3, 7, and 16 decimal digits of precision. fp64 is used here only as the gold reference.
sm_86 / consumer Ampere: The compute capability of GeForce Ampere GPUs, including the RTX 3080, where fp16-accumulate tensor math runs at twice the fp32-accumulate rate.
RTX 3080: A consumer NVIDIA GPU (Ampere, sm_86, 10 GB). The single card every number here was measured on.
cuBLAS: NVIDIA's dense linear-algebra library, the GEMM backend under PyTorch. It defaults to fp32 accumulation for safety.
allow_fp16_accumulation: The PyTorch flag (torch.backends.cuda.matmul.allow_fp16_accumulation) that switches cuBLAS to the full-rate fp16 accumulator. Requires PyTorch 2.7+.
TF32: A 19-bit tensor-core format Ampere can substitute for fp32. Disabled in these benchmarks so the context baseline is a true fp32 CUDA-core solve.
Right-hand side (RHS): A column of the system A X = B. Solving for many at once turns the apply into a wide, compute-bound GEMM.
Dense linear solve: Finding X in A X = B for a dense A. Factor or invert once, then apply to every block of right-hand sides.
Condition number: How much a matrix amplifies input error. Low means well-conditioned and safe for fp16; the operators here are 101 and 146.
Covariance / kernel matrix: A dense n×n matrix of pairwise similarities between data points. The object solved against in GP and kernel methods.
GP / RBF kernel: Gaussian process with a radial-basis-function kernel, K_ij = exp(−‖x_i−x_j‖²/2ℓ²). The dense system in the inverse benchmark.
Kernel ridge regression: Regression that solves against K + λI, one of the methods whose predict step is the block apply measured here.
Reduced-order model (ROM): A simulation compressed onto a small basis, then advanced as a small dense system. Each step is a matmul over many parameter settings.
Green's function: The solution operator of a linear PDE: apply it to the sources and get the field. Discretized, it is a dense matrix.
Screened-Poisson: The operator (I − α∇²), also called modified-Helmholtz. Symmetric positive-definite and well-conditioned, so its inverse-apply is numerically clean.
Arithmetic intensity: Math done per byte moved from memory. High intensity (many RHS) makes an apply compute-bound, so tensor-core rate is the ceiling.
Compute-bound / memory-bound: Limited by math throughput versus by memory bandwidth. The block applies here are compute-bound, so the accumulator rate sets the speed.
Relative error: ‖X − X_gold‖ / ‖X_gold‖ against an fp64 solve. The accuracy figure traded for speed: ~1e-3 for fp16-acc, ~1e-6 for true fp32.
Residual: ‖A X − B‖ / ‖B‖, how well X actually satisfies the system. The fp16 path's residual is 7.9e-3 on the inverse apply.
Iterative refinement: Correcting a solve with extra steps. For a dense explicit inverse each step is a full GEMM, so it cannot cheaply pass the fp16 floor.
Newton–Schulz iteration: Building an inverse by repeated matmuls, X ← X(2I − AX). Its small early corrections round to zero in fp16 and stall.
TFLOP/s: Trillions of floating-point operations per second. Achieved throughput, the natural speed axis for a compute-bound GEMM.
nerf-frontier / search-frontier: Sibling projects that apply the same fp16-accumulate flag to neural-field rendering and to vector search, respectively.