← Writing
kernelsJune 17, 2026 · 18 min read

One cuBLAS flag makes GPU retrieval scoring 1.59× faster — and nothing for sparse search

The full path to 1.59× on a consumer RTX 3080: every retrieval and recsys pattern collapses to one dense matmul that cuBLAS runs at half rate, one flag flips it, and recall@10 holds at 0.987 against the fp32 ranking. Two numbers reported, 1.59× on the matmul and 1.25× end to end, plus the sparse follow-on that was measured and refuted on purpose.

1.59×faster retrieval scoringfp32-accfp16-acc

Strip a search engine or a recommender down to its hot loop and one operation is left: score a batch of query embeddings against a corpus of item embeddings, then keep the top few per query. The scoring is a single matrix multiply, Q [n_q×d] @ Cᵀ [d×n_c]. On a consumer RTX 3080 that matmul runs at half the card’s tensor-core rate by default, for a numerical safety margin retrieval never asked for. One line of PyTorch takes the other half back.

Two numbers come out of it, not one, because the fair report of a pipeline is the whole pass and not its fastest slice: 1.59× on the scoring matmul, 1.25× once the top-k merge is included. The ranking check holds at recall@10 0.987 against the fp32 ranking of the same embeddings. Then the obvious next conjecture, that the same flag would also speed up sparse search, was measured and killed on purpose. This is the reasoning in order, the dead-end included.

1.59×
faster scoring matmul
91.3 vs 57.5 TFLOP/s
1.25×
faster end to end
485 vs 608 ms, incl. top-k
0.987
recall@10 vs fp32 ranking
9.87 of 10 results survive

The whole stack collapses to one matmul

The starting move is not to optimize a search index. It is to notice that retrieval and recommendation come in what look like three different architectures, and all three bottom out in the same operation.

  • Brute-force maximum inner-product search scores every query against every item and takes the exact top-k.
  • Two-tower recommenders run a user tower and an item tower, then score user vectors against item vectors.
  • ANN-then-rerank narrows to a candidate shortlist with an approximate index, then scores that shortlist exactly.

Different names, one kernel. Each ends in a dense general matrix multiply, Q [n_q×d] @ Cᵀ [d×n_c], followed by a top-k selection. There is no fourth case to special-case, and that reduction is what makes one lever general. Optimize the GEMM and every retrieval pattern inherits the gain at once.

Three architectures funnel into one dense GEMM; a separate top-k stage runs downstream.
Brute-force MIPSevery query against every item, exact
Two-tower recsysuser vectors against item vectors
ANN + rerankquery against a candidate shortlist, exact
Q·Cᵀ dense GEMMtensor cores · the flag acts here

The pipeline is a dense GEMM plus a separate top-k, and that second half matters as much as the first. Top-k is not a matrix multiply. It does not touch the tensor cores, and nothing done to the accumulator can speed it up. Recognizing that split before running anything is why the benchmark commits in advance to reporting two numbers, the matmul alone and the whole pass. Quoting only the first would flatter the lever; the gap between them measures how much top-k overhead a real pipeline carries.

Why all three architectures reduce to the same GEMM

Brute-force maximum inner-product search forms the full score matrix Q @ Cᵀ and reads off the largest entries per row. A two-tower recommender is the same product, with the query rows produced by a user tower and the columns by an item tower. ANN-then-rerank shrinks n_c from the whole corpus to a shortlist first, but the exact rerank on that shortlist is still Q @ Cᵀ at a smaller width.

In every case the arithmetic is the dot products, 2 · n_q · n_c · d floating-point operations, and the selection that follows is a comparison, not a multiply. The accumulator lives inside the dot products, so that is the only place a flag on the accumulator can act. The shortlist width changes how much work there is; it does not change what the work is.

Why the default leaves half the card idle

With the target fixed on the scoring matmul, the next question is why it 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 products into an accumulator that can be fp16 or fp32. On consumer Ampere (sm_86, the RTX 3080’s architecture) 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 on every input, because it is the numerically safe choice on any input. Datacenter Ampere carries no such cap, so this is specifically a GeForce ceiling: anyone serving retrieval on consumer cards is leaving roughly half the tensor throughput on the floor without knowing it.

