When the GPU wins: mandelbrot

gpu/mandelbrot renders the Mandelbrot set at 4096² with histogram equalisation. Its main.bend is byte-identical to upstream’s benchmark, and it is the mirror image of the previous chapter: the same balanced fork tree, but the work at the leaves is uniform, and that turns out to be the whole difference. The code below is that file, quoted by line number.

Why the leaves are uniform

The set is drawn by iterating z ↦ z² + c a fixed number of times per pixel and recording when z escapes. The obvious implementation breaks out of the loop when it escapes — and that is a branch, which makes pixels take different amounts of time, which is the thing that kills the GPU.

So it does not break out. Upstream’s own header says it in a line:

every pixel runs ITERS iterations with no branches (once escaped, sel freezes z rather than jumping out), so the instruction stream is identical for every pixel

The loop that does that is mit — one iteration per match step, ending when the counter runs out rather than when the point escapes:

def mit(n: Nat, +cr: U32, +ci: U32, +zr: U32, +zi: U32, esc: U32,
        it: U32) -> U32:
  match n:
    case 0n:
      it
    case 1n+p:
      +r2 = asr8(U32.mul(zr, zr))
      +i2 = asr8(U32.mul(zi, zi))
      +e2 = U32.or(esc, b2u(U32.is_gt(U32.add(r2, i2), 1024)))
      nzr = U32.add(U32.sub(r2, i2), cr)
      nzi = U32.add(asr8(U32.mul(2, U32.mul(zr, zi))), ci)
      sr = sel(e2, nzr, zr)
      si = sel(e2, nzi, zi)
      mit(p, cr, ci, sr, si, e2, U32.add(it, b2u(U32.is_zero(e2))))

Follow e2: once the escape test |z|² > 1024 has fired, esc stays set forever, because it is or-ed into itself on the next iteration. And e2 never appears in a case, so nothing branches on it. It only decides values:

def sel.go(x: U32, y: U32, z: Bool) -> U32:
  match z:
    case True{}:
      x
    case False{}:
      y

def sel(t: U32, x: U32, y: U32) -> U32:
  sel.go(x, y, U32.is_zero(t))

sel is a select — a conditional value, not a conditional jump. sr and si are the new or the old z, chosen without a jump, and a frozen z stops changing. Every pixel takes the same path through the same instructions; only the data differs.

The escape count is still correct, because the iterations after the escape are wasted work — and wasted work costs nothing to a machine built to run many lanes in lockstep. That is the exact inverse of n-queens. There, the search tree could not be made uniform. Here, it can be made uniform by spending more arithmetic.

The two passes

The render is not one fork. It is two, plus something in between, and each pass has a type of its own:

type Hs is Data:
  Hl{h0: U32, h1: U32, h2: U32, h3: U32, h4: U32, h5: U32, h6: U32, h7: U32}

type Lut is Data:
  Lut{k0: U32, k1: U32, k2: U32, k3: U32, k4: U32, k5: U32, k6: U32, k7: U32}

Hs is one histogram — eight scalar buckets, one per eighth of the iteration budget:

# escape time -> one of 8 buckets, knob-stable: (it*8)/ITERS clamped
def bkt(it: U32, tn: Nat) -> U32:
  +bq = U32.div(U32.mul(it, 8), U32.from_nat(tn))
  sel(b2u(U32.is_gt(bq, 7)), bq, 7)

1. Histogram. Fork 2^18 blocks of 64 pixels; each block sorts its escape counts into those 8 buckets and the buckets are merged pairwise up the fork tree:

def hfold(+d: Nat, +i: U32, +tn: Nat) -> Hs:
  match d:
    case 0n:
      hchunk(64n, i, tn, 0, 0, 0, 0, 0, 0, 0, 0)
    case 1n+p:
      a b = hfold(p, i, tn) hfold(p, U32.add(i, U32.shln(64, p)), tn)
      hzip(a, b)

hfold is the same depth-d fork as batch in the previous chapter, and hzip is its smerge: one bucket set added to another, pairwise, all the way up.

