Dense linear algebra is the hot loop of computational physics, and the same shape hides inside problems that look unrelated. A dense linear solve, a Gaussian-process regression, a reduced-order model, a PDE Green’s function: strip each to the operation it repeats and one thing is left, a large dense matrix multiply. On a consumer RTX 3080 that multiply runs at half the rate the tensor cores can sustain, for a numerical safety margin these workloads never asked for. One line of PyTorch takes the other half back.
The gain is 1.91× on a dense inverse and 1.84× on a PDE, at a relative error near 1e-3 against a double-precision reference. The same measurement can be made to read 5.3×, and refusing to quote that larger number is the spine of the whole piece. The number that survives is the one that holds the kernel fixed and changes only the thing under test.
Four workloads, one shape
The starting move is not to optimize any one solver. It is to notice that four workloads that look different collapse to the same loop. Each forms an operator once, then applies it to a block of many right-hand sides. The formation is a one-time cost. The apply is a large dense GEMM, and it is the part that repeats.
That reframing is what makes a single lever worth chasing. If the recurring cost of a dense solve, a kernel regression, and a Green’s-function evaluation are literally the same matmul, then one change to that matmul pays four times.
Why do four different workloads reduce to the same matmul?
Strip each candidate to its recurring hot loop and the loop is identical: form the operator once, apply it to many columns.
- Dense linear solve. Factor the operator once, then apply the factorization to many right-hand sides. The factorization is amortized; the triangular applies are the repeated cost.
- Gaussian-process / kernel-ridge regression. Form and invert a dense covariance matrix once, then apply that inverse to many query vectors. Every prediction is an apply.
- Reduced-order models. Project onto a reduced basis, then apply the reduced operator repeatedly across parameter samples. (Described here by structure, not benchmarked in this work.)
- PDE Green’s functions. Form the dense operator inverse once, then apply it to many source fields. Each source is a column of the right-hand side.
The unifying statement is short: form the operator once, apply it to many RHS. The formation happens once; the apply happens per batch, per ensemble member, per parameter sample. So the apply is the thing to optimize, and the apply is a dense matrix multiply.
The apply is the recurring cost
The split is not rhetorical. It is measured. Forming the 8,192² inverse costs 202 ms once, computed in
fp32 through cuSOLVER via torch.linalg.inv. Each 2,048-column apply after that costs 2.34 ms. The
formation is expensive and happens a single time; the apply is cheap and happens without end. Optimizing
the 202 ms buys nothing after the first call. Optimizing the 2.34 ms compounds.
Where does the 202 ms go, and why optimize the 2.34 ms instead?
The 202 ms is a one-time factorization or explicit inversion of the operator. In a real run it is paid once at setup, then never again while the same operator is reused across a batch of right-hand sides, an uncertainty-quantification ensemble, or a sweep of parameter samples.
The 2.34 ms apply is paid on every one of those. For any run with more than roughly a hundred applies per formation, total time is dominated by the applies, and the formation rounds to a constant. This is the amortization argument that decides where to look: a 2× on the 2.34 ms is a near-2× on the whole run; a 2× on the 202 ms vanishes into the first iteration.
One wrinkle: the 202 ms and the 2.34 ms are timed on different clocks. The formation uses a
single wall-clock call around one torch.linalg.inv plus a device synchronize; the applies use paired
CUDA events and a median of ten. Both are internally consistent, and they are only ever compared as
orders of magnitude, not subtracted.
Why the apply runs at half rate
With the target fixed on the apply, the next question is why that matmul is slow, and the answer is not an algorithm to invent. It is a default to override.
Tensor cores take fp16 inputs and sum the running dot product into an accumulator that can be fp16 or fp32. On consumer Ampere (sm_86) the fp32-accumulate path issues at half the rate of the fp16-accumulate path. cuBLAS, the GEMM backend PyTorch calls, picks fp32 accumulation by default, because it is safe on any input regardless of whether the workload needs that much precision. Half the tensor-core throughput sits idle on every dense solve, paying for accuracy that the surrogate-modeling workloads this is for never read.
That leaves one loose end. Raising the accumulator’s issue rate only helps if the apply is compute-bound. If it were memory-bound, the faster accumulator would stall waiting on HBM and the wall-clock would not move. So the compute-bound claim has to be checked, not assumed.
What is an accumulator, and why does fp32 cost double on sm_86?
A tensor-core matmul multiplies fp16 (or bf16) input tiles and adds each product into a running sum. That running sum is the accumulator, and its precision is a separate choice from the input precision. Inputs stay fp16; the accumulator can be fp16 or fp32.
On sm_86 GeForce silicon, the fp32-accumulate mode issues at half the throughput of the fp16-accumulate mode. cuBLAS defaults to fp32 accumulation because it is the numerically conservative choice that is correct on any input. The lever holds the fp16 inputs fixed and changes only the accumulator, so nothing about the data or the kernel shape moves; only the internal summation precision does. That is what makes it one flag rather than a rewrite.
The scope of this half-rate default is consumer Ampere specifically. Whether datacenter Ampere or newer architectures carry the same default is not measured in this work, so nothing here should be read as an A100 or Hopper number.
Is the apply really compute-bound? The arithmetic-intensity budget.
For the inverse apply with n=8192 and R=2048 columns, the work is 2·n²·R ≈ 2.7×10¹¹ floating-point operations. The memory traffic is the fp16 operands and result, roughly 0.2 GB. That puts the arithmetic intensity around 1,400 operations per byte (an approximation, from the roughly-2.7×10¹¹ FLOP over roughly-0.2 GB estimate above), far above the roofline ridge where a workload flips from memory-bound to compute-bound.
Compute-bound is exactly the regime where the accumulator’s issue rate is the binding constraint. The apply is not waiting on memory, so doubling the rate at which the tensor cores can retire fused multiply-adds converts directly into wall-clock, instead of hiding behind HBM traffic. This is the check that says the lever will actually land, and it is why the same flag does almost nothing for the sparse, memory-bound workloads discussed near the end.
Why the accumulator, and why it is safe to turn
A matmul running at half its achievable rate usually invites a menu of fixes: retile for better occupancy, batch more right-hand sides per launch, move to a cheaper input format, or hand-write a fused kernel. The diagnosis crosses most of that menu out before it starts. The apply is compute-bound, so there is no memory-traffic problem to rearrange around. cuBLAS already autotunes the tile for a GEMM this large, so there is little occupancy left to win. The inputs are already fp16, so there is no input precision left to shed. The one quantity the half-rate fact actually points at is the accumulator issue rate, and the accumulator is the single knob sitting on it. That is why the fix is one line and not a kernel project. The bottleneck and the knob turn out to be the same object.
Turning that knob only helps if the precision it spends is precision the workload was not using. fp16 accumulation sums the running dot product in half precision, and the reflex worry is that a physics solve cannot afford it. Reasoning from what the error actually is settles the worry. The accumulator adds a bounded error floor, set by how many terms get summed and how well-conditioned the operator is, and for a well-conditioned dense apply that floor lands near 1e-3. The targets here, surrogate models, uncertainty-quantification ensembles, and kernel regressions, already carry input noise from discretization, model form, and measurement that sits well above 1e-3. When the floor the accumulator adds is smaller than the noise already in the inputs, holding it any lower buys nothing the problem can read. That is a claim about a specific class of well-conditioned operators, not about fp16 in general, and the fp64 gold solve later in the piece is what turns it from an argument into a measured bound.
One decision is left, and it is about measurement. The speedup is being credited to a single variable, the accumulator, so the comparison has to move that variable alone. The tempting baseline is the naive fp32 path a default run would use, but it differs from the lever in two ways at once, the compute path and the accumulator, so any speedup read against it credits the accumulator for a change the accumulator never made. The comparison that isolates the knob runs the same tensor-core kernel on the same fp16 inputs and changes only the accumulator. Choosing that harder baseline over the flattering one is what makes the resulting number attributable to the thing under test, and it is why the figure that follows lands well below what a looser setup would print.
The lever, measured three ways
The intervention is a runtime toggle, flipped right before the apply:
torch.backends.cuda.matmul.allow_fp16_accumulation = True # PyTorch >= 2.7
X = Ainv @ B # apply the inverse to a block of RHS, at full rateNo custom CUDA, no rebuild, no change to the operator or the data. It requires PyTorch 2.7 or newer, and every number here is from torch 2.12 with CUDA 13 on a single RTX 3080 (10 GB).
To attribute the gain to the accumulator and nothing else, each workload runs three ways, holding the fp16 inputs fixed and moving only the accumulator:
- True fp32 on the CUDA cores, with TF32 explicitly disabled. The accurate reference.
- fp32-accumulate on the tensor cores, the cuBLAS default. The fair baseline.
- fp16-accumulate on the tensor cores. The lever.
How the three configurations are pinned to one variable
The whole comparison depends on changing exactly one thing between rungs. The harness pins the other variables explicitly:
torch.manual_seed(0) # same operator, same RHS every run
# rung 1: true fp32 on the CUDA cores, TF32 off so it is genuinely fp32
torch.backends.cuda.matmul.allow_tf32 = False
X_fp32 = A_fp32 @ B_fp32
# rung 2: fp32-accumulate on the tensor cores (cuBLAS default), fp16 inputs
torch.backends.cuda.matmul.allow_fp16_accumulation = False
X_acc32 = A_fp16 @ B_fp16
# rung 3: fp16-accumulate on the tensor cores, same fp16 inputs
torch.backends.cuda.matmul.allow_fp16_accumulation = True
X_acc16 = A_fp16 @ B_fp16Inputs are fp16 for both tensor-core rungs, so the only difference between rung 2 and rung 3 is the accumulator. Rung 1 exists to give an accurate reference at a known cost, and TF32 is switched off so it is a true CUDA-core fp32 path rather than a TF32-accelerated one. Each rung is timed with paired CUDA events, one untimed warmup, then the median of ten repetitions, so Python-side dispatch overhead does not pollute the GPU number.
Same kernel, only the accumulator changed: 1.91× on the inverse apply and 1.84× on the PDE apply, over the fair baseline.
Watch one apply race the other
The ladders show the endpoints. The gap is easier to feel as a race, one lane per accumulator, each filling at its measured wall-clock so the fp16 lane crosses the line first.
The fair number, not the flattering one
A bigger number is available, and it is wrong to quote. Compare fp16-accumulate against the true-fp32 CUDA-core path and the inverse apply reads 5.35×, the PDE 5.19×. That figure moves two things at once: CUDA cores to tensor cores, and fp32 accumulation to fp16 accumulation. Attributing all of it to the accumulator overstates the lever by roughly 3×.
The fair baseline exists to strip out the first change. It runs on the same tensor cores with the same fp16 inputs, so the only difference from the lever is the accumulator. Against it, the number is 1.91× and 1.84×. That is the figure that survives a like-for-like comparison, so that is the headline.
Wouldn't real-world fp32 be faster anyway, with TF32 on?
Yes, and it makes the 5.3× look even worse, which is another reason not to lead with it. The true-fp32 reference here is measured with TF32 deliberately switched off, so it is a genuine CUDA-core fp32 path. Many real deployments leave TF32 on by default, which would make an out-of-the-box “naive fp32” faster than the 12.51 ms and 30.09 ms measured here.
A faster fp32 reference means the cross-path ratio in practice would be smaller than 5.3×, not larger. So the already-discounted 5.3× is if anything an overstatement of the overstatement. The fair baseline sidesteps the whole issue by never involving the fp32 path in the headline at all.
The accuracy trade, in three tiers
Speed is only half the trade. The other half is measured against an fp64 gold solve, so the cost is bounded rather than hidden. Each rung buys a known amount of accuracy for a known amount of throughput:
| Path | Throughput | Rel. error vs fp64 | When to reach for it |
|---|---|---|---|
| True fp32 (CUDA core) | 22.0 TFLOP/s | 1.3e-6 | Stiff or near-singular operators |
| fp32-accumulate (tensor) | 61.6 TFLOP/s | 3.9e-4 | Accurate middle, ~2.8× for near-free |
| fp16-accumulate (tensor) | 117.5 TFLOP/s | 3.4e-3 | Noise-tolerant surrogates, ~1e-3 fine |
Errors are relative to an fp64 gold solve. The PDE workload lands at 1.4e-6 / 3.5e-4 / 2.3e-3 on the same three rungs. The lever's 2.3e-3 to 3.4e-3 sits below the input noise of the surrogate-modeling problems this targets.
The middle rung is a finding in its own right, not just a waypoint. fp32-accumulate on the tensor cores runs about 2.8× the true-fp32 CUDA-core speed while holding error near 3.5e-4, three orders of magnitude tighter than the lever. When 1e-3 is too loose but the full fp32 path is too slow, that middle setting is the answer, and it costs almost no accuracy relative to true fp32.
Where the lever holds
The most load-bearing decision in the whole build was drawing the boundary before claiming the win. The PDE benchmark uses the screened Poisson operator, I − α∇² with α=2e-3, which is well-conditioned at cond ≈ 146. It does not use plain Poisson, whose inverse is ill-conditioned and would amplify the fp16 error until the result fell apart. That is not a favorable-case cherry-pick. It is a statement of where the lever holds, and it comes with the three places it does not.
Why screened Poisson and not plain Poisson?
Plain Poisson, ∇²u = f, has a near-singular discrete operator: its smallest eigenvalues sit close to zero, so the inverse has enormous entries and any rounding error in the apply is multiplied by a large condition number. Running that inverse through fp16 accumulation would amplify the error well past 1e-3.
The screened operator I − α∇² adds a small α=2e-3 shift that lifts the spectrum away from zero, giving cond ≈ 146. At that conditioning the fp16 error stays bounded near 2.3e-3, which is why the benchmark uses it. Choosing the well-conditioned operator up front, and saying so, is the difference between reporting where a method works and hiding where it does not.
Three conditions mark the boundary, each one routing the workload back to the fp32 path:
- Ill-conditioned operators amplify the error. A near-singular system turns 1e-3 into something much larger, because the inverse magnifies whatever the accumulator rounds off. On a stiff or near-singular operator, stay on the fp32 path. This is exactly why the screened operator was chosen over plain Poisson.
- Long time-marching drifts. Stepping a PDE forward in fp16 accumulates error along the trajectory. The screened operator’s smoothest eigenmode has eigenvalue near 1, so its per-step rounding error is never damped and compounds across steps. The lever is validated for independent, one-shot applies, not for iterative time-stepping.
- A ~1e-3 floor that refinement cannot cheaply fix. The fully-fp16 path (fp16-stored inverse plus fp16 accumulation) floors at 3.4e-3 relative error and 7.9e-3 residual. Recovering more digits is not cheap, for the reason in the aside below. If a result needs better than 1e-3, stay on the fp32 path.
Why iterative refinement doesn't cheaply buy back the last digit
The usual way to rescue a low-precision solve is iterative refinement: compute a cheap approximate solution, measure the residual, and correct. It works when the correction is cheap relative to the solve. Here it is not, for two separate reasons.
First, for a dense explicit inverse, each refinement step is itself a full GEMM, the same 2·n²·R work as the original apply. One refinement step roughly doubles the cost, which erases most of the speedup it was meant to protect.
Second, Newton–Schulz iterative inversion, X ← X(2I − AX), stalls in fp16. Its early corrections are small, and below fp16’s resolution they round to zero, so the iteration never gets moving. The workable recipes are to run the inversion in fp32, or to fp16-warm-start and then refine in fp32; a pure-fp16 refinement loop does not converge to a useful floor. This is a negative result kept because it was hard-won, not trimmed because it is inconvenient.
How every number was checked
The reason a modest 1.9× is credible where a flashy 5.3× would not be is that every figure hangs off a reproducible chain, and the chain is deliberately conservative at each step.
- fp64 gold, not fp32, as the accuracy reference. Error is measured against a double-precision solve, so the fp32 and fp16 rungs are both scored against the same ground truth.
- The fair baseline holds the kernel fixed. The headline compares like-for-like; the flattering cross-path number is shown and then demoted; the middle rung is reported too, so nothing is hidden between the extremes.
- Timing is GPU-side. Paired
torch.cuda.Eventtimestamps, one untimed warmup, median of ten repetitions, so Python dispatch overhead does not leak into the number. - TF32 is controlled.
allow_tf32 = Falsebefore the true-fp32 path, so the baseline is a genuine CUDA-core fp32 rather than a TF32-accelerated one that would flatter the lever. - The artifact is the claim. A fixed seed feeds committed
results_solve.jsonandresults_pde.json, which are the evidence, and the figure is rendered from those same JSONs.
The same lever, three domains
The last reason this matters is that it is not a physics trick. It is a property of a library default, so it fires anywhere the hot loop is a dense tensor-core GEMM whose output tolerates roughly 1e-3 of noise. The identical one-line flag carries to GPU retrieval scoring, where it is 1.59× on the scoring matmul with recall@10 holding at 0.987 against the fp32-accumulate ranking of the same embeddings. It carries to neural-field rendering, re-validated there against PSNR rather than relative error.
| Domain | Custom CUDA | Lever speedup | Quality gate |
|---|---|---|---|
| Dense solve / PDE (this work) | none, ~231 LOC bench | 1.84–1.91× | rel-err 2–3e-3 vs fp64 gold |
| Retrieval / recsys | none, ~246 LOC bench | 1.59× scoring | recall@10 = 0.987 |
| Neural-field render | 3 .cu, fusion on top | 1.69–1.81× render | PSNR 35.70 |
The retrieval and render numbers are measured in their own repos against their own metrics; they are not asserted to equal the physics 1.84–1.91×. In all three the lever is the same one-line toggle. Only the render repo added custom CUDA, and that is fusion layered on top of the flag, not the flag itself.
Why the same flag does nothing for sparse search
The retrieval repo went looking for where the lever stops, and found it in sparse workloads. Unstructured sparse matrix multiply (SpMM through cuSPARSE) reached only 0.93 TFLOP/s, roughly 120× below the dense fp16-accumulate ceiling near 111.7 TFLOP/s on the same card. Structured 2:4 sparsity through CUTLASS, the one sparse pattern that does reach the tensor cores, came in at just 1.04× over its matched dense baseline.
The reason is the arithmetic-intensity check from earlier, run in reverse. Sparse multiply is memory-bound and mostly runs on the CUDA cores, so there is no half-rate accumulator to unlock; the faster accumulator has nothing to accelerate. This is the mirror image of the physics boundaries. Each domain maps where its lever works and where it stops, and for this lever the line is dense tensor-core GEMM on one side, memory-bound sparse on the other.
What we learned, and where it goes next
On a consumer RTX 3080, the recurring cost of a dense physics solve is a compute-bound tensor-core GEMM that cuBLAS runs at half rate for a precision the workload never reads. One line switches the accumulator to full rate and returns 1.91× on the inverse apply and 1.84× on the PDE apply, at 2.3e-3 to 3.4e-3 relative error against an fp64 gold solve.
The lasting result is the method as much as the number. Crediting a speedup to one variable means measuring against the baseline that moves only that variable, which is why the figure is 1.91× over a matched tensor-core kernel and not the 5.3× a cross-path comparison would print. That discipline is what makes the number portable. It states exactly what changed and exactly what it cost, so it holds up when carried into a new setting. The boundary is part of the same result, and it reads as a routing map rather than a caveat. Ill-conditioned operators, long time-marching, and anything needing better than 1e-3 belong on the fp32 path, because for those the floor the accumulator adds is no longer below the error already in the problem. Everywhere the floor stays under the noise, the throughput is there to take.
Where this goes next is already in view. The same one-line toggle carried unchanged from a dense physics solve to retrieval scoring and neural-field rendering, each gating the trade on its own quality metric, and the pattern holds for any hot loop that reduces to a large dense tensor-core GEMM whose output tolerates roughly 1e-3, a wide class of surrogate, ranking, and rendering work. Onboarding each new domain has cost about a page of benchmark code and, so far, zero custom kernels, so the next domain is a measurement away, not a rebuild.