Computing Ext to stem 400 on an HPC cluster: a worklog

• Row reduction. The classical bottleneck. Nassau's algorithm shrinks the matrices, but they are still the biggest objects we build.
• Milnor multiplication. The bottleneck since sseq adopted Nassau's algorithm in 2022: 93% of the time at stem 120.

In 2022 I asked Christian Nassau about GPUs, and guessed the gains were “potentially several orders of magnitude”. This year I worked on both at once.

-- forward pass over F₂, b columns at a time
for each panel of b columns, left to right:
  for each column j of the panel:   -- serial
    pick a row below the pivots with a 1 in j
    none? j is a free column: skip it
    swap it up; for each row r below it with a 1 in j:
      row r ← row r + pivot row;  L[r, j] ← 1
  -- k ≤ b pivots
        ┌ P  X ┐   P: the k pivot rows
        └ 0  Y ┘   Y: everything else

  Y ← Y + L · X                    -- one product
then back-substitute, right to left

And NVIDIA had its own reason to build a 1-bit matrix multiplier.

Neural networks quantize their weights and inputs to make matrix products cheap, all the way down to one bit.

Quantizing to 0/1 is reducing mod 2: our use case exactly.

  cvt.u64.u32  %rd32, %r2663;
  or.b64   %rd13, %rd32, 4611686293305360384;
  mad.lo.s32  %r2664, %r3606, 24576, %r2657;
  and.b32   %r2665, %r2664, 262016;
  shr.u32  %r2666, %r2665, 4;
  cvt.u64.u32  %rd33, %r2666;
  or.b64   %rd8, %rd33, 4611686293305360384;
  wgmma.mma_async.sync.aligned.m64n128k256.row.col.s32.b1.b1.and.popc {%r3395,%r3394,%r3393,%r3392,%r3391,%r3390,%r3389,%r3388,%r3387,%r3386,%r3385,%r3384,%r3383,%r3382,%r3381,%r3380,%r3379,%r3378,%r3377,%r3376,%r3375,%r3374,%r3373,%r3372,%r3371,%r3370,%r3369,%r3368,%r3367,%r3366,%r3365,%r3364,%r3363,%r3362,%r3361,%r3360,%r3359,%r3358,%r3357,%r3356,%r3355,%r3354,%r3353,%r3352,%r3351,%r3350,%r3349,%r3348,%r3347,%r3346,%r3345,%r3344,%r3343,%r3342,%r3341,%r3340,%r3339,%r3338,%r3337,%r3336,%r3335,%r3334,%r3333,%r3332},%rd8,%rd13, 1;
  add.s32  %r2667, %r2664, 8192;
  and.b32   %r2668, %r2667, 262016;
  shr.u32  %r2669, %r2668, 4;
  cvt.u64.u32  %rd34, %r2669;
  or.b64   %rd10, %rd34, 4611686293305360384;

wgmma

warpgroup MMA: 128 threads issue it together

mma_async

runs while the threads do other work

sync.aligned

every thread arrives at the same instruction

m64n128k256

a $64 \times 128$ output tile, $256$ bits deep

row.col

operand layouts

s32

32-bit counts

b1.b1

one-bit inputs

and.popc

AND, then count the ones

{%r3395, …}

64 registers per thread: 128 threads hold the $64 \times 128$ tile

One instruction: $2 \cdot 64 \cdot 128 \cdot 256 \approx 4.2$ million bit operations.

One bit operation every 100 attoseconds, on average. In that time light travels about 30 nanometres: a hundred-odd silicon atoms.

In September's production run, row reduction was under 0.5% of the CPU time. ✓

• Nassau. Compute each term's position from a few small tables of dimensions of $\mathbb{F}_2[\xi_1, \dots, \xi_k]$.
• sseq. A hash map: $O(1)$, and faster than the tables at small degrees.
• The GPU port. The kernel needed a lookup it computes, not a hash map, so the port used Nassau's tables anyway. PR #270 brought them to the CPU.