2. The CDF. A serial pass turning the root histogram into an equalisation lookup table. Small, sequential, unavoidable — and worth reading, because it is the one piece of this program with no fork in it at all:

# histogram equalization: LUT k = cdf(k) scaled to 0..255
def cdf(h: Hs) -> Lut:
  match h:
    case Hl{+h0, h1, h2, h3, h4, h5, h6, h7}:
      +c1 = U32.add(h0, h1)
      +c2 = U32.add(c1, h2)
      +c3 = U32.add(c2, h3)
      +c4 = U32.add(c3, h4)
      +c5 = U32.add(c4, h5)
      +c6 = U32.add(c5, h6)
      +c7 = U32.add(c6, h7)
      Lut{U32.div(U32.mul(h0, 255), c7), U32.div(U32.mul(c1, 255), c7), U32.div(U32.mul(c2, 255), c7), U32.div(U32.mul(c3, 255), c7),
          U32.div(U32.mul(c4, 255), c7), U32.div(U32.mul(c5, 255), c7), U32.div(U32.mul(c6, 255), c7), U32.div(U32.mul(c7, 255), c7)}

3. Recolour. Fork again, one leaf per pixel this time, running each pixel back through the lookup table and summing by position:

def rcol(+d: Nat, +i: U32, +tn: Nat, +k0: U32, +k1: U32, +k2: U32, +k3: U32,
         +k4: U32, +k5: U32, +k6: U32, +k7: U32) -> U32:
  match d:
    case 0n:
      rpix(i, tn, k0, k1, k2, k3, k4, k5, k6, k7)
    case 1n+p:
      a b = rcol(p, i, tn, k0, k1, k2, k3, k4, k5, k6, k7) rcol(p, U32.add(i, U32.shln(1, p)), tn, k0, k1, k2, k3, k4, k5, k6, k7)
      U32.add(a, b)

The checksum mixes the lookup table and the recoloured result, so it is sensitive to both passes. That matters — it means a GPU build that silently skipped the recolour would not produce the recorded value.

The numbers

Three runs each, checksum 3101455856 everywhere.

Wall clock says the GPU is about 6× faster than ten cores. That understates it badly, because roughly 85 ms of every GPU run is the entry cost measured in the GPU chapter, and it does not shrink when the work does.

Subtract it and the picture is:

This is the largest GPU win in the repository, and it is worth being precise about why it is not reported as 20×: because a wall-clock number for a program that runs once has to include the door. A real renderer loops over many frames or many scenes and pays it once — at which point this becomes the 20× number.

Two things to know before timing it yourself

The first run is slow. The first measurement of this binary here was 0.226 s, which settled to 0.108 s on repetition — a factor of two. Every GPU number in this book is three runs, and the first is discarded in spirit even where it is printed.

The stderr line is not an error, and you must not filter it out.

bend: compiling the GPU program (... .gpu is missing or stale)

The host binary checks for a kernel and rebuilds it when it cannot find one. If you grep -v that line away while timing, you are timing a different program. The cost is real and it belongs in the number.

Sizes

The source records two configurations:

Both are one call apart:

# small: rend!(2n, 7n) expect 887240761; big rend!(18n, its()) expect 3101455856
def main() -> IO(Unit):
  IO.print(U32.show(rend!(18n, its())))

The two binaries in this directory are the big one. The small one is useful for checking a rebuild — it finishes quickly and still verifies.

The files, and how to build them

bend main.bend -o gpu            # the GPU path, plus a .gpu Metal kernel
bend main.bend -o main.c
clang -O2 main.c -o cpu -lm      # CPU only

A useful sanity check: build the CPU binary yourself and compare. The result here was 70,696 bytes with the same checksum and the same 5.15 s as the committed cpu, which is a decent sign that the generated C is deterministic and that nothing about the measurement depends on which clang invocation you used.

That closes the performance half of the book. What is left is a worked example — Conway’s Life, built slowly enough to see every piece of the machinery, affinity, kinds and the fork included, earn its place — and then the same machinery pointed at a different question: not how fast, but how do you know it is right, without reading the code.

Next: Life the obvious way, and the trap in it.