That default is priced for the worst input, and retrieval is not the worst input. This is the deduction the whole lever rests on, and it is a property of the workload, not of the hardware. The consumed output of retrieval is a ranking, not raw scores. A user sees an ordered list; nobody reads the third decimal of a dot product. So a little numerical noise is tolerable as long as the ordering survives. fp16 accumulation trades on the order of 1e-3 relative per-score error (an estimate from the format’s ten-bit mantissa, not a number this benchmark records) for throughput. For retrieval, whose inputs are already lossy embeddings sitting behind a quantized index, that is well inside the noise the system already tolerates. The trade is only conditionally safe, which is exactly why it forces the verification step that comes later.

The lever is a runtime toggle, flipped right before the scoring matmul and off right after:

torch.backends.cuda.matmul.allow_fp16_accumulation = True   # PyTorch >= 2.7
scores = queries @ corpus.T        # now runs at the full tensor-core rate
torch.backends.cuda.matmul.allow_fp16_accumulation = False

No custom CUDA, no rebuild, no change to the model or the embeddings. The work here is not the kernel; it is knowing which workloads can take the lower-precision path and proving retrieval is one of them. That proof is the rest of this post.

What is the accumulator, and why does fp32 cost double on sm_86?

A tensor-core mma instruction multiplies fp16 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, a segmentation choice in the consumer part. cuBLAS defaults to fp32 accumulation because it is correct on any input, including matrices with a large dynamic range where fp16 accumulation would overflow or drop small terms. Retrieval scores are bounded dot products of unit-scale embeddings, so neither failure mode applies, and the safe default is simply leaving performance unused. The lever holds the fp16 inputs fixed and moves only the accumulator, so nothing about the data or the kernel shape changes; only the internal summation precision does.

What the flag buys, measured two ways

The benchmark is a realistic scoring pass: a 2.1M-item corpus of fp16 embeddings at d=512, 8,192 queries, top-10 per query, the corpus streamed through the GPU in 64 MiB tiles of 65,536 items so memory stays flat no matter how large the corpus grows. Two measurements matter, and the split is the same one drawn at the start: the scoring matmul on its own, which is exactly what the flag touches, and the end-to-end pass including the top-k merge, which is what a pipeline actually pays.

Scoring matmul only · 2.1M-item corpus, d=512, 8,192 queries · RTX 3080 · median of 5
fp32-accumulate (cuBLAS default)57.5 TFLOP/s · 306 ms
fp16-accumulate (the flag)91.3 TFLOP/s · 193 ms · 1.59×
End-to-end pass incl. top-10 merge · same corpus and hardware · median of 5
fp32-accumulate28.3 G-scores/s · 608 ms · computed
fp16-accumulate35.4 G-scores/s · 485 ms · 1.25×

The scoring matmul goes 1.59× faster. End to end the win shrinks to 1.25×, because selecting the top 10 per query is not a matmul. It runs on the CUDA cores at the same speed regardless of the accumulator, and it takes a fixed slice of wall-clock that dilutes the lever. The 1.25× is the number to quote for a full pipeline, and the gap between it and 1.59× points straight at the next optimization: fuse the top-k into the scoring kernel so more of the pass rides the faster path.

RTX 3080 · torch 2.12 / CUDA 13 · median of 5 runs
PathScoring matmulThroughputEnd-to-endSpeedup
fp32-accumulate (default)306 ms57.5 TFLOP/s608 ms1.00×
fp16-accumulate (the flag)193 ms91.3 TFLOP/s485 ms1.59× / 1.25×

Exact factors: 1.588 scoring, 1.254 end to end. Scoring-only is faster because the top-k merge costs the same either way and dilutes the end-to-end number.

Streaming a 2 GiB corpus, and why the scoring lever is 1.59× and not 2×

The corpus is 2,097,152 items at d=512 in fp16, exactly 2 GiB. It would fit in the RTX 3080’s 10 GB, but the benchmark streams it through the GPU in 64 MiB tiles of 65,536 items, 32 tiles per pass. That keeps the resident footprint flat and constant regardless of corpus size, so the same loop scales to corpora far larger than VRAM. Each tile is scored against all 8,192 queries, and for the end-to-end measurement each tile’s local top-10 is merged into a running buffer of shape (8192, 10).