At degree 400 the tables look up 2.7× faster, and the index on a stem-400 worker shrinks from 12.4 GB of hash maps to 2.1 GB of tables. Below degree ~300 the hash map wins: it falls out of cache only as the basis grows.

Call for call, one H200 does the work of 79 server cores. Multiplication: good enough, for now ✓

The reference for everything after it.

• Order the signatures colexicographically, $\sigma_1 < \sigma_2 < \cdots < \sigma_N$. Step $i$ reads the $\sigma_i$-coordinates of one shared vector, the differential $dx$ being fixed, solves for a correction $c_i$, and writes it back: $dx \leftarrow dx + d(c_i)$.
• $d(c_i)$ has coordinates only in signatures $\sigma_j \ge \sigma_i$. So step $j$ reads what steps $i < j$ wrote: a forward substitution, not $N$ independent solves.
• Checked at stem 120: 33% of the entries a step reads had been changed by earlier steps (up to 44%), in 86% of bidegrees.

The signatures are a filtration, not a grading: corrections only go up. So which signatures feed which?

And it grows: $A(n)$ has $2^{(n+1)(n+2)/2}$ signatures, but its DAG is only $\binom{n+3}{3}$ deep. The width is exponential, the depth is a cubic: the further out we go, the more concurrency there is.

Need to scale horizontally.

Another node sees only what we send it.

So we need self-contained, well-defined tasks: pieces of work we can hand over whole.

We have them: the signature steps of Act III, each answered by a correction we add back.

Is the network fast enough? Who hands out the tasks? And in what order?

$\dfrac{m}{\beta} \ll t$

$m$ bytes over a link of bandwidth $\beta$, against $t$, the time to compute the task.

Our H200 nodes' link is 25 Gb/s Ethernet: $\beta \approx 3$ GB/s.

Both comfortably. The row reductions are about 97% of the bytes. ✓

• The coordinator holds the resolution and the bidegree and signature DAGs, sends tasks, and merges the answers by XOR. It never computes.
• A task is one signature step: a bidegree, a signature, and each generator's $d(x)$ there: a few hundred KB. Its correction comes back: a few tens of MB.
• A worker keeps its own copy of the resolution: the coordinator sends every new generator to every worker.
• An rref server is an H200 that only row-reduces. Workers send it their large matrices; cards without Hopper send all of them.

On the same four GPUs this is no faster. It is what lets any GPU on the cluster join.

• The problem. $Q \mid \mathit{prec} \mid C_{\max}$: machines of different speeds, tasks in a DAG, finish as early as possible.
• Plus memory. A card runs several tasks at once. Card $i$ has memory $M_i$, task $j$ needs $m_j$, and at every moment $\sum_{j \text{ on } i} m_j \le M_i$.
• The easy case. Identical machines, no memory (Graham, 1966): any schedule that never idles finishes by $W/P + D$, at most twice the lower bound $\max(W/P, D)$.
• Ours. Speeds differ by card and by kind of task ($R$, not $Q$), and every time is an estimate. A greedy scheduler can hand the critical task to the slowest card.

A GPU runs 32 threads in lockstep: a thread that hits a dead end idles until its warp is done. Building only matrices that contribute keeps the warps busy on useful work, and needs no table of matrices in memory.

The GPU kept its lead: one H200 went from ≈ 79 server cores to ≈ 250–300 (estimates). Milnor ✓

• Kernels are compiled at run time, for whatever card the worker lands on.
• Only Hopper has the warpgroup 1-bit instruction our kernels use: the others send their large row reductions to an H200 over the network.
• Under 2 GB of device memory is what lets a 16 GB V100 or a 24 GB A30 join at all.
• Each new kind of card first redoes stems 0–200 and must match bit for bit: 0 mismatches.

The CPU code stays the reference for every new kernel.

$A(5)$ only becomes usable at stem 585, so the signature DAGs stay at most 36 deep. At 512 the obstacle is the volume of work, not the dependencies.