Skip to main content

Command Palette

Search for a command to run...

Argmin in C#: from 3 to 36 billion elements per second

Updated
19 min readView as Markdown
N
A young .NET Developer who has spent half his life working with and on .NET

What porting Algorithmica's "Argmin" case study to .NET taught us about SIMD, RyuJIT and Zen 4, and why it ended in a patch for TensorPrimitives.IndexOfMin.


What this is about

Argmin is the simplest task imaginable: return the index of the smallest element in an array. In Algorithmica's Algorithms for Modern Hardware it serves as a case study in how far you can push a loop with a data-dependent branch once you take the branch predictor, SIMD and memory bandwidth seriously. The final resulting Algorithm will run at the speed of the L2 cache on my machine.

We ported the book's full progression to C# (.NET 11), three times over, once in each of the vector programming models .NET offers:

Model Namespace What it is
Vector<T> System.Numerics Width-agnostic, 256 bits on our machine, no lane permutes
Vector256<T> System.Runtime.Intrinsics Fixed 256 bits with ISA-neutral helpers (Min, ConditionalSelect, Shuffle, ...)
Avx2 / Avx / Sse2 System.Runtime.Intrinsics.X86 Raw instructions, mirroring _mm256_* calls

Everything was measured with BenchmarkDotNet on a Ryzen 7 7800X3D (Zen 4, AVX-512 capable), N = 65536 ints, i.e. 256 KB, which sits comfortably in L2. The road from the first working version to the final table was much longer than expected, and almost every detour taught something that applies beyond argmin.

The journey in short:

Stage GElem/s (Random)
Scalar loop 3.1
Index-vector loop, first version (both vectors blended with the same mask) 6.6
Index-vector loop as in the book (vpminsd + one blend) 20
Vector<T> with a stackalloc index helper 3.5
Vector<T> with Vector<int>.Indices 14.7
Vector<T> with Vector.LoadUnsafe 19.2
Branch-on-mask 24
Blocks of 32 35
Blocks of 256 36
TensorPrimitives.IndexOfMin<int> from the BCL 5.6
Our block algorithm as TensorPrimitives.IndexOfMin<int> 36

Stage 1: the scalar baseline

The naive loop is the baseline:

public static int Simple(ReadOnlySpan<int> a)
{
    int k = 0;
    for (int i = 1; i < a.Length; i++)
        if (a[i] < a[k])
            k = i;
    return k;
}

3.1 GElem/s on random data. The book's next step is a variant with [[unlikely]], which tells the compiler the branch is almost never taken and thereby steers the code layout. The JIT has no such attribute, but the order of the condition does the same job: writing the compare the other way round (> instead of <) flips which path becomes the fall-through, so the direction of the condition works as a likely hint.

The scalar loop is as fast as it is because of the branch predictor. On random data the running minimum almost never drops after the first few hundred elements, so the branch is nearly perfectly predictable and costs nothing. That observation, that the hit rate of the running minimum decides what a branch costs, runs through every branch-based variant below.


Stage 2: the index-vector loop and the bug we wanted to pin on the JIT

The book's first SIMD variant keeps two vectors: eight running minima and their eight indices. Each iteration loads a block, compares, and updates the indices with a blend.

Our first version looked the way you would write it "generically": compare, then select both the minimum and the index with the same mask. This incidentally also mirrored the implementation pattern found in the current version of IndexOfMin from the BCL TensorHelpers.

var mask = Avx2.CompareGreaterThan(min, v);   // v < min
min = Avx2.BlendVariable(min, v, mask);
idx = Avx2.BlendVariable(idx, cur, mask);

6.6 GElem/s. The portable Vector256.ConditionalSelect in the same form reached 17.3. A look at the disassembly showed that on our AVX-512 machine the JIT does not turn Avx2.BlendVariable into vpblendvb but into vpmovb2m plus vpblendmb. After some further digging we found that:

  1. vpblendvb has no EVEX encoding. On AVX-512 hardware the JIT's importer (hwintrinsicxarch.cpp, case NI_AVX2_BlendVariable) therefore rewrites every BlendVariable into the masked form vpblendmb, converting the vector mask into a k-register with vpmovb2m.
  2. If the mask is used once, lowering folds the whole thing back into a plain vpblendvb. If it is used twice, the vpmovb2m is shared by CSE and stays.
  3. On Zen 4, vpmovb2m is a vector-to-mask domain crossing with 4 cycles of latency. And because we selected the minimum through the mask as well, that crossing sat on the loop-carried chain of the minimum.

So the BCL apparently has supoptimal code written for IndexOfMin. The book advances the minimum with _mm256_min_epi32 and blends only the index:

var mask = Avx2.CompareGreaterThan(min, v);   // v < min
idx = Avx2.BlendVariable(idx, cur, mask);     // the mask's only use: stays vpblendvb
min = Avx2.Min(min, v);                       // vpminsd, a 1-cycle chain

Now the mask is single-use, the blend stays vpblendvb, and the minimum chain is a single vpminsd. All three models now compile to the same vpcmpgtd / vpblendvb / vpminsd loop and run at about 20 GElem/s.

The lesson we learned today is that on my hardware, mask is only cheap if it has exactly one use.


Stage 3: Vector<T> and two more traps

The Vector<T> variant of the same loop initially sat at 3.5 GElem/s, a sixth of the Vector256<T> version. Two causes, both at source level.

Trap 1: a stackalloc helper is never inlined

Vector<T> has no constructor from eight literals, so we had built the starting vector {0, 1, ..., W-1} in a small helper:

static Vector<int> MakeIndices()
{
    Span<int> tmp = stackalloc int[Vector<int>.Count];
    for (int i = 0; i < tmp.Length; i++) tmp[i] = i;
    return new Vector<int>(tmp);
}

RyuJIT inlines a method containing localloc only when the importer can turn the allocation into a fixed-size local. That requires the size to be a constant in the IL stack expression and at most 32 bytes (DEFAULT_MAX_LOCALLOC_TO_LOCAL_SIZE in importer.cpp). Otherwise the inline attempt is fatal (CALLSITE_LOCALLOC_SIZE_UNKNOWN), with or without AggressiveInlining.

Roslyn evaluates every non-literal stackalloc size once into an IL local, because the Span constructor needs it too. Vector<int>.Count, Vector256<int>.Count, sizeof(...): all of them arrive as ldloc, which the importer does not fold. Only a true C# compile-time literal such as stackalloc int[8] qualifies, and even that only with AggressiveInlining.

The consequence of the failed inline: the helper returns the vector through a hidden buffer, the caller's copy becomes address-exposed, and the loop-carried blend round-trips through the stack every iteration:

vmovups  ymm4, [rsp+60]
vpblendvb ymm3, ymm4, ymm2, ymm1
vmovups  [rsp+60], ymm3

The fix is Vector<int>.Indices, which the JIT materialises as a constant (gtNewSimdGetIndicesNode). That alone lifted the variant to 14.7 GElem/s.

Trap 2: new Vector<int>(span.Slice(i)) keeps two range checks

The rest of the gap was the load. new Vector<int>(a.Slice(i)) produces two bounds checks per iteration (one for Slice, one for the constructor's length check) that the JIT does not eliminate from the loop condition i + w <= n. The fixed-width variants load with Vector256.LoadUnsafe instead.

Vector.LoadUnsafe(ref MemoryMarshal.GetReference(a), i) is its width-agnostic twin. With it, every Vector<T> argmin variant compiles to the same loop as Vector256<T>, and the index-vector loop reaches 19.2 GElem/s.


Stage 4: small things you only see in the disassembly

Two details that each touch a handful of instructions but are measurable.

A nuint index with a precomputed last start. All loops were written as for (; i + w <= n; i += w) with an int index and loaded through (nuint)i. That costs a movsxd and a lea per iteration. With nuint i and last = n - w the index-vector loop shrinks from nine to seven macro-ops but does not get faster (19.9 → 20.1): it was never dispatch-bound. The variants that carry more scalar work per vector operation do gain 8 to 21 percent, e.g. BranchOnMask from 19.3 to 23.5 and Blocks256 in Vector256<T> from 33.6 to 36.3.

Rescanning over a slice instead of [from, to). The scalar rescan helper that several variants call on a hit was written as for (j = from; j < to) over the whole span. The JIT cannot prove to <= a.Length at the call sites, so a jae RNGCHKFAIL per element stayed in. Iterating over a.Slice(from, to - from) instead lets the JIT recognise the j < s.Length with s[j] pattern and drop the check. On decreasing data, where the rescan dominates, that took BranchOnMask from 24.8 to 19.3 µs.


Stage 5: why 20 and not 35?

The index-vector loop stalled at 20 GElem/s while the book reaches the sane with its block variants on a Zen2. The progression there:

BranchOnMask compares a block against the broadcast running minimum, tests the mask with vptest and only does a scalar rescan of the block on a hit. On random data an 8-block only hits when it holds a new global minimum, about 12 times per 65536 elements. The hot loop is thus six macro-ops without any loop-carried vector operation: 24 GElem/s.

Blocks32 loads four vectors, folds them with three vpminsd and does one compare per 32 elements: 35 GElem/s. That is the L2 fill throughput for a 256 KB working set. A probe loop that only loads and does vpminsd reads 36.

Blocks256 first computes the minimum of each 256-element block in vector code, remembers the winning block and then searches it for the first match by equality compare: 36 GElem/s, and unlike Blocks32 also on decreasing data, because the second phase never looks at more than 256 elements.

Why does the index-vector loop stay at 20? Not because of the select. Zen 4 has no per-instruction bottleneck in this loop: vpminsd, vpcmpgtd, vpaddd, vpternlogd and vpblendmd are all one-cycle operations on four pipes, vpblendvb one cycle on two. But the probe shows: the identical instruction sequence with the compare reading a constant instead of the minimum runs at 25.5, and unrolling by two with independent (minimum, index) pairs lifts the book loop to 24.7. Every further instruction that consumes the loaded vector costs throughput, because the core cannot keep enough iterations in flight ahead of L2. The lever for the index-vector loop is independent accumulators, not a different select instruction.


Stage 6: the generic template, and what Zen 4 has to say about it

Anyone writing an argmin without the book almost inevitably lands on the template "compare, select the minimum and the index with one mask, bump the index vector". That is exactly the form the BCL uses too (more on that below). We wanted to know whether any formulation of this template can beat the book loop, and used the probe workflow to try every select variant the JIT can emit.

Formulation (N = 65536, Random) Loop instructions GElem/s
Book: vpcmpgtd, vpblendvb (index), vpminsd (minimum) 7 20.5
Vector256.ConditionalSelect ×2 with a shared mask → vpternlogd ×2 8 18.1
Avx512F.VL.TernaryLogic for the index select → JIT emits vpcmpgtd k1 + vpblendmd 7 21.4 (needs AVX-512)
Avx512F.VL.CompareLessThan + Avx512F.VL.BlendVariable ×2 → vpcmpltd k1, vpblendmd ×2 8 8.0
Avx2.BlendVariable ×2 with a shared mask → vpmovb2m + vpblendmb ×2 8 6.6
Two separate compares, each mask single-use → vpblendvb ×2 9 13.1

The uops.info numbers explain the order: on Zen 4, a compare that writes a k-register (vpcmpgtd k, ymm, ymm) has 4 cycles of latency, against 1 cycle for the vector-to-vector variant. vpmovb2m is likewise a 4-cycle crossing. As soon as the running minimum passes through a mask register, every step costs 4 cycles. And even the best branchless form of the template, two vpternlogd, puts the compare on the minimum chain (2 instructions per step instead of one vpminsd). The book says "min is faster", and that holds here.

Two JIT quirks worth knowing in passing: Avx512F.VL.BlendVariable with a single-use k-mask is lowered back to the VEX vpblendvb, and Avx2.And/AndNot/Or as a hand-built select becomes a four-instruction vpternlogd/vpand mess.

The branchy reading beats the book, but only once tuned properly

One formulation of the template did beat the book loop: the branchy one. Compare, vptest, continue if no lane was smaller, otherwise both selects and build the index vector from the loop counter.

var v = Avx.LoadVector256(p + i);
var mask = Avx2.CompareGreaterThan(g, v);           // v < running minimum (g: in every lane)
if (Avx.TestZ(mask, mask)) continue;                // predicted: no new minimum in this block
var cur = Avx2.Add(Vector256.Create((int)i), lanes);
min = Vector256.ConditionalSelect(mask, v, min);
idx = Vector256.ConditionalSelect(mask, cur, idx);
g = BroadcastMin(v);                                // vpermq/vpshufd tree, no trip through a GPR

What makes the difference is the hit count, because every hit is a branch mispredict. Comparing against the per-lane running minimum, a block hits as soon as any of its eight lanes sees a lane record: 68 times per 65536 random elements. Comparing against the broadcast global minimum, it hits only on a new global minimum: 13 times. In the tiered probe the per-lane form read 24.1, the broadcast form 27.8 to 29.2, BranchOnMask 26.4. BenchmarkDotNet confirms it: 23.5 GElem/s on Random against 20.7 for the book loop.

The price: on decreasing data every block hits, the loop runs the full hit path per eight elements and drops to 9.3. The branchless template form is better there at 18.1, because it has no branch to pay for, and BranchOnMask with its scalar rescan drops to 3.4.

The same hit path under Blocks32

The nicest finding of this stage: put the same vectorized hit path under the Blocks32 loop (four loads, three vpminsd, one compare against the broadcast minimum, one vptest per 32 elements) and you keep Blocks32's throughput on random data (34.8 against 35.1) while replacing the 32-iteration rescan that collapses Blocks32 on decreasing data: 14.2 instead of 3.6 GElem/s.

One subtlety: the compares inside the hit block have to be against the per-lane minimum, not the old global one. Otherwise a later vector of the block can overwrite a smaller value already recorded in the same lane.


Stage 7: how you measure something like this at all

A large share of the time went not into code but into the question of how to compare formulations quickly and reliably. What proved itself:

A probe project with DOTNET_JitDisasm. A throwaway console app with one [MethodImpl(NoInlining)] method per candidate, each an exact mirror of the repo method. The release runtime is enough:

DOTNET_TieredCompilation=0 DOTNET_JitDisasm='Cand*' DOTNET_JitDisasmDiffable=1 DOTNET_JitStdOutFile=disasm.txt dotnet JitProbe.dll

Then print only the loop block per method. The authoritative check remains BenchmarkDotNet's --disasm, but the probe takes seconds instead of minutes.

Timing in the probe is only right with TieredCompilation=0, otherwise a 500-call warmup ends inside the 100 ms call-counting delay and you measure half tier-0/OSR code (the index-vector loop read 7 instead of 20 GElem/s).


Stage 7: 512 bits, and why wider is not always better

DOTNET_PreferredVectorBitWidth=512 leaves Vector<int>.Count at 8 on Zen 4, because Vector512.IsHardwareAccelerated is already true there and that variable only raises the preferred width. Vector<T> follows DOTNET_MaxVectorTBitWidth. With the right variable the index-vector loop becomes vmovups zmm, vpaddd, vpcmpltd k1, vpblendmd zmm{k1}, vpminsd and runs 29 percent faster from unchanged source (19.8 → 25.5): the blend now goes through a mask register and escapes the vpblendvb register constraint.

Blocks32, on the other hand, loses 7 percent and Blocks256 30 percent (30.7 → 21.6). Zen 4 executes 512-bit operations as two 256-bit halves, so the Min loops do not get faster, while the horizontal step through the Vector<T> indexer now walks 16 instead of 8 lanes per block, and Blocks256 does that 256 times per call.


Stage 8: TensorPrimitives.IndexOfMin and the road into the BCL

As a reference we had included TensorPrimitives.IndexOfMin<int> from System.Numerics.Tensors: 5.6 GElem/s. A sixth of Blocks256, a quarter of the index-vector loop. This is the moment where all the earlier findings show up in a single code path.

The BCL (.NET 10) uses the index-vector algorithm, but its operator keeps result and resultIndex with a first-index tie-break: LessThan, an EqualsAny branch that ORs in Equals & IndexLessThan on equality, then BlendVariable(left, right, ~mask) for both vectors. Two things follow:

  1. The minimum is produced by the blend instead of by vpminsd, so the chain is compare, invert, blend. The EqualsAny branch and the ~mask are dead weight: a strict compare already keeps the earlier index within a lane, and the horizontal step picks the lowest index among equal lanes.
  2. The mask is a vector-typed value that the JIT materialises as such. On AVX-512 hardware the chain therefore crosses between mask and vector registers twice per iteration. The 512-bit loop has 15 instructions, five of them on the result chain, and takes about 15 cycles per 16 elements.

The proof that it is the mask crossings: with DOTNET_EnableAVX512=0 the same code produces the plain AVX2 loop with vpblendvb ×2, a three-instruction chain, and runs at 11.1 GElem/s, twice as fast. Same binary, AVX-512 merely switched off.

In dotnet/runtime main (PR #127454, .NET 11) the core shared by IndexOfMin, IndexOfMax, IndexOfMinMagnitude and IndexOfMaxMagnitude had already been rewritten to the strict index-vector template, but it still blends the minimum with the same mask as the index: 5.1 GElem/s here, no faster than our implementation.

The port

We wrote our Vector256_Blocks256 against the BCL's operator interface and dropped it in as the shared IndexOfMinMaxCore<T, TOperator>, for all four searches:

  • Pass 1 reduces each block of 32 vectors with two accumulators of the operator's lane-wise Reduce (Vector.Min, Vector.Max, Vector.MinMagnitude, Vector.MaxMagnitude, the same aggregation operators the old core already relied on for its final match) and keeps the first block whose result wins under the operator's strict scalar Compare.
  • Pass 2 scans only that block for the first element the result does not beat under Compare. That is the operator's own tie rule (first equal element, -0 before +0 for min, +1 over -1 for max-magnitude) without any bitwise equality.
  • NaN falls out of the reductions being IEEE minimum/maximum: a block with a NaN reduces to NaN, then the first NaN of that block is returned.

IndexOfMin<int> now runs at 35.6 to 36.7 GElem/s on all four input patterns at 512 bits, identical to Vector256_Blocks256. The block loop is two vpminsd zmm with memory operands, the only call is the final scan.

A first cut was 12 percent slower: two inlined block reductions plus the scan helpers had exhausted the JIT's inlining budget, and it emitted a call HorizontalAggregate per block. One block loop and NoInlining on the two helpers that run once per search fixed that.

The runtime's own tests pass for all four searches and all 17 element types on the 512-, 256- and 128-bit and scalar paths, extended with tests for multi-block lengths, NaN and signed zero. The PR is up as dotnet/runtime#133969.


What remains

The findings that carry beyond argmin:

  1. The book is right about min instead of blend. The running minimum belongs on a 1-cycle vpminsd chain. Every formulation that puts the compare on that chain loses, and every one that routes it through a mask register loses four cycles per step on Zen 4.
  2. Avx2.BlendVariable is only cheap on AVX-512 hardware with a single-use mask. A shared mask becomes vpmovb2m. Vector256.ConditionalSelect with a shared mask becomes vpternlogd and is the better choice when you genuinely need two selects.
  3. Vector<T> needs Vector<T>.Indices and Vector.LoadUnsafe. A stackalloc helper with a non-literal size is never inlined, and new Vector<T>(span.Slice(i)) keeps two range checks. With both fixes Vector<T> is on par for vertical arithmetic; for anything involving lane movement it remains unusable.
  4. Bounds checks stay wherever the JIT cannot prove the length. for (j = from; j < to) over a span keeps the check; a.Slice(from, to - from) with j < s.Length does not.
  5. The index-vector loop is L2-bound, not select-bound. The lever is independent accumulator pairs, not the select instruction. The block variants reach the bandwidth because they execute only one vpminsd per loaded vector.
  6. Hit count equals mispredicts. For every branchy variant, what matters is how often the running minimum drops. Between "never" and "ln N" hits the difference is invisible; all separation between variants happens between Random and Decreasing. Comparing against the global instead of the per-lane minimum cuts the hit count from 68 to 13 per 65536 elements.

Appendix: the final table

BenchmarkDotNet --job short, N = 65536 ints, Ryzen 7 7800X3D, .NET 10. GElem/s, higher is better.

Variant Increasing HitOneOverN Random Decreasing
Scalar_Simple 3.0 3.1 3.1 3.1
Vector256_IndexVector 19.8 20.6 20.3 19.9
Avx2_IndexVectorBlendBoth 17.1 18.1 17.6 17.5
Avx2_IndexVectorBlendBothOnHit 22.0 23.9 22.1 8.8
Vector256_BranchOnMask 22.3 24.2 24.2 3.3
Vector256_Blocks32 35.0 35.2 33.3 2.2
Avx2_Blocks32BlendBothOnHit 34.7 35.7 34.8 14.2
Vector256_MinThenFind 35.3 32.4 16.9 15.8
Vector256_Blocks256 34.8 36.6 35.9 35.3
TensorPrimitives_IndexOfMin (BCL, .NET 10) 4.9 5.4 5.6 5.6
TensorPrimitives_IndexOfMin (our port) 36 36 36 36

The input patterns span the hit rate: Increasing (a[i] = i, the minimum is the first element and never moves), HitOneOverN (increasing baseline, block k holds a new global minimum with probability 1/k, 9 block hits per 65536 elements), Random (Random(42).Next() per element, 12 block hits and 51 lane records) and Decreasing (every element is a new minimum). All variants return the first index of the minimum.

MinThenFind is the only variant whose cost depends not on the hit count but on the position of the minimum (0, 7077, 59735, 65535 in the four patterns): its second pass scans up to the first occurrence.

The code lives in https://github.com/Mrnikbobjeff/HPC, the full tables and disassembly findings are in the repo README.

Algorithmic Improvements

Part 1 of 1

Working on improving the .NET BCL