That streaming is also why the scoring lever lands at 1.59× rather than at the roughly 2:1 the hardware issue rate would allow. The pass moves a real 2 GiB of corpus through HBM and emits scores, so it is not perfectly compute-bound; the accumulator flag speeds up the arithmetic but not that memory traffic. The 2× is a hardware ceiling for a purely compute-bound GEMM. 1.59× is what a corpus-streaming scoring pass actually reaches.

Does the ranking survive?

Half-rate accumulation changes the scores, so the real question is whether it changes the ranking. The answer has to isolate a single variable, and the way to do that is to compare against the ranking the same fp16 embeddings produce under fp32 accumulation. Same corpus, same queries, same code path; the only thing that differs between the two lists is allow_fp16_accumulation.

recall@10 came in at 0.987. On average 9.87 of each query’s true top-10 results still appear after the switch. For approximate nearest-neighbor retrieval, whose inputs are already lossy embeddings behind a quantized index, that is well inside the noise the system already tolerates.

Why recall, and not raw per-score error

The consumed output of retrieval is an ordered list, so the metric that matters is how much of the list survives, not how far individual scores move. Per-score error under fp16 accumulation is on the order of 1e-3 relative, an estimate from the format’s ten-bit mantissa rather than a recorded measurement. That error only matters where it flips an ordering, and near the top of a ranked list the gaps between neighboring scores are usually larger than the noise.

recall@10 = 0.987158203125 says the flips are rare: 9.87 of every 10 top results are unchanged. A score error that never crosses a ranking boundary is invisible to the user, and this is the check that turns the conditional licensing argument (noise is fine if the ranking holds) into a measured fact rather than a hope.

The obvious next conjecture is that if one flag unlocks the dense half of a retrieval stack, it should also help the sparse half. Lexical matching, co-occurrence tables, and graph adjacency are all sparse-matrix multiplies, and they sit right next to the dense scoring in a real system. The reasoning treated this as a claim to test and possibly kill, not an assumption to carry forward, so it got its own benchmark with a dense-GEMM reference ceiling built in. It was measured, and it was refuted.

Sparse hypothesis · n=8192, 1% density · RTX 3080 · median of 10
Dense fp16-acc GEMM111.7 TFLOP/s · square 8192³ · tensor cores
Unstructured SpMM, cuSPARSE0.93 TFLOP/s · 120× below · CUDA cores

Unstructured sparse multiply lands 120× below the dense tensor-core ceiling (119.8× exactly) because it is memory-bound and dispatches to cuSPARSE on the general-purpose CUDA cores. It never reaches a tensor core, so there is no half-rate accumulator to unlock. The flag touches a machine this workload never runs on. That is the same arithmetic-intensity argument that made the dense lever work, run in reverse.

One sparse pattern does reach the tensor cores: 2:4 structured sparsity, where exactly two of every four values are nonzero, run through CUTLASS. On this shape it came in at 1.04× over its matched dense baseline, with fp16 accumulation already on for both sides so the number reflects the sparsity pattern alone. The literature quotes 1.3–1.7× for 2:4 in principle; 1.04× is what this shape measured on a 3080, and the gap between principle and measurement is the finding, not a footnote.

2:4 structured vs matched dense baseline · 8192×8192×4096 · fp16-acc both sides · median of 10
Dense baseline (fp16-acc)114.6 TFLOP/s · 4.80 ms
2:4 structured (CUTLASS)119.2 TFLOP/s · 1.04× · marginal
Sparse operations · RTX 3080 · median of 10 runs
OperationResultvs dense tensor-core ceiling
Dense fp16-acc GEMM (square 8192³)111.7 TFLOP/s1.0× (the ceiling)
Unstructured SpMM (sparse × dense), cuSPARSE0.93 TFLOP/s120× slower
SpGEMM (sparse × sparse), cuSPARSE14.8 msno FLOP rate defined
2:4 structured (sparse × dense), CUTLASS1.04× vs matched densereaches tensor cores

Unstructured paths run on the CUDA cores and are memory-bound, a different machine than the tensor cores the flag touches. SpGEMM output has 37,274,926 nonzeros.

Why unstructured SpMM never touches a tensor core

Tensor cores consume dense tiles. An unstructured sparse matrix at 1% density has no dense tiles to feed them, so cuSPARSE runs the multiply on the general-purpose CUDA cores as a gather-and-accumulate over irregular indices. That work is memory-bound: it moves index and value arrays and does little arithmetic per byte, sitting far below the roofline ridge where the tensor cores live.

The accumulator width the flag controls belongs to a unit this kernel never dispatches to, so there is nothing for it to accelerate. At 0.93 TFLOP/s the SpMM is 120× under the 111.7 TFLOP/s dense ceiling, which is the arithmetic-intensity gap made concrete: dense scoring is compute-bound on the tensor cores, unstructured sparse is memory-bound on the CUDA cores, and no accumulator flag moves work from one machine to the other.

2:4 structured sparsity, and why its 1.04× is the pattern alone

2:4 structured sparsity keeps exactly two nonzeros in every group of four, a layout the Ampere tensor cores can accept directly through CUTLASS. To isolate the effect of the pattern, fp16 accumulation is turned on for both the dense baseline and the 2:4 side, so the accumulator is held fixed and only the sparsity changes.

The matched dense baseline runs at 114.6 TFLOP/s (4.80 ms); the 2:4 kernel reaches 119.2 TFLOP/s, a 1.04× gain. The literature cites 1.3–1.7× for 2:4 in principle, so this shape on a 3080 underperforms the theory, which is why it is reported as marginal rather than as a win. The isolation matters: because fp16 accumulation is on for both sides, the 1.04× is attributable to the structured pattern and not to the accumulator question at all.

Why SpGEMM has no throughput number

Sparse-times-sparse multiply (SpGEMM) has no well-defined FLOP rate, because the number of multiply-adds depends on the sparsity pattern of the product, which is not known until it is computed. Fabricating a TFLOP/s figure would mean dividing by an operation count that is itself an artifact of the inputs.

The benchmark therefore records wall-time and output size rather than a made-up rate: 14.8 ms, with 37,274,926 nonzeros in the result. Like SpMM it runs on the CUDA cores and never reaches a tensor core, so the accumulator flag has nothing to act on here either.

Turn the refutation into a design

The refutation is a design, not a dead end. Once the boundary is drawn between the memory-bound sparse machine and the compute-bound dense machine, the pipeline that follows falls out directly. Keep candidate generation and re-ranking as separate stages, and apply the flag only where it lands.

120×
unstructured SpMM below dense
0.93 vs 111.7 TFLOP/s
1.04×
2:4 structured, this shape
literature claims 1.3–1.7×
two-stage
the design it prescribes
sparse gen, then dense re-rank
The flag applies only to the dense re-rank stage, and that is by design.
query
sparse candidate-genCUDA cores · flag: no effect
dense re-ranktensor cores · flag: 1.59×
top-k
results
  • Sparse candidate generation (inverted index, graph walk) is memory-bound on the CUDA cores. The flag does nothing here, and that is fine. This stage was never on the tensor cores.
  • Dense re-ranking (score the shortlist against embeddings) is the GEMM. This is where the 1.59× lives.

The end-to-end 1.25× says the same thing from the other direction: the top-k merge is now the tail worth chasing. Fuse it into the scoring kernel and the end-to-end number climbs from 1.25× back toward the 1.59× the matmul already gets.

How every number was checked

A one-line flag is easy to oversell, so every figure hangs off a discipline that is deliberately conservative:

  • recall@k isolation. Ground truth for recall is the fp32-accumulate ranking of the same fp16 embeddings, not a separate higher-precision model or human relevance labels. That leaves the accumulator as the only variable, so recall@10 measures the accumulator’s effect and nothing else.
  • Two numbers, not the flattering one. The scoring-only 1.59× and the end-to-end 1.25× are both reported, because the pass a pipeline runs includes the top-k merge the flag cannot touch. Quoting only 1.59× would overstate what a system gains.
  • A reference ceiling for the negative result. The sparse conjecture was given a purpose-built benchmark with a dense-GEMM ceiling in it, so “sparse is slow” is a measured 120× ratio, not a vibe, and the 2:4 test holds fp16 accumulation on for both sides so its 1.04× is the pattern alone.
  • GPU-side timing. Every figure is timed with paired CUDA events, median of 5 runs for retrieval and 10 for the sparse tests, so Python dispatch overhead does not leak into the number.
  • The artifact is the claim. The numbers here come from committed results_retrieval.json and results_sparse.json produced by the two benchmark scripts; the tables and ladders read from those same files.

What the numbers cover, and what’s next

Every number is one RTX 3080 (sm_86, consumer Ampere, 10 GB), torch 2.12 with CUDA 13. The claims hold inside that boundary. Two of the questions it raises, how close the lever can get to the hardware ceiling and whether a wider embedding gains more, are the next measurements to run.

The ~2× ceiling is not the 1.59× result

The hardware fp16-to-fp32 accumulate issue-rate ratio is about 2:1, which sets a theoretical ceiling near 2× for a perfectly compute-bound GEMM. That ceiling is asserted from the hardware, and it is not what any benchmark here measured.

The measured scoring lever is 1.59×, and it sits below the ceiling because a real scoring pass streams a 2 GiB corpus through HBM and emits scores, spending time on memory traffic and tiled-loop tails that the ideal ratio ignores. The 2× is the roofline; 1.59× is what this workload reached, and reading the ceiling as the result would overstate the lever. The distance between them is not slack in the claim but headroom with a direction. Anything that makes the pass more compute-bound, fusing the top-k or widening the embedding, moves the measured number toward the roofline, and those are the next two levers.

Would a wider embedding dimension gain more? Untested.

There is a plausible argument that a larger d would raise the speedup, because wider embeddings mean more arithmetic per byte of score traffic, which pushes the matmul further into the compute-bound regime where the accumulator rate dominates. It is the same mechanism that explains why 1.59× sits below 2×, run forward instead of backward.

Every retrieval number here is d=512, so the artifact holds one point on that curve, not the two a trend needs. The dimension sweep is the clean next experiment. Hold the corpus and the flag fixed, vary d, and watch whether the speedup climbs toward the roofline the way the compute-bound argument predicts. Wider embeddings are the hypothesis to test next, and it is a one-variable test.

The same flag, other domains

This is one line, not a kernel, and it is a property of a library default rather than a retrieval trick. It fires anywhere the hot loop is a dense tensor-core GEMM whose consumed output tolerates roughly 1e-3 of noise, at a fixed cost each domain re-validates with its own quality metric.

One lever, more than one domain · each speedup and quality gate measured in its own repo
DomainCustom CUDALever speedupQuality gate
Retrieval / recsys (this work)none, runtime toggle1.59× scoring, 1.25× end to endrecall@10 = 0.987
Dense solve / PDEnone, runtime toggle1.84–1.91× applyrel-err 2–3e-3 vs fp64 gold
Neural-field renderfusion on top of the flag1.69–1.81× renderPSNR 35.70

The physics and render numbers are measured in their own repos against their own metrics; they are not asserted to equal the retrieval 1.59×.

The identical flag carries to dense physics solves, where it measures 1.91× on a matrix-inverse apply and 1.84× on a PDE Green’s-function apply at about 1e-3 relative error against a double-precision reference, and to neural-field rendering, re-validated there against PSNR. One default, switched off by libraries for safety, is worth roughly 1.6× to 1.9× wherever the hot loop is a dense tensor-core GEMM whose output tolerates about 1e-3 of noise.

1 line
of PyTorch
no custom kernel, no rebuild
0.987
recall@10 quality gate
ranking survives the switch
sm_86
consumer Ampere only
A100 not measured here

On this hardware, measured against the ranking you would actually ship, one flag buys 1.59× on the scoring matmul and a true 1.25× on the pass that includes the top-k the flag cannot touch. The same benchmark that proved it also mapped what comes next. Unstructured sparse never reaches the tensor cores, so the pipeline splits into two stages and the lever rides only the dense re-rank. The distance from 1.25× to 1.59×, then from 1.59× toward the hardware’s 2×, is headroom with a known path. Fuse the top-k into the scoring kernel, then sweep the embedding width the compute-bound argument says should gain more.