WebGPU · WGSL · compute shaders

Thread by Thread

An introduction to GPU programming in 17 exercises. Every example runs live on your GPU: edit the code on the right, and the picture on the left recompiles as you type.

By Zoravur Singh and Claude (Anthropic) · Source on GitHub

Start here How this page works

Each exercise on this page is a small live GPU program. The picture on the left is drawn by the code on the right, and the code runs on your graphics card, once per pixel, every frame.

The strip under each exercise lists what the GPU runs every frame: each entry point, in order, with how many workgroups it launches. You’ll learn what that means in a minute.

The playground follows the conventions of compute.toys, a public WebGPU shader playground. Almost everything you write here can be pasted into compute.toys and will run unchanged. When you outgrow this page, that’s the place to go next.

Why this feels hard

If shaders feel harder than they look, that’s normal. The syntax of WGSL, the WebGPU Shading Language, looks like a small, stripped-down Rust or C. But the language isn’t the hard part. The execution model is:

People who say shaders are easy have usually internalized the parallel mental model and forgotten how long it took. The exercises below build that model one idea at a time.

Part 0 The mental model

One thread per pixel

A compute shader is a function that the GPU runs many times at once. Each running copy is called an invocation (or a thread). You don’t call it in a loop. Instead you dispatch it: you tell the GPU how many copies to launch, and it runs them in parallel.

Invocations are launched in fixed-size batches called workgroups. The line

@compute @workgroup_size(16, 16)

says each workgroup is a 16 × 16 square of 256 threads. To cover a 640 × 480 screen, the playground launches 40 × 30 workgroups. Each thread learns its position through built-in inputs:

Built-in Type Meaning
global_invocation_id vec3u This thread’s position in the whole dispatch. For a screen, it’s the pixel.
local_invocation_id vec3u Position inside its own workgroup, from (0, 0) to (15, 15).
workgroup_id vec3u Which workgroup this thread belongs to.

They’re related by one formula:

the dispatch: 7 × 4 workgroupsscreenedge workgroups hang past the screen,which is why every shader starts with a bounds checkworkgroup (3, 1): 16 × 16 threads, 8 × 8 drawnthread with local_invocation_id (2, 1)global_invocation_id = workgroup_id × 16 + local_invocation_id → (3, 1) × 16 + (2, 1) = (50, 17)
Each workgroup covers a 16 × 16 tile of pixels. The global ID is the workgroup’s corner plus the thread’s offset inside it.

For the first half of this page, only global_invocation_id matters: it tells each thread which pixel it owns. The local and workgroup IDs become important in Part VI, when threads in a workgroup start cooperating.

Anatomy of the starting shader

Here is the program every exercise starts from. It’s the same default shader compute.toys gives you. Read through it, then try changing the numbers on the col line.

@compute @workgroup_size(16, 16)
fn main_image(@builtin(global_invocation_id) id: vec3u) {
    // Viewport resolution (in pixels)
    let screen_size = textureDimensions(screen);

    // Prevent overdraw for workgroups on the edge of the viewport
    if (id.x >= screen_size.x || id.y >= screen_size.y) { return; }

    // Pixel coordinates (centre of pixel, origin at bottom left)
    let fragCoord = vec2f(f32(id.x) + .5, f32(screen_size.y - id.y) - .5);

    // Normalised pixel coordinates (from 0 to 1)
    let uv = fragCoord / vec2f(screen_size);

    // Time varying pixel colour
    var col = .5 + .5 * cos(time.elapsed + uv.xyx + vec3f(0., 2., 4.));

    // Convert from gamma-encoded to linear colour space
    col = pow(col, vec3f(2.2));

    // Output to screen (linear colour space)
    textureStore(screen, id.xy, vec4f(col, 1.));
}

Line by line:

fn main_image(@builtin(global_invocation_id) id: vec3u) declares the entry point. @builtin(global_invocation_id) asks the GPU to fill in id with this thread’s global position. vec3u is a vector of three unsigned 32-bit integers. Only .x and .y matter for a 2D image.

textureDimensions(screen) returns the size of the output image in pixels. screen is a storage texture the playground provides. Each thread writes its pixel into it with textureStore.

The bounds check exists because the dispatch is rounded up to whole workgroups. A screen 650 pixels wide needs 41 workgroups of 16, which is 656 threads per row, so six threads per row have no pixel. Returning early keeps them from writing outside the image.

fragCoord rebuilds the pixel coordinate from the thread ID. The name comes from Shadertoy’s GLSL, where fragment shaders get fragCoord for free. Three things happen in that line:

  1. id.x and id.y say which pixel this thread owns.
  2. The .5 moves to the center of the pixel. Pixel 0 covers the interval from 0 to 1, so its center is 0.5. Fragment shaders use the same convention, so this avoids half-pixel errors later.
  3. screen_size.y - id.y flips y. In textures, row 0 is at the top and y grows downward. In math and in Shadertoy, y grows upward. After the flip, the top row gets fragCoord.y = H - 0.5 and the bottom row gets 0.5.

The f32(...) conversions are needed because WGSL never mixes integers and floats silently. You always convert explicitly.

uv divides by the screen size, turning pixel units into the range 0 to 1, with (0, 0) at the bottom left. uv means the same thing at any resolution, which is why most effects are written in terms of it.

col is the color. uv.xyx is a swizzle: it builds a vec3f from the components x, y, x. Adding a scalar to a vector adds it to every component.

pow(col, vec3f(2.2)) converts the color into linear light. The screen texture stores linear values, and the display converts them back. If you want a value of 0.5 to look like middle grey, apply this line. Every exercise keeps it.

textureStore(screen, id.xy, vec4f(col, 1.)) writes the pixel. The fourth component is alpha. Notice that storing uses id.xy, the raw texture coordinate with y pointing down, not fragCoord.

What the playground gives you

These names are always available. They mirror the compute.toys prelude.

Name What it is
screen The output image, texture_storage_2d<rgba16float, write>. Write with textureStore(screen, coord, color).
time.elapsed Seconds since start (f32).
time.frame Frame counter (u32). It’s 0 on the first frame after a restart.
time.delta Seconds since the last frame.
mouse.pos Mouse position in pixels (vec2i), with the origin at the top left like id.
mouse.click 1 while the button is held on the canvas, else 0.
passLoad, passStore Read and write scratch images between entry points. See Part V.
#storage name type Declares a storage buffer, such as #storage data array<f32>. See Part VI.
#workgroup_count name x y z Sets how many workgroups an entry point launches.
#dispatch_count name n Runs an entry point n times in a row each frame. dispatch.id tells you which run it is.

Every @compute function in your file is an entry point, and all of them run each frame in the order they appear. If you don’t give one a #workgroup_count, it gets enough workgroups to cover the screen.

Part I Coordinates and color

1 Gradient

A vector constructor can take a mix of smaller vectors and scalars, as long as the component count adds up. vec3f(uv, 0.0) takes the two components of uv and appends a zero:

let a = vec3f(uv, 0.0);          // (uv.x, uv.y, 0.0)
let b = vec3f(0.2);              // (0.2, 0.2, 0.2): one value fills every slot
let c = vec4f(a.zy, 1.0, 1.0);   // swizzles work anywhere

Since uv runs from 0 to 1 across the screen, using it as a color is the simplest way to see a coordinate system. This is the first debugging tool you’ll use over and over: when you’re not sure what a value is doing, draw it.

Make red increase from left to right and green increase from bottom to top. The bottom-left corner should be black and the top-right corner yellow. If the green is brightest at the bottom, your y is flipped.

@compute @workgroup_size(16, 16)
fn main_image(@builtin(global_invocation_id) id: vec3u) {
    let screen_size = textureDimensions(screen);
    if (id.x >= screen_size.x || id.y >= screen_size.y) { return; }

    let fragCoord = vec2f(f32(id.x) + .5, f32(screen_size.y - id.y) - .5);
    let uv = fragCoord / vec2f(screen_size);

    // TODO: red from uv.x, green from uv.y, no blue
    var col = vec3f(0.0);

    col = pow(col, vec3f(2.2));
    textureStore(screen, id.xy, vec4f(col, 1.));
}
@compute @workgroup_size(16, 16)
fn main_image(@builtin(global_invocation_id) id: vec3u) {
    let screen_size = textureDimensions(screen);
    if (id.x >= screen_size.x || id.y >= screen_size.y) { return; }

    let fragCoord = vec2f(f32(id.x) + .5, f32(screen_size.y - id.y) - .5);
    let uv = fragCoord / vec2f(screen_size);

    // uv.x drives red, uv.y drives green
    var col = vec3f(uv, 0.0);

    col = pow(col, vec3f(2.2));
    textureStore(screen, id.xy, vec4f(col, 1.));
}

Try removing the pow line. The gradient gets brighter and the midpoint looks washed out. That’s the difference between linear light and the way screens encode brightness.

2 Hard split

On a CPU, you’d write “if x is past the middle, white; otherwise black.” On a GPU, you usually write it as math instead.

The reason is how GPUs execute. Threads run in lockstep groups (often 32 at a time) that share one instruction stream. When threads in the same group take different sides of an if, the group runs both sides and masks out the results each thread shouldn’t keep. This is called divergence. A small if costs almost nothing, but the habit of expressing choices as arithmetic keeps code fast and branch-free. WGSL gives you tools for it:

Function Result
step(edge, x) 0.0 if x < edge, else 1.0
select(f, t, cond) t if cond is true, else f. Note the order: false value first.
mix(a, b, t) Linear blend a + (b - a) * t
clamp(x, lo, hi) x limited to the range

Make the left half of the screen black and the right half white, using step instead of if. Then try moving the split with the mouse: mouse.pos.x is in pixels, so divide by the screen width to get the same 0-to-1 scale as uv.

@compute @workgroup_size(16, 16)
fn main_image(@builtin(global_invocation_id) id: vec3u) {
    let screen_size = textureDimensions(screen);
    if (id.x >= screen_size.x || id.y >= screen_size.y) { return; }

    let fragCoord = vec2f(f32(id.x) + .5, f32(screen_size.y - id.y) - .5);
    let uv = fragCoord / vec2f(screen_size);

    // TODO: 0.0 on the left half, 1.0 on the right half, without `if`
    let v = 0.0;
    var col = vec3f(v);

    col = pow(col, vec3f(2.2));
    textureStore(screen, id.xy, vec4f(col, 1.));
}
@compute @workgroup_size(16, 16)
fn main_image(@builtin(global_invocation_id) id: vec3u) {
    let screen_size = textureDimensions(screen);
    if (id.x >= screen_size.x || id.y >= screen_size.y) { return; }

    let fragCoord = vec2f(f32(id.x) + .5, f32(screen_size.y - id.y) - .5);
    let uv = fragCoord / vec2f(screen_size);

    // Split at the middle, or at the mouse once it has moved over the canvas
    var edge = 0.5;
    if (mouse.pos.x > 0) { edge = f32(mouse.pos.x) / f32(screen_size.x); }

    let v = step(edge, uv.x);
    var col = vec3f(v);

    col = pow(col, vec3f(2.2));
    textureStore(screen, id.xy, vec4f(col, 1.));
}

3 Checkerboard

Most patterns come from tiling space: split the plane into cells, then decide what happens in each one.

Multiply uv by 8 and you get coordinates that run from 0 to 8. floor of that is a cell index from 0 to 7 on each axis. For a checkerboard, a cell is white when the sum of its x and y index is odd. The % operator works on floats in WGSL, so (cell.x + cell.y) % 2.0 is either 0 or 1.

Draw an 8 × 8 checkerboard. Then try fract(uv * 8.0) as a color to see the coordinates inside each cell.

@compute @workgroup_size(16, 16)
fn main_image(@builtin(global_invocation_id) id: vec3u) {
    let screen_size = textureDimensions(screen);
    if (id.x >= screen_size.x || id.y >= screen_size.y) { return; }

    let fragCoord = vec2f(f32(id.x) + .5, f32(screen_size.y - id.y) - .5);
    let uv = fragCoord / vec2f(screen_size);

    // TODO: cell index with floor(), then even/odd with %
    let cell = uv * 8.0;
    var col = vec3f(cell / 8.0, 0.0);

    col = pow(col, vec3f(2.2));
    textureStore(screen, id.xy, vec4f(col, 1.));
}
@compute @workgroup_size(16, 16)
fn main_image(@builtin(global_invocation_id) id: vec3u) {
    let screen_size = textureDimensions(screen);
    if (id.x >= screen_size.x || id.y >= screen_size.y) { return; }

    let fragCoord = vec2f(f32(id.x) + .5, f32(screen_size.y - id.y) - .5);
    let uv = fragCoord / vec2f(screen_size);

    let cell = floor(uv * 8.0);            // which cell, 0..7 on each axis
    let odd = (cell.x + cell.y) % 2.0;     // 0.0 or 1.0
    var col = vec3f(odd);

    col = pow(col, vec3f(2.2));
    textureStore(screen, id.xy, vec4f(col, 1.));
}

The cells aren’t square, because uv stretches to fit the screen’s shape. The next part fixes that.

Part II Shapes with distance

4 Circle

Neither fragCoord nor uv puts the origin in the middle of the screen. You may have met the range from −1 to 1 as normalized device coordinates, but that’s the space a vertex shader outputs in the graphics pipeline, and a compute shader skips that pipeline entirely. So shader code builds its own centered space, almost always like this:

let p = (2.0 * fragCoord - vec2f(screen_size)) / f32(screen_size.y);

This puts (0, 0) at the center with y running from −1 to 1. The key detail is dividing by the height only. That keeps x and y on the same scale, so a circle stays round instead of stretching into an ellipse. On a 16:9 screen, x runs from about −1.78 to 1.78.

With a centered space, a circle is one line: a point is inside a circle of radius r when length(p) < r. Better still, think of length(p) - r as a signed distance: negative inside the circle, zero on its edge, positive outside, and in each case equal to how far you are from the edge. Distance functions like this are the core of 2D shader art.

Draw a white circle of radius 0.5 in the center. Resize the page (or rotate your phone) and check that it stays round.

@compute @workgroup_size(16, 16)
fn main_image(@builtin(global_invocation_id) id: vec3u) {
    let screen_size = textureDimensions(screen);
    if (id.x >= screen_size.x || id.y >= screen_size.y) { return; }

    let fragCoord = vec2f(f32(id.x) + .5, f32(screen_size.y - id.y) - .5);

    // Centered, aspect-correct: (0,0) in the middle, y from -1 to 1
    let p = (2.0 * fragCoord - vec2f(screen_size)) / f32(screen_size.y);

    // TODO: white where length(p) < 0.5
    var col = vec3f(abs(p), 0.0);

    col = pow(col, vec3f(2.2));
    textureStore(screen, id.xy, vec4f(col, 1.));
}
@compute @workgroup_size(16, 16)
fn main_image(@builtin(global_invocation_id) id: vec3u) {
    let screen_size = textureDimensions(screen);
    if (id.x >= screen_size.x || id.y >= screen_size.y) { return; }

    let fragCoord = vec2f(f32(id.x) + .5, f32(screen_size.y - id.y) - .5);
    let p = (2.0 * fragCoord - vec2f(screen_size)) / f32(screen_size.y);

    let d = length(p) - 0.5;               // signed distance to the circle
    var col = vec3f(select(0.0, 1.0, d < 0.0));

    col = pow(col, vec3f(2.2));
    textureStore(screen, id.xy, vec4f(col, 1.));
}

5 Soft circle

Zoom into the edge of that circle and you’ll see stair steps. Each pixel is either fully in or fully out, a problem called aliasing. The fix is to fade across the edge over about one pixel.

smoothstep(a, b, x) returns 0 when x ≤ a, 1 when x ≥ b, and a smooth S-curve in between. Feed it the signed distance and fade across a band one pixel wide on each side of the edge.

How wide is one pixel in p units? p.y spans 2 units over the screen’s height, so one pixel is 2.0 / f32(screen_size.y). (Fragment shaders can compute this with fwidth, but compute shaders have no neighboring pixels to compare against, so you work it out from the resolution.)

let px = 2.0 / f32(screen_size.y);
let inside = 1.0 - smoothstep(-px, px, d);   // 1 inside, 0 outside, soft edge

Once you have a mask from 0 to 1, mix(background, foreground, mask) paints with it. This mask-and-mix pattern is how nearly all 2D shader art is layered.

Give the circle an anti-aliased edge and paint it in a color of your choice on a dark background. Then make px ten times wider and see what happens.

@compute @workgroup_size(16, 16)
fn main_image(@builtin(global_invocation_id) id: vec3u) {
    let screen_size = textureDimensions(screen);
    if (id.x >= screen_size.x || id.y >= screen_size.y) { return; }

    let fragCoord = vec2f(f32(id.x) + .5, f32(screen_size.y - id.y) - .5);
    let p = (2.0 * fragCoord - vec2f(screen_size)) / f32(screen_size.y);

    let d = length(p) - 0.5;

    // TODO: replace the hard edge with smoothstep over about one pixel,
    //       then mix a background and a foreground color with the mask
    let inside = select(0.0, 1.0, d < 0.0);
    var col = vec3f(inside);

    col = pow(col, vec3f(2.2));
    textureStore(screen, id.xy, vec4f(col, 1.));
}
@compute @workgroup_size(16, 16)
fn main_image(@builtin(global_invocation_id) id: vec3u) {
    let screen_size = textureDimensions(screen);
    if (id.x >= screen_size.x || id.y >= screen_size.y) { return; }

    let fragCoord = vec2f(f32(id.x) + .5, f32(screen_size.y - id.y) - .5);
    let p = (2.0 * fragCoord - vec2f(screen_size)) / f32(screen_size.y);

    let d = length(p) - 0.5;
    let px = 2.0 / f32(screen_size.y);            // one pixel, in p units
    let inside = 1.0 - smoothstep(-px, px, d);

    let bg = vec3f(0.05, 0.06, 0.1);
    let fg = vec3f(1.0, 0.75, 0.3);
    var col = mix(bg, fg, inside);

    col = pow(col, vec3f(2.2));
    textureStore(screen, id.xy, vec4f(col, 1.));
}

6 Rings and unions

Because a signed distance is a real distance, you can combine shapes with simple arithmetic:

Operation Code Why it works
Union min(a, b) You’re inside either shape if either distance is negative.
Intersection max(a, b) Inside both only if both are negative.
Subtraction max(a, -b) Inside a and outside b.
Outline (ring) abs(d) - w Distance to the edge itself, thickened by w.
Move sd(p - c) Evaluate the shape in coordinates centered at c.

The last row is worth pausing on. To move a shape, you don’t move the shape: you shift the coordinates you evaluate it in. That idea returns in Part III.

The starter code draws a single circle’s distance field as contour lines, so you can see that d really is a distance everywhere on screen, not only near the edge.

Write a helper fn sdCircle(p: vec2f, r: f32) -> f32. Draw two overlapping circles as one merged shape using min, and a thin ring around both. Keep the edges anti-aliased.

@compute @workgroup_size(16, 16)
fn main_image(@builtin(global_invocation_id) id: vec3u) {
    let screen_size = textureDimensions(screen);
    if (id.x >= screen_size.x || id.y >= screen_size.y) { return; }

    let fragCoord = vec2f(f32(id.x) + .5, f32(screen_size.y - id.y) - .5);
    let p = (2.0 * fragCoord - vec2f(screen_size)) / f32(screen_size.y);

    // The distance field of one circle, drawn as contour lines.
    // Blue outside, orange inside, white on the edge.
    let d = length(p) - 0.5;
    var col = select(vec3f(0.3, 0.55, 0.95), vec3f(0.95, 0.6, 0.25), d < 0.0);
    col *= 0.75 + 0.25 * cos(d * 60.0);
    col = mix(col, vec3f(1.0), 1.0 - smoothstep(0.0, 0.015, abs(d)));

    // TODO: an sdCircle helper, two circles merged with min(), and a ring

    col = pow(col, vec3f(2.2));
    textureStore(screen, id.xy, vec4f(col, 1.));
}
fn sdCircle(p: vec2f, r: f32) -> f32 {
    return length(p) - r;
}

@compute @workgroup_size(16, 16)
fn main_image(@builtin(global_invocation_id) id: vec3u) {
    let screen_size = textureDimensions(screen);
    if (id.x >= screen_size.x || id.y >= screen_size.y) { return; }

    let fragCoord = vec2f(f32(id.x) + .5, f32(screen_size.y - id.y) - .5);
    let p = (2.0 * fragCoord - vec2f(screen_size)) / f32(screen_size.y);
    let px = 2.0 / f32(screen_size.y);

    // Two circles, each evaluated in shifted coordinates, merged with min
    let a = sdCircle(p - vec2f(-0.3, 0.0), 0.4);
    let b = sdCircle(p - vec2f( 0.3, 0.0), 0.4);
    let blobs = min(a, b);

    // A ring: distance to a circle's edge, thickened by 0.015
    let ring = abs(sdCircle(p, 0.85)) - 0.015;

    var col = vec3f(0.05, 0.06, 0.1);
    col = mix(col, vec3f(0.3, 0.55, 0.95), 1.0 - smoothstep(-px, px, blobs));
    col = mix(col, vec3f(1.0), 1.0 - smoothstep(-px, px, ring));

    col = pow(col, vec3f(2.2));
    textureStore(screen, id.xy, vec4f(col, 1.));
}

Part III Time and motion

7 Moving circle

time.elapsed is the number of seconds since the shader started. Anything that depends on it animates, because every frame the whole image is recomputed with a slightly larger value.

To move a shape, shift the coordinates, as in the last table. A center that follows (sin t, cos t) traces a circle, since sine and cosine of the same angle are the coordinates of a point on the unit circle.

time.frame counts frames instead. Prefer elapsed for motion, since frame rates vary between machines, and use frame when you need something to happen exactly once or to seed randomness.

Make a circle of radius 0.2 orbit the center at radius 0.5. Then add a second circle that orbits twice as fast in the other direction.

fn sdCircle(p: vec2f, r: f32) -> f32 {
    return length(p) - r;
}

@compute @workgroup_size(16, 16)
fn main_image(@builtin(global_invocation_id) id: vec3u) {
    let screen_size = textureDimensions(screen);
    if (id.x >= screen_size.x || id.y >= screen_size.y) { return; }

    let fragCoord = vec2f(f32(id.x) + .5, f32(screen_size.y - id.y) - .5);
    let p = (2.0 * fragCoord - vec2f(screen_size)) / f32(screen_size.y);
    let px = 2.0 / f32(screen_size.y);

    // TODO: make the center move with time.elapsed
    let center = vec2f(0.5, 0.0);
    let d = sdCircle(p - center, 0.2);

    var col = vec3f(0.05, 0.06, 0.1);
    col = mix(col, vec3f(1.0, 0.75, 0.3), 1.0 - smoothstep(-px, px, d));

    col = pow(col, vec3f(2.2));
    textureStore(screen, id.xy, vec4f(col, 1.));
}
fn sdCircle(p: vec2f, r: f32) -> f32 {
    return length(p) - r;
}

@compute @workgroup_size(16, 16)
fn main_image(@builtin(global_invocation_id) id: vec3u) {
    let screen_size = textureDimensions(screen);
    if (id.x >= screen_size.x || id.y >= screen_size.y) { return; }

    let fragCoord = vec2f(f32(id.x) + .5, f32(screen_size.y - id.y) - .5);
    let p = (2.0 * fragCoord - vec2f(screen_size)) / f32(screen_size.y);
    let px = 2.0 / f32(screen_size.y);
    let t = time.elapsed;

    let c1 = 0.5 * vec2f(sin(t), cos(t));                  // orbit clockwise
    let c2 = 0.5 * vec2f(sin(-2.0 * t), cos(-2.0 * t));    // twice as fast, reversed
    let d1 = sdCircle(p - c1, 0.2);
    let d2 = sdCircle(p - c2, 0.12);

    var col = vec3f(0.05, 0.06, 0.1);
    col = mix(col, vec3f(1.0, 0.75, 0.3), 1.0 - smoothstep(-px, px, d1));
    col = mix(col, vec3f(0.35, 0.6, 1.0), 1.0 - smoothstep(-px, px, d2));

    col = pow(col, vec3f(2.2));
    textureStore(screen, id.xy, vec4f(col, 1.));
}

8 Ripples

A wave is a sine of position minus time. sin(k * x - w * t) is a pattern with k controlling how tightly the crests are packed and w controlling how fast they move. Because the crest sits where k x - w t is constant, it travels at speed w / k.

Use distance from the center as the position and you get concentric rings moving outward. sin returns −1 to 1, so remap it with 0.5 + 0.5 * sin(...) before using it as a blend factor.

Color each pixel by sin(length(p) * 20.0 - time.elapsed * 3.0), remapped to 0–1, and blend between two colors of your choice. Fade the waves out toward the edges. Then flip the sign of the time term.

@compute @workgroup_size(16, 16)
fn main_image(@builtin(global_invocation_id) id: vec3u) {
    let screen_size = textureDimensions(screen);
    if (id.x >= screen_size.x || id.y >= screen_size.y) { return; }

    let fragCoord = vec2f(f32(id.x) + .5, f32(screen_size.y - id.y) - .5);
    let p = (2.0 * fragCoord - vec2f(screen_size)) / f32(screen_size.y);

    let d = length(p);

    // TODO: a wave that depends on d and time.elapsed, remapped to 0..1
    let wave = 0.0;
    var col = mix(vec3f(0.1, 0.2, 0.5), vec3f(0.9, 0.95, 1.0), wave);

    col = pow(col, vec3f(2.2));
    textureStore(screen, id.xy, vec4f(col, 1.));
}
@compute @workgroup_size(16, 16)
fn main_image(@builtin(global_invocation_id) id: vec3u) {
    let screen_size = textureDimensions(screen);
    if (id.x >= screen_size.x || id.y >= screen_size.y) { return; }

    let fragCoord = vec2f(f32(id.x) + .5, f32(screen_size.y - id.y) - .5);
    let p = (2.0 * fragCoord - vec2f(screen_size)) / f32(screen_size.y);

    let d = length(p);
    let wave = 0.5 + 0.5 * sin(d * 20.0 - time.elapsed * 3.0);
    var col = mix(vec3f(0.1, 0.2, 0.5), vec3f(0.9, 0.95, 1.0), wave);
    col *= smoothstep(1.5, 0.0, d);        // fade out toward the edges

    col = pow(col, vec3f(2.2));
    textureStore(screen, id.xy, vec4f(col, 1.));
}

9 Rotation

Rotating a point by angle a around the origin is a matrix multiply:

fn rot(a: f32) -> mat2x2f {
    let c = cos(a);
    let s = sin(a);
    return mat2x2f(c, s, -s, c);
}

WGSL matrices are column-major: mat2x2f(c, s, -s, c) has first column (c, s) and second column (−s, c). Multiplying rot(a) * p turns p counter-clockwise by a.

As with moving, you rotate a shape by rotating the coordinates the other way. To spin a square by +t, evaluate the square at rot(-t) * p. This is the most common source of “why is it spinning backwards?”

For the square itself, the starter includes the standard box distance function. You don’t need to derive it, but it works by measuring how far abs(p) sticks out past the box’s corner.

Spin the square around the center. Then add a second, smaller square that orbits the first while spinning around its own center. You’ll need to translate first, then rotate.

// Signed distance to an axis-aligned box with half-size b
fn sdBox(p: vec2f, b: vec2f) -> f32 {
    let q = abs(p) - b;
    return length(max(q, vec2f(0.0))) + min(max(q.x, q.y), 0.0);
}

// TODO: fn rot(a: f32) -> mat2x2f

@compute @workgroup_size(16, 16)
fn main_image(@builtin(global_invocation_id) id: vec3u) {
    let screen_size = textureDimensions(screen);
    if (id.x >= screen_size.x || id.y >= screen_size.y) { return; }

    let fragCoord = vec2f(f32(id.x) + .5, f32(screen_size.y - id.y) - .5);
    let p = (2.0 * fragCoord - vec2f(screen_size)) / f32(screen_size.y);
    let px = 2.0 / f32(screen_size.y);

    // TODO: rotate the coordinates before measuring the box
    let q = p;
    let d = sdBox(q, vec2f(0.35));

    var col = vec3f(0.05, 0.06, 0.1);
    col = mix(col, vec3f(0.35, 0.6, 1.0), 1.0 - smoothstep(-px, px, d));

    col = pow(col, vec3f(2.2));
    textureStore(screen, id.xy, vec4f(col, 1.));
}
fn sdBox(p: vec2f, b: vec2f) -> f32 {
    let q = abs(p) - b;
    return length(max(q, vec2f(0.0))) + min(max(q.x, q.y), 0.0);
}

// Counter-clockwise rotation. Columns are (c, s) and (-s, c).
fn rot(a: f32) -> mat2x2f {
    let c = cos(a);
    let s = sin(a);
    return mat2x2f(c, s, -s, c);
}

@compute @workgroup_size(16, 16)
fn main_image(@builtin(global_invocation_id) id: vec3u) {
    let screen_size = textureDimensions(screen);
    if (id.x >= screen_size.x || id.y >= screen_size.y) { return; }

    let fragCoord = vec2f(f32(id.x) + .5, f32(screen_size.y - id.y) - .5);
    let p = (2.0 * fragCoord - vec2f(screen_size)) / f32(screen_size.y);
    let px = 2.0 / f32(screen_size.y);
    let t = time.elapsed;

    // Big square: rotate the space by -t so the square turns by +t
    let d1 = sdBox(rot(-t) * p, vec2f(0.3));

    // Small square: move to its orbit position, then spin about its own center
    let center = rot(t * 0.5) * vec2f(0.75, 0.0);
    let d2 = sdBox(rot(-3.0 * t) * (p - center), vec2f(0.1));

    var col = vec3f(0.05, 0.06, 0.1);
    col = mix(col, vec3f(0.35, 0.6, 1.0), 1.0 - smoothstep(-px, px, d1));
    col = mix(col, vec3f(1.0, 0.75, 0.3), 1.0 - smoothstep(-px, px, d2));

    col = pow(col, vec3f(2.2));
    textureStore(screen, id.xy, vec4f(col, 1.));
}

Part IV Randomness

10 Static

WGSL has no random() function, and neither does any other shading language. A shader invocation has no hidden state that survives between calls, and thousands of threads can’t share one global generator without stepping on each other.

Instead you use a hash function: a function that scrambles an integer so thoroughly that its outputs look random. Feed in something unique to this thread, such as the pixel coordinate, and you get a number that looks unrelated to its neighbors. Feed in the same input and you get the same output, which is useful: randomness in shaders is reproducible by default.

A good, fast choice is the PCG hash:

fn pcg(v: u32) -> u32 {
    let state = v * 747796405u + 2891336453u;
    let word = ((state >> ((state >> 28u) + 4u)) ^ state) * 277803737u;
    return (word >> 22u) ^ word;
}

// A float in [0, 1). Keeping the top 24 bits fits exactly in an f32.
fn u2f(h: u32) -> f32 {
    return f32(h >> 8u) / 16777216.0;
}

The u suffix marks unsigned integer literals. Integer arithmetic on u32 wraps around on overflow, which is exactly what a hash wants.

To combine several inputs, nest the hash: pcg(x ^ pcg(y ^ pcg(frame))). Each layer mixes the previous result with a new input. To get several independent numbers for one pixel, keep hashing: h, pcg(h), pcg(pcg(h)).

If your seed is a float, convert it with bitcast<u32>(x), which reinterprets the float’s raw bits. Don’t use u32(x): that’s a numeric conversion that truncates, so 0.2, 0.7 and 0.99 all become 0 and give identical “random” values.

Fill the screen with colored TV static that changes every frame. Then make it hold still by removing the frame from the seed.

fn pcg(v: u32) -> u32 {
    let state = v * 747796405u + 2891336453u;
    let word = ((state >> ((state >> 28u) + 4u)) ^ state) * 277803737u;
    return (word >> 22u) ^ word;
}

fn u2f(h: u32) -> f32 {
    return f32(h >> 8u) / 16777216.0;
}

@compute @workgroup_size(16, 16)
fn main_image(@builtin(global_invocation_id) id: vec3u) {
    let screen_size = textureDimensions(screen);
    if (id.x >= screen_size.x || id.y >= screen_size.y) { return; }

    // TODO: hash id.x, id.y and time.frame together, then take three
    //       numbers from the hash for red, green and blue
    var col = vec3f(u2f(pcg(id.x)));

    col = pow(col, vec3f(2.2));
    textureStore(screen, id.xy, vec4f(col, 1.));
}
fn pcg(v: u32) -> u32 {
    let state = v * 747796405u + 2891336453u;
    let word = ((state >> ((state >> 28u) + 4u)) ^ state) * 277803737u;
    return (word >> 22u) ^ word;
}

fn u2f(h: u32) -> f32 {
    return f32(h >> 8u) / 16777216.0;
}

@compute @workgroup_size(16, 16)
fn main_image(@builtin(global_invocation_id) id: vec3u) {
    let screen_size = textureDimensions(screen);
    if (id.x >= screen_size.x || id.y >= screen_size.y) { return; }

    // Nest the hash to mix in each input. Drop time.frame for still noise.
    let h = pcg(id.x ^ pcg(id.y ^ pcg(time.frame)));
    let r = pcg(h);
    let g = pcg(r);
    let b = pcg(g);
    var col = vec3f(u2f(r), u2f(g), u2f(b));

    col = pow(col, vec3f(2.2));
    textureStore(screen, id.xy, vec4f(col, 1.));
}

The solution chains pcg to get three numbers. Why not pcg(h), pcg(h + 1u), pcg(h + 2u)? Because the pixel to the right probably has a seed one larger, so its red would equal your green. Chaining avoids that kind of accidental correlation.

11 Random tiles

Here is the key trick behind procedural textures: hash the cell, not the pixel. Every pixel in a cell computes the same cell index, so every pixel gets the same random value, and the cell appears as one solid random color.

Cell indices from floor are floats. Convert them to integers first with vec2i(...). Negative indices are fine, because bitcast<u32> of a negative i32 is just a different bit pattern.

Split the screen into 8 × 8 cells like the checkerboard, and give each cell its own random color. Then make the pattern reshuffle once per second by mixing u32(time.elapsed) into the hash.

fn pcg(v: u32) -> u32 {
    let state = v * 747796405u + 2891336453u;
    let word = ((state >> ((state >> 28u) + 4u)) ^ state) * 277803737u;
    return (word >> 22u) ^ word;
}

fn u2f(h: u32) -> f32 {
    return f32(h >> 8u) / 16777216.0;
}

@compute @workgroup_size(16, 16)
fn main_image(@builtin(global_invocation_id) id: vec3u) {
    let screen_size = textureDimensions(screen);
    if (id.x >= screen_size.x || id.y >= screen_size.y) { return; }

    let fragCoord = vec2f(f32(id.x) + .5, f32(screen_size.y - id.y) - .5);
    let uv = fragCoord / vec2f(screen_size);

    let cell = vec2i(floor(uv * 8.0));

    // TODO: hash the cell (not the pixel) into a color
    var col = vec3f(vec2f(cell) / 8.0, 0.0);

    col = pow(col, vec3f(2.2));
    textureStore(screen, id.xy, vec4f(col, 1.));
}
fn pcg(v: u32) -> u32 {
    let state = v * 747796405u + 2891336453u;
    let word = ((state >> ((state >> 28u) + 4u)) ^ state) * 277803737u;
    return (word >> 22u) ^ word;
}

fn u2f(h: u32) -> f32 {
    return f32(h >> 8u) / 16777216.0;
}

fn hash2(c: vec2i, seed: u32) -> u32 {
    return pcg(bitcast<u32>(c.x) ^ pcg(bitcast<u32>(c.y) ^ pcg(seed)));
}

@compute @workgroup_size(16, 16)
fn main_image(@builtin(global_invocation_id) id: vec3u) {
    let screen_size = textureDimensions(screen);
    if (id.x >= screen_size.x || id.y >= screen_size.y) { return; }

    let fragCoord = vec2f(f32(id.x) + .5, f32(screen_size.y - id.y) - .5);
    let uv = fragCoord / vec2f(screen_size);

    let cell = vec2i(floor(uv * 8.0));
    let h = hash2(cell, u32(time.elapsed));     // new pattern every second
    var col = vec3f(u2f(h), u2f(pcg(h)), u2f(pcg(pcg(h))));

    col = pow(col, vec3f(2.2));
    textureStore(screen, id.xy, vec4f(col, 1.));
}

12 Value noise

Random tiles are blocky. For clouds, terrain or fire you want randomness that varies smoothly. Value noise gets it by putting a random value at every corner of the grid and blending between the four corners around each point:

  1. i = floor(x) is the cell. Its corners are i, i + (1,0), i + (0,1) and i + (1,1). Hash each one.
  2. f = fract(x) is the position inside the cell.
  3. Blend the bottom two corners by f.x, the top two by f.x, then blend those two results by f.y.

Blending with f directly leaves visible creases at cell borders, because the slope jumps there. Using the curve u = f * f * (3.0 - 2.0 * f) instead (the same curve smoothstep uses) flattens the slope at each border and hides the grid.

A single layer of noise looks like soft blobs. For natural detail, add several layers, each at double the frequency and half the strength. This is called fractal Brownian motion (fBm), and it’s how the header of this page is drawn.

Write noise(x: vec2f) -> f32 and display noise(p * 4.0). Then write fbm with five octaves and animate it by adding time to the coordinates.

fn pcg(v: u32) -> u32 {
    let state = v * 747796405u + 2891336453u;
    let word = ((state >> ((state >> 28u) + 4u)) ^ state) * 277803737u;
    return (word >> 22u) ^ word;
}

// A random value in [0, 1) for each integer grid point
fn hash(c: vec2i) -> f32 {
    return f32(pcg(bitcast<u32>(c.x) ^ pcg(bitcast<u32>(c.y))) >> 8u) / 16777216.0;
}

fn noise(x: vec2f) -> f32 {
    let i = vec2i(floor(x));
    let f = fract(x);
    // TODO: hash the four corners of cell i and blend them smoothly
    return hash(i);
}

@compute @workgroup_size(16, 16)
fn main_image(@builtin(global_invocation_id) id: vec3u) {
    let screen_size = textureDimensions(screen);
    if (id.x >= screen_size.x || id.y >= screen_size.y) { return; }

    let fragCoord = vec2f(f32(id.x) + .5, f32(screen_size.y - id.y) - .5);
    let p = (2.0 * fragCoord - vec2f(screen_size)) / f32(screen_size.y);

    var col = vec3f(noise(p * 4.0));

    col = pow(col, vec3f(2.2));
    textureStore(screen, id.xy, vec4f(col, 1.));
}
fn pcg(v: u32) -> u32 {
    let state = v * 747796405u + 2891336453u;
    let word = ((state >> ((state >> 28u) + 4u)) ^ state) * 277803737u;
    return (word >> 22u) ^ word;
}

fn hash(c: vec2i) -> f32 {
    return f32(pcg(bitcast<u32>(c.x) ^ pcg(bitcast<u32>(c.y))) >> 8u) / 16777216.0;
}

fn noise(x: vec2f) -> f32 {
    let i = vec2i(floor(x));
    let f = fract(x);
    let u = f * f * (3.0 - 2.0 * f);          // smooth blend weights

    let a = hash(i);
    let b = hash(i + vec2i(1, 0));
    let c = hash(i + vec2i(0, 1));
    let d = hash(i + vec2i(1, 1));
    return mix(mix(a, b, u.x), mix(c, d, u.x), u.y);
}

// Five octaves: each doubles the frequency and halves the amplitude
fn fbm(x0: vec2f) -> f32 {
    var x = x0;
    var amp = 0.5;
    var sum = 0.0;
    for (var i = 0; i < 5; i++) {
        sum += amp * noise(x);
        x = x * 2.0 + vec2f(3.1, 1.7);         // offset so octaves don't align
        amp *= 0.5;
    }
    return sum;
}

@compute @workgroup_size(16, 16)
fn main_image(@builtin(global_invocation_id) id: vec3u) {
    let screen_size = textureDimensions(screen);
    if (id.x >= screen_size.x || id.y >= screen_size.y) { return; }

    let fragCoord = vec2f(f32(id.x) + .5, f32(screen_size.y - id.y) - .5);
    let p = (2.0 * fragCoord - vec2f(screen_size)) / f32(screen_size.y);

    let n = fbm(p * 2.5 + vec2f(time.elapsed * 0.3, 0.0));
    var col = mix(vec3f(0.1, 0.25, 0.6), vec3f(1.0), smoothstep(0.3, 0.8, n));

    col = pow(col, vec3f(2.2));
    textureStore(screen, id.xy, vec4f(col, 1.));
}

Seeded generators

A hash gives one number per input. When one pixel needs many random numbers, it’s more convenient to wrap the hash in a tiny generator: keep a u32 state, and advance it each time you draw a number.

// A starting state from the pixel, the frame, and your own seed
fn rng_init(pixel: vec2u, frame: u32, seed: u32) -> u32 {
    return pcg(pixel.x ^ pcg(pixel.y ^ pcg(frame ^ pcg(seed))));
}

// The next float in [0, 1); updates the state through the pointer
fn rng_next(state: ptr<function, u32>) -> f32 {
    *state = pcg(*state);
    return f32(*state >> 8u) / 16777216.0;
}

var rng = rng_init(id.xy, time.frame, 1234u);
let a = rng_next(&rng);
let b = rng_next(&rng);    // a different number

The state must be a var, since it changes, and &rng passes a pointer so the function can update it.

Give each part of your shader its own stream with a different seed, for example rng_init(id.xy, time.frame, 1u) for color and 2u for shape. With one shared generator, inserting one extra rng_next call shifts every number after it, which makes tweaking confusing. This is the same idea behind JAX’s split and fold_in: random streams are derived from keys, not from a shared, mutating global.

Part V More than one pass

So far, every shader has had one entry point that computes each pixel from scratch. Fragment shaders on Shadertoy work the same way. Compute shaders can do more: run several entry points per frame, pass images from one to the next, and keep data from one frame to the following one.

13 Blur pass

A blur needs each pixel to read its neighbors. Within one pass, you can’t read what other threads are computing: there’s no order between them, so a neighbor may not have written yet. The fix is to split the work. One entry point draws the image; a second one reads the finished image and blurs it. The boundary between entry points guarantees that everything from the first is done before the second starts.

In this playground (as in compute.toys), entry points hand images to each other through pass textures: four scratch images, each the size of the screen, with four channels per pixel.

passStore(layer, coord, value)   // write a vec4f to pixel coord of scratch image `layer` (0–3)
passLoad(layer, coord, lod)      // read it back; `lod` is the mip level, always 0 here

A few details:

A box blur replaces each pixel with the average of the square of pixels around it. With radius R, that’s (2R + 1)² samples. Pixels near the edge would read outside the image, so clamp the coordinate to the valid range.

The scene pass draws a sharp picture into layer 0. In main_image, replace each pixel with the average of the 13 × 13 square around it. Show the blurred version only to the right of the mouse (or of the middle, before the mouse moves), so you can compare.

fn sdCircle(p: vec2f, r: f32) -> f32 {
    return length(p) - r;
}

// Pass 1: draw a sharp scene into scratch layer 0
@compute @workgroup_size(16, 16)
fn scene(@builtin(global_invocation_id) id: vec3u) {
    let screen_size = textureDimensions(screen);
    if (id.x >= screen_size.x || id.y >= screen_size.y) { return; }

    let fragCoord = vec2f(f32(id.x) + .5, f32(screen_size.y - id.y) - .5);
    let p = (2.0 * fragCoord - vec2f(screen_size)) / f32(screen_size.y);

    // A checkerboard with orbiting circles: lots of hard edges
    let cell = floor(p * 6.0);
    let checker = abs((cell.x + cell.y) % 2.0);
    var col = mix(vec3f(0.06, 0.07, 0.12), vec3f(0.2, 0.22, 0.32), checker);
    for (var i = 0; i < 3; i++) {
        let a = time.elapsed * 0.7 + f32(i) * 2.094;
        let d = sdCircle(p - 0.55 * vec2f(cos(a), sin(a)), 0.22);
        let tint = 0.55 + 0.45 * cos(vec3f(0.0, 2.0, 4.0) + f32(i) * 2.0);
        col = select(col, tint, d < 0.0);
    }

    // Store in linear light so the blur averages real brightness
    passStore(0, vec2i(id.xy), vec4f(pow(col, vec3f(2.2)), 1.0));
}

// Pass 2: read layer 0 and write the screen
@compute @workgroup_size(16, 16)
fn main_image(@builtin(global_invocation_id) id: vec3u) {
    let screen_size = textureDimensions(screen);
    if (id.x >= screen_size.x || id.y >= screen_size.y) { return; }

    let here = vec2i(id.xy);

    // TODO: average the (2R+1) x (2R+1) square around `here`,
    //       clamping coordinates to the screen
    let col = passLoad(0, here, 0).rgb;

    textureStore(screen, id.xy, vec4f(col, 1.));
}
fn sdCircle(p: vec2f, r: f32) -> f32 {
    return length(p) - r;
}

@compute @workgroup_size(16, 16)
fn scene(@builtin(global_invocation_id) id: vec3u) {
    let screen_size = textureDimensions(screen);
    if (id.x >= screen_size.x || id.y >= screen_size.y) { return; }

    let fragCoord = vec2f(f32(id.x) + .5, f32(screen_size.y - id.y) - .5);
    let p = (2.0 * fragCoord - vec2f(screen_size)) / f32(screen_size.y);

    let cell = floor(p * 6.0);
    let checker = abs((cell.x + cell.y) % 2.0);
    var col = mix(vec3f(0.06, 0.07, 0.12), vec3f(0.2, 0.22, 0.32), checker);
    for (var i = 0; i < 3; i++) {
        let a = time.elapsed * 0.7 + f32(i) * 2.094;
        let d = sdCircle(p - 0.55 * vec2f(cos(a), sin(a)), 0.22);
        let tint = 0.55 + 0.45 * cos(vec3f(0.0, 2.0, 4.0) + f32(i) * 2.0);
        col = select(col, tint, d < 0.0);
    }
    passStore(0, vec2i(id.xy), vec4f(pow(col, vec3f(2.2)), 1.0));
}

const R = 6;

@compute @workgroup_size(16, 16)
fn main_image(@builtin(global_invocation_id) id: vec3u) {
    let screen_size = textureDimensions(screen);
    if (id.x >= screen_size.x || id.y >= screen_size.y) { return; }

    let here = vec2i(id.xy);
    let last = vec2i(screen_size) - 1;

    var sum = vec3f(0.0);
    for (var dy = -R; dy <= R; dy++) {
        for (var dx = -R; dx <= R; dx++) {
            let q = clamp(here + vec2i(dx, dy), vec2i(0), last);
            sum += passLoad(0, q, 0).rgb;
        }
    }
    let blurred = sum / f32((2 * R + 1) * (2 * R + 1));
    let sharp = passLoad(0, here, 0).rgb;

    // Compare: sharp on the left of the split, blurred on the right
    var split = i32(screen_size.x) / 2;
    if (mouse.pos.x > 0) { split = mouse.pos.x; }
    var col = select(sharp, blurred, here.x > split);
    if (abs(here.x - split) < 1) { col = vec3f(1.0); }

    textureStore(screen, id.xy, vec4f(col, 1.));
}

A box blur of radius R costs (2R + 1)² reads per pixel: 169 here, and 1,681 at radius 20. The standard improvement is a separable blur. Blur horizontally into layer 1, then blur that vertically into the screen. The result is identical, and it costs 2 × (2R + 1) reads instead. Try it: add a third entry point between scene and main_image.

14 Game of Life

Pass textures keep their contents between frames. That makes them memory: a simulation can read last frame’s state and write this frame’s state.

Conway’s Game of Life is the classic first stateful compute shader. The grid is made of cells that are alive or dead. Every step, each cell counts its eight neighbors:

You can’t update the grid in place, because a cell’s neighbors would see a mix of old and new values depending on which threads happened to run first. Every cell must read the old generation and write the new one. The usual pattern keeps two buffers and swaps their roles each frame, known as ping-pong. Here the runtime does it for you, since passLoad always reads the previous state and passStore writes the next one.

Three practical details:

Implement the rules in life. Right now it copies each cell unchanged, so the random soup just sits there. Hold the mouse button on the canvas to draw new live cells.

const CELL = 4u;   // each cell is drawn as a CELL x CELL block of pixels

fn pcg(v: u32) -> u32 {
    let state = v * 747796405u + 2891336453u;
    let word = ((state >> ((state >> 28u) + 4u)) ^ state) * 277803737u;
    return (word >> 22u) ^ word;
}

// State of a cell (1.0 alive, 0.0 dead), wrapping around the edges
fn alive(c: vec2i, grid: vec2i) -> f32 {
    let w = (c + grid) % grid;
    return passLoad(0, w, 0).x;
}

// Step: one thread per cell. x = alive, y = a fading trail for display
@compute @workgroup_size(16, 16)
fn life(@builtin(global_invocation_id) id: vec3u) {
    let grid = vec2i(textureDimensions(screen) / CELL);
    let c = vec2i(id.xy);
    if (c.x >= grid.x || c.y >= grid.y) { return; }

    let prev = passLoad(0, c, 0);

    // First frame: random soup, about 30% alive
    if (time.frame == 0u) {
        let r = f32(pcg(id.x ^ pcg(id.y ^ 12345u)) >> 8u) / 16777216.0;
        let a = select(0.0, 1.0, r < 0.3);
        passStore(0, c, vec4f(a, a, 0.0, 1.0));
        return;
    }

    let me = alive(c, grid);

    // TODO: count the 8 neighbors with alive(), then apply the rules
    var next = me;

    // Hold the mouse to draw live cells
    let m = mouse.pos / i32(CELL);
    if (mouse.click == 1 && length(vec2f(c - m)) < 3.0) { next = 1.0; }

    passStore(0, c, vec4f(next, max(next, prev.y * 0.92), 0.0, 1.0));
}

// Display: read the cell under this pixel
@compute @workgroup_size(16, 16)
fn main_image(@builtin(global_invocation_id) id: vec3u) {
    let screen_size = textureDimensions(screen);
    if (id.x >= screen_size.x || id.y >= screen_size.y) { return; }

    let s = passLoad(0, vec2i(id.xy / CELL), 0);
    var col = mix(vec3f(0.03, 0.035, 0.07), vec3f(0.2, 0.28, 0.75), s.y);
    col = mix(col, vec3f(1.0, 0.9, 0.7), s.x);

    col = pow(col, vec3f(2.2));
    textureStore(screen, id.xy, vec4f(col, 1.));
}
const CELL = 4u;

fn pcg(v: u32) -> u32 {
    let state = v * 747796405u + 2891336453u;
    let word = ((state >> ((state >> 28u) + 4u)) ^ state) * 277803737u;
    return (word >> 22u) ^ word;
}

fn alive(c: vec2i, grid: vec2i) -> f32 {
    let w = (c + grid) % grid;
    return passLoad(0, w, 0).x;
}

@compute @workgroup_size(16, 16)
fn life(@builtin(global_invocation_id) id: vec3u) {
    let grid = vec2i(textureDimensions(screen) / CELL);
    let c = vec2i(id.xy);
    if (c.x >= grid.x || c.y >= grid.y) { return; }

    let prev = passLoad(0, c, 0);

    if (time.frame == 0u) {
        let r = f32(pcg(id.x ^ pcg(id.y ^ 12345u)) >> 8u) / 16777216.0;
        let a = select(0.0, 1.0, r < 0.3);
        passStore(0, c, vec4f(a, a, 0.0, 1.0));
        return;
    }

    let me = alive(c, grid);

    var n = 0.0;
    for (var dy = -1; dy <= 1; dy++) {
        for (var dx = -1; dx <= 1; dx++) {
            if (dx == 0 && dy == 0) { continue; }
            n += alive(c + vec2i(dx, dy), grid);
        }
    }
    // Birth on exactly 3; survival on 2 or 3
    var next = select(0.0, 1.0, n == 3.0 || (me == 1.0 && n == 2.0));

    let m = mouse.pos / i32(CELL);
    if (mouse.click == 1 && length(vec2f(c - m)) < 3.0) { next = 1.0; }

    passStore(0, c, vec4f(next, max(next, prev.y * 0.92), 0.0, 1.0));
}

@compute @workgroup_size(16, 16)
fn main_image(@builtin(global_invocation_id) id: vec3u) {
    let screen_size = textureDimensions(screen);
    if (id.x >= screen_size.x || id.y >= screen_size.y) { return; }

    let s = passLoad(0, vec2i(id.xy / CELL), 0);
    var col = mix(vec3f(0.03, 0.035, 0.07), vec3f(0.2, 0.28, 0.75), s.y);
    col = mix(col, vec3f(1.0, 0.9, 0.7), s.x);

    col = pow(col, vec3f(2.2));
    textureStore(screen, id.xy, vec4f(col, 1.));
}

Part VI Threads that cooperate

Everything so far has been embarrassingly parallel: each thread’s output depends only on its inputs, never on other threads in the same pass. The last three exercises aren’t like that. They need to combine thousands of values into one, which is the heart of GPU computing for machine learning.

15 Parallel sum

The task: add up 16,384 numbers. One thread doing it alone would take 16,384 steps while thousands of other threads sit idle. The parallel approach uses three new tools.

Storage buffers

A storage buffer is a plain array in GPU memory that any thread can read and write. Declare one with a directive at the top of the file:

#storage data array<f32>
#storage partials array<f32>

The runtime turns each line into a buffer binding, so data[i] works like any array. Buffers here hold 2 million floats each. compute.toys allows two storage buffers, so the exercises stick to two, and pack several things into one buffer when they need more.

Dispatch sizes

By default, an entry point gets enough workgroups to cover the screen. That’s wrong for a 1D array, so set it explicitly:

#workgroup_count fill 64 1 1     // 64 workgroups × 256 threads = 16,384 threads

The directives only configure entry points; you still write fill, reduce and finish as ordinary @compute functions with exactly those names.

Workgroup memory and barriers

Threads in the same workgroup can share a small, fast scratchpad:

var<workgroup> scratch: array<f32, 256>;

Every workgroup gets its own copy, visible to all 256 of its threads. To use it safely, threads must wait for each other. workgroupBarrier() makes each thread pause until every thread in its workgroup has reached the same point, and makes their writes to workgroup memory visible.

The tree reduction

Inside one workgroup, the 256 values are summed in halving steps:

  1. Each thread loads one value into scratch[lid.x]. Barrier.
  2. With stride = 128, each thread with lid.x < stride adds scratch[lid.x + stride] into its own slot. Barrier.
  3. Repeat with stride 64, 32, … down to 1. Barrier after every step.
  4. scratch[0] now holds the workgroup’s total. Thread 0 writes it out.
load31415926stride 481067stride 21417stride 131workgroupBarrier()
Eight values summed in three steps instead of seven. With 256 values, it’s 8 steps. Each dashed line is a barrier.

There’s no barrier across workgroups: workgroups may run in any order, and not necessarily at the same time. So each workgroup writes its total to partials[workgroup_id], and a second entry point, finish, runs a single workgroup that reduces those 64 partials the same way. The boundary between entry points is the global synchronization you don’t otherwise have.

Every thread must reach every barrier. Don’t put workgroupBarrier() inside if (lid.x < stride), and don’t return early before one. WGSL’s uniformity analysis rejects code where a barrier might be skipped by some threads. Keep the barrier inside the loop but outside the if.

lid versus gid. “Once per workgroup” is always lid.x == 0u. gid.x == 0u is true for exactly one thread in the entire dispatch, so only partials[0] would ever be written. This mix-up is probably the most common compute shader bug there is.

When the data doesn’t divide evenly into workgroups, don’t skip threads, since that would break the barrier rule. Load the identity element instead: 0 for a sum, -inf for a max, 1 for a product. Here 64 × 256 is exactly 16,384, so it doesn’t come up.

Fill in reduce and finish. The display shows the total as a bar, green when it equals 16,384 exactly and red otherwise. The strip underneath shows the 64 partial sums: each should be equally bright. If only the first one lights up, check lid versus gid. If the total flickers between frames, a barrier is missing.

#storage data array<f32>
#storage partials array<f32>
#workgroup_count fill 64 1 1
#workgroup_count reduce 64 1 1
#workgroup_count finish 1 1 1

const N = 16384u;      // 64 workgroups x 256 threads
const RESULT = 64u;    // partials[0..64) = per-workgroup sums, partials[64] = total

var<workgroup> scratch: array<f32, 256>;

@compute @workgroup_size(256)
fn fill(@builtin(global_invocation_id) gid: vec3u) {
    data[gid.x] = 1.0;   // all ones, so the correct total is 16384
}

@compute @workgroup_size(256)
fn reduce(@builtin(global_invocation_id) gid: vec3u,
          @builtin(local_invocation_id) lid: vec3u,
          @builtin(workgroup_id) wid: vec3u) {
    // TODO: load data[gid.x] into scratch[lid.x], then barrier
    // TODO: tree reduction, stride 128, 64, ..., 1, barrier after each step
    // TODO: thread 0 of each workgroup writes scratch[0] to partials[wid.x]
}

@compute @workgroup_size(64)
fn finish(@builtin(local_invocation_id) lid: vec3u) {
    // TODO: the same reduction over the 64 partials, with one workgroup.
    //       Reuse scratch (first 64 slots). Thread 0 writes partials[RESULT].
}

// ---------- display: total as a bar, partials as a strip ----------
@compute @workgroup_size(16, 16)
fn main_image(@builtin(global_invocation_id) id: vec3u) {
    let size = textureDimensions(screen);
    if (id.x >= size.x || id.y >= size.y) { return; }
    let uv = (vec2f(id.xy) + 0.5) / vec2f(size);   // y points down here

    let total = partials[RESULT];
    var col = vec3f(0.06, 0.07, 0.12);
    if (uv.y > 0.3 && uv.y < 0.5) {
        col = vec3f(0.14, 0.16, 0.26);
        if (uv.x < total / f32(N)) {
            col = select(vec3f(0.95, 0.35, 0.3), vec3f(0.3, 0.85, 0.5), abs(total - f32(N)) < 0.5);
        }
    }
    let quarter = abs(fract(uv.x * 4.0 + 0.5) - 0.5) * f32(size.x) / 4.0;
    if (quarter < 1.0 && uv.y > 0.27 && uv.y < 0.53) { col = vec3f(0.55); }
    if (uv.y > 0.62 && uv.y < 0.72) {
        let i = min(u32(uv.x * 64.0), 63u);
        col = vec3f(0.35, 0.6, 1.0) * clamp(partials[i] / 256.0, 0.0, 1.0);
    }
    textureStore(screen, id.xy, vec4f(pow(col, vec3f(2.2)), 1.));
}
#storage data array<f32>
#storage partials array<f32>
#workgroup_count fill 64 1 1
#workgroup_count reduce 64 1 1
#workgroup_count finish 1 1 1

const N = 16384u;
const RESULT = 64u;

var<workgroup> scratch: array<f32, 256>;

@compute @workgroup_size(256)
fn fill(@builtin(global_invocation_id) gid: vec3u) {
    data[gid.x] = 1.0;
}

@compute @workgroup_size(256)
fn reduce(@builtin(global_invocation_id) gid: vec3u,
          @builtin(local_invocation_id) lid: vec3u,
          @builtin(workgroup_id) wid: vec3u) {
    scratch[lid.x] = data[gid.x];   // gid picks the element, lid picks the slot
    workgroupBarrier();

    for (var stride = 128u; stride > 0u; stride >>= 1u) {
        if (lid.x < stride) {
            scratch[lid.x] += scratch[lid.x + stride];
        }
        workgroupBarrier();          // outside the if: every thread reaches it
    }

    if (lid.x == 0u) {               // once per workgroup: lid, not gid
        partials[wid.x] = scratch[0];
    }
}

@compute @workgroup_size(64)
fn finish(@builtin(local_invocation_id) lid: vec3u) {
    scratch[lid.x] = partials[lid.x];
    workgroupBarrier();

    for (var stride = 32u; stride > 0u; stride >>= 1u) {
        if (lid.x < stride) {
            scratch[lid.x] += scratch[lid.x + stride];
        }
        workgroupBarrier();
    }

    if (lid.x == 0u) {
        partials[RESULT] = scratch[0];
    }
}

@compute @workgroup_size(16, 16)
fn main_image(@builtin(global_invocation_id) id: vec3u) {
    let size = textureDimensions(screen);
    if (id.x >= size.x || id.y >= size.y) { return; }
    let uv = (vec2f(id.xy) + 0.5) / vec2f(size);

    let total = partials[RESULT];
    var col = vec3f(0.06, 0.07, 0.12);
    if (uv.y > 0.3 && uv.y < 0.5) {
        col = vec3f(0.14, 0.16, 0.26);
        if (uv.x < total / f32(N)) {
            col = select(vec3f(0.95, 0.35, 0.3), vec3f(0.3, 0.85, 0.5), abs(total - f32(N)) < 0.5);
        }
    }
    let quarter = abs(fract(uv.x * 4.0 + 0.5) - 0.5) * f32(size.x) / 4.0;
    if (quarter < 1.0 && uv.y > 0.27 && uv.y < 0.53) { col = vec3f(0.55); }
    if (uv.y > 0.62 && uv.y < 0.72) {
        let i = min(u32(uv.x * 64.0), 63u);
        col = vec3f(0.35, 0.6, 1.0) * clamp(partials[i] / 256.0, 0.0, 1.0);
    }
    textureStore(screen, id.xy, vec4f(pow(col, vec3f(2.2)), 1.));
}

Once it’s green, change fill to random values and verify the sum a second, dumb way: one entry point with #workgroup_count check 1 1 1 in which a single thread loops over all 16,384 values. The two results won’t match to the last bit, because floating-point addition isn’t associative and the two methods add in different orders. That difference is normal, and worth seeing once.

This reduction is the single most important pattern on this page. Softmax, layer norm and attention all contain it.

16 Softmax

Softmax turns a list of scores into probabilities that are positive and sum to 1:

softmax(x)i=exi−m∑jexj−m\operatorname{softmax}(x)_i = \frac{e^{x_i - m}}{\sum_j e^{x_j - m}}

where mm is the largest of the xjx_j. Subtracting the maximum m doesn’t change the answer, since it scales the top and bottom by the same factor. But it keeps exp from overflowing: exp(89.0) is already too large for an f32. Every real implementation does this.

On the GPU, softmax is two reductions and an elementwise pass:

  1. Max reduction: the same tree as the sum, with max instead of +.
  2. Sum reduction of exp(x - m).
  3. Normalize: every element becomes exp(x - m) / sum, written back in place.

Each reduction is two entry points (per-workgroup, then finish), so you’ll write five new entry points. They run in the order they appear in the file.

Two things to plan yourself:

Implement steps 1–3 between fill and check_reduce. The top bar is the sum of your outputs and turns green when it’s 1. The bottom bar is the probability of element 5,000, which fill sets to 20.0, the largest value in the array. It should end up close to 0.96.

#storage data array<f32>
#storage stats array<f32>
#workgroup_count fill 64 1 1
#workgroup_count check_reduce 64 1 1
#workgroup_count check_finish 1 1 1

const N = 16384u;
const PLANT = 5000u;           // fill puts the largest value here

// Layout of the `stats` buffer
const MAX_PARTIALS = 0u;       // [0, 64)
const MAX = 64u;
const SUM_PARTIALS = 128u;     // [128, 192)
const SUM = 192u;
const CHECK_PARTIALS = 256u;   // [256, 320)
const CHECK = 320u;

var<workgroup> scratch: array<f32, 256>;

fn pcg(v: u32) -> u32 {
    let state = v * 747796405u + 2891336453u;
    let word = ((state >> ((state >> 28u) + 4u)) ^ state) * 277803737u;
    return (word >> 22u) ^ word;
}

// Tree-sum scratch[0..width) into scratch[0]. Call from every thread.
fn tree_sum(lid: u32, width: u32) {
    for (var s = width / 2u; s > 0u; s >>= 1u) {
        if (lid < s) { scratch[lid] += scratch[lid + s]; }
        workgroupBarrier();
    }
}

// Random scores in [-10, 10), with one planted maximum of 20
@compute @workgroup_size(256)
fn fill(@builtin(global_invocation_id) gid: vec3u) {
    let r = f32(pcg(gid.x) >> 8u) / 16777216.0;
    data[gid.x] = select(r * 20.0 - 10.0, 20.0, gid.x == PLANT);
}

// TODO: your entry points go here, in order, each with a #workgroup_count:
//   reduce_max, finish_max  ->  stats[MAX]
//   reduce_sum, finish_sum  ->  stats[SUM]   (sum of exp(x - max))
//   normalize               ->  data[i] = exp(data[i] - max) / sum

// ---------- check: sum the outputs, which should be 1 ----------
@compute @workgroup_size(256)
fn check_reduce(@builtin(global_invocation_id) gid: vec3u,
                @builtin(local_invocation_id) lid: vec3u,
                @builtin(workgroup_id) wid: vec3u) {
    scratch[lid.x] = data[gid.x];
    workgroupBarrier();
    tree_sum(lid.x, 256u);
    if (lid.x == 0u) { stats[CHECK_PARTIALS + wid.x] = scratch[0]; }
}

@compute @workgroup_size(64)
fn check_finish(@builtin(local_invocation_id) lid: vec3u) {
    scratch[lid.x] = stats[CHECK_PARTIALS + lid.x];
    workgroupBarrier();
    tree_sum(lid.x, 64u);
    if (lid.x == 0u) { stats[CHECK] = scratch[0]; }
}

@compute @workgroup_size(16, 16)
fn main_image(@builtin(global_invocation_id) id: vec3u) {
    let size = textureDimensions(screen);
    if (id.x >= size.x || id.y >= size.y) { return; }
    let uv = (vec2f(id.xy) + 0.5) / vec2f(size);

    let total = stats[CHECK];
    let planted = data[PLANT];
    var col = vec3f(0.06, 0.07, 0.12);
    if (uv.y > 0.25 && uv.y < 0.42) {
        col = vec3f(0.14, 0.16, 0.26);
        if (uv.x < total) {
            col = select(vec3f(0.95, 0.35, 0.3), vec3f(0.3, 0.85, 0.5), abs(total - 1.0) < 1e-3);
        }
    }
    if (uv.y > 0.58 && uv.y < 0.75) {
        col = vec3f(0.14, 0.16, 0.26);
        if (uv.x < planted) { col = vec3f(0.35, 0.6, 1.0); }
    }
    textureStore(screen, id.xy, vec4f(pow(col, vec3f(2.2)), 1.));
}
#storage data array<f32>
#storage stats array<f32>
#workgroup_count fill 64 1 1
#workgroup_count reduce_max 64 1 1
#workgroup_count finish_max 1 1 1
#workgroup_count reduce_sum 64 1 1
#workgroup_count finish_sum 1 1 1
#workgroup_count normalize 64 1 1
#workgroup_count check_reduce 64 1 1
#workgroup_count check_finish 1 1 1

const N = 16384u;
const PLANT = 5000u;

const MAX_PARTIALS = 0u;
const MAX = 64u;
const SUM_PARTIALS = 128u;
const SUM = 192u;
const CHECK_PARTIALS = 256u;
const CHECK = 320u;

var<workgroup> scratch: array<f32, 256>;

fn pcg(v: u32) -> u32 {
    let state = v * 747796405u + 2891336453u;
    let word = ((state >> ((state >> 28u) + 4u)) ^ state) * 277803737u;
    return (word >> 22u) ^ word;
}

fn tree_sum(lid: u32, width: u32) {
    for (var s = width / 2u; s > 0u; s >>= 1u) {
        if (lid < s) { scratch[lid] += scratch[lid + s]; }
        workgroupBarrier();
    }
}

// Same tree, different combining operation
fn tree_max(lid: u32, width: u32) {
    for (var s = width / 2u; s > 0u; s >>= 1u) {
        if (lid < s) { scratch[lid] = max(scratch[lid], scratch[lid + s]); }
        workgroupBarrier();
    }
}

@compute @workgroup_size(256)
fn fill(@builtin(global_invocation_id) gid: vec3u) {
    let r = f32(pcg(gid.x) >> 8u) / 16777216.0;
    data[gid.x] = select(r * 20.0 - 10.0, 20.0, gid.x == PLANT);
}

// ---- 1. max ----
@compute @workgroup_size(256)
fn reduce_max(@builtin(global_invocation_id) gid: vec3u,
              @builtin(local_invocation_id) lid: vec3u,
              @builtin(workgroup_id) wid: vec3u) {
    scratch[lid.x] = data[gid.x];
    workgroupBarrier();
    tree_max(lid.x, 256u);
    if (lid.x == 0u) { stats[MAX_PARTIALS + wid.x] = scratch[0]; }
}

@compute @workgroup_size(64)
fn finish_max(@builtin(local_invocation_id) lid: vec3u) {
    scratch[lid.x] = stats[MAX_PARTIALS + lid.x];
    workgroupBarrier();
    tree_max(lid.x, 64u);
    if (lid.x == 0u) { stats[MAX] = scratch[0]; }
}

// ---- 2. sum of exp(x - max) ----
@compute @workgroup_size(256)
fn reduce_sum(@builtin(global_invocation_id) gid: vec3u,
              @builtin(local_invocation_id) lid: vec3u,
              @builtin(workgroup_id) wid: vec3u) {
    scratch[lid.x] = exp(data[gid.x] - stats[MAX]);
    workgroupBarrier();
    tree_sum(lid.x, 256u);
    if (lid.x == 0u) { stats[SUM_PARTIALS + wid.x] = scratch[0]; }
}

@compute @workgroup_size(64)
fn finish_sum(@builtin(local_invocation_id) lid: vec3u) {
    scratch[lid.x] = stats[SUM_PARTIALS + lid.x];
    workgroupBarrier();
    tree_sum(lid.x, 64u);
    if (lid.x == 0u) { stats[SUM] = scratch[0]; }
}

// ---- 3. normalize in place ----
@compute @workgroup_size(256)
fn normalize(@builtin(global_invocation_id) gid: vec3u) {
    data[gid.x] = exp(data[gid.x] - stats[MAX]) / stats[SUM];
}

// ---------- check ----------
@compute @workgroup_size(256)
fn check_reduce(@builtin(global_invocation_id) gid: vec3u,
                @builtin(local_invocation_id) lid: vec3u,
                @builtin(workgroup_id) wid: vec3u) {
    scratch[lid.x] = data[gid.x];
    workgroupBarrier();
    tree_sum(lid.x, 256u);
    if (lid.x == 0u) { stats[CHECK_PARTIALS + wid.x] = scratch[0]; }
}

@compute @workgroup_size(64)
fn check_finish(@builtin(local_invocation_id) lid: vec3u) {
    scratch[lid.x] = stats[CHECK_PARTIALS + lid.x];
    workgroupBarrier();
    tree_sum(lid.x, 64u);
    if (lid.x == 0u) { stats[CHECK] = scratch[0]; }
}

@compute @workgroup_size(16, 16)
fn main_image(@builtin(global_invocation_id) id: vec3u) {
    let size = textureDimensions(screen);
    if (id.x >= size.x || id.y >= size.y) { return; }
    let uv = (vec2f(id.xy) + 0.5) / vec2f(size);

    let total = stats[CHECK];
    let planted = data[PLANT];
    var col = vec3f(0.06, 0.07, 0.12);
    if (uv.y > 0.25 && uv.y < 0.42) {
        col = vec3f(0.14, 0.16, 0.26);
        if (uv.x < total) {
            col = select(vec3f(0.95, 0.35, 0.3), vec3f(0.3, 0.85, 0.5), abs(total - 1.0) < 1e-3);
        }
    }
    if (uv.y > 0.58 && uv.y < 0.75) {
        col = vec3f(0.14, 0.16, 0.26);
        if (uv.x < planted) { col = vec3f(0.35, 0.6, 1.0); }
    }
    textureStore(screen, id.xy, vec4f(pow(col, vec3f(2.2)), 1.));
}

Why is the planted element’s probability about 0.96? Its term is e^0 = 1 after subtracting the max. The other 16,383 values are spread evenly over [−10, 10), so each contributes on average about e−20×1101≈2.3×10−6e^{-20} \times 1101 \approx 2.3 \times 10^{-6}, for a combined total near 0.037. That gives 1 / 1.037 ≈ 0.964. Checking the GPU’s answer against a number you can work out by hand is the habit that matters most for ML kernels, where a slightly wrong result produces no visible symptom.

Part VII Matrix multiply

17 Matrix multiply

Matrix multiplication is the operation that most of machine learning runs on. For square N × N matrices, C = A × B means

Cij=∑k=0N−1AikBkj.C_{ij} = \sum_{k=0}^{N-1} A_{ik} \, B_{kj}.

The straightforward GPU version gives one thread per output element: thread (j, i) loops over k and computes C[i][j]. That’s a 2D dispatch: @workgroup_size(16, 16) with #workgroup_count matmul 16 16 1 covers a 256 × 256 output.

2D data in flat buffers

Storage buffers are one-dimensional, so a matrix is stored row by row. Element (row, col) of an N-wide matrix lives at index row * N + col. With only two buffers, pack A and B into one buffer at different offsets: A starts at 0 and B starts at N * N.

Conversions, again

WGSL’s conversion functions match the shape of what they produce. f32(x) converts a scalar; for vectors you use the vector constructor, which converts every component:

let ij = vec2u(uv * vec2f(DIMS));   // not u32(...): that only takes a scalar

Converting a float to u32 truncates toward zero, so for non-negative values it’s the same as floor.

Knowing it’s right

Fill A with the identity matrix (1 on the diagonal, 0 elsewhere) and B with random values. Then C must equal B exactly, since each output is one 1 × b plus zeros. The display shows C on the left, and on the right a map of where C matches B: green where it does, red where it doesn’t.

Write matmul. With the starter’s placeholder, C is all zeros, so the left panel is black and the right one is red. When it’s correct, the left panel looks like static and the right one is solid green.

#storage inputs array<f32>
#storage outputs array<f32>
#workgroup_count fill 16 16 1
#workgroup_count matmul 16 16 1

const N = 256u;          // matrices are N x N; 16 x 16 workgroups of 16 x 16
const A = 0u;            // A starts here in `inputs`
const B = N * N;         // B starts here in `inputs`

fn pcg(v: u32) -> u32 {
    let state = v * 747796405u + 2891336453u;
    let word = ((state >> ((state >> 28u) + 4u)) ^ state) * 277803737u;
    return (word >> 22u) ^ word;
}

// A = identity, B = random. gid.y is the row, gid.x the column.
@compute @workgroup_size(16, 16)
fn fill(@builtin(global_invocation_id) gid: vec3u) {
    let index = gid.y * N + gid.x;
    inputs[A + index] = select(0.0, 1.0, gid.x == gid.y);
    inputs[B + index] = f32(pcg(index) >> 8u) / 16777216.0;
}

@compute @workgroup_size(16, 16)
fn matmul(@builtin(global_invocation_id) gid: vec3u) {
    let row = gid.y;
    let col = gid.x;
    // TODO: outputs[row * N + col] = sum over k of A[row][k] * B[k][col]
    outputs[row * N + col] = 0.0;
}

// ---------- display: C on the left, C == B check on the right ----------
@compute @workgroup_size(16, 16)
fn main_image(@builtin(global_invocation_id) id: vec3u) {
    let size = textureDimensions(screen);
    if (id.x >= size.x || id.y >= size.y) { return; }

    let side = min(size.x / 2u, size.y);
    let panel = id.x / side;
    let local = vec2u(id.x % side, id.y);
    if (panel > 1u || local.y >= side) {
        textureStore(screen, id.xy, vec4f(0.004, 0.004, 0.008, 1.));
        return;
    }
    let ij = min(local * N / side, vec2u(N - 1u));    // (column, row)
    let index = ij.y * N + ij.x;
    let c = outputs[index];
    let ok = abs(c - inputs[B + index]) < 1e-5;

    var col = vec3f(c);
    if (panel == 1u) { col = select(vec3f(0.9, 0.2, 0.15), vec3f(0.2, 0.75, 0.4), ok); }
    textureStore(screen, id.xy, vec4f(pow(col, vec3f(2.2)), 1.));
}
#storage inputs array<f32>
#storage outputs array<f32>
#workgroup_count fill 16 16 1
#workgroup_count matmul 16 16 1

const N = 256u;
const A = 0u;
const B = N * N;

fn pcg(v: u32) -> u32 {
    let state = v * 747796405u + 2891336453u;
    let word = ((state >> ((state >> 28u) + 4u)) ^ state) * 277803737u;
    return (word >> 22u) ^ word;
}

@compute @workgroup_size(16, 16)
fn fill(@builtin(global_invocation_id) gid: vec3u) {
    let index = gid.y * N + gid.x;
    inputs[A + index] = select(0.0, 1.0, gid.x == gid.y);
    inputs[B + index] = f32(pcg(index) >> 8u) / 16777216.0;
}

@compute @workgroup_size(16, 16)
fn matmul(@builtin(global_invocation_id) gid: vec3u) {
    let row = gid.y;
    let col = gid.x;
    var acc = 0.0;                               // accumulate in a register
    for (var k = 0u; k < N; k++) {
        acc += inputs[A + row * N + k] * inputs[B + k * N + col];
    }
    outputs[row * N + col] = acc;                // one write at the end
}

@compute @workgroup_size(16, 16)
fn main_image(@builtin(global_invocation_id) id: vec3u) {
    let size = textureDimensions(screen);
    if (id.x >= size.x || id.y >= size.y) { return; }

    let side = min(size.x / 2u, size.y);
    let panel = id.x / side;
    let local = vec2u(id.x % side, id.y);
    if (panel > 1u || local.y >= side) {
        textureStore(screen, id.xy, vec4f(0.004, 0.004, 0.008, 1.));
        return;
    }
    let ij = min(local * N / side, vec2u(N - 1u));
    let index = ij.y * N + ij.x;
    let c = outputs[index];
    let ok = abs(c - inputs[B + index]) < 1e-5;

    var col = vec3f(c);
    if (panel == 1u) { col = select(vec3f(0.9, 0.2, 0.15), vec3f(0.2, 0.75, 0.4), ok); }
    textureStore(screen, id.xy, vec4f(pow(col, vec3f(2.2)), 1.));
}

The solution accumulates into a local variable and writes the output once. Writing outputs[index] += ... inside the loop also gives the right answer, but it reads and writes global memory 256 times instead of once. On a GPU, that difference is usually what separates fast code from slow code.

Stretch: tiling

The naive kernel does 2 reads for every multiply-add. Every thread in a row of the output reads the same row of A from global memory, and every thread in a column reads the same column of B: each value is fetched 256 times. GPUs can do arithmetic far faster than they can fetch from memory, so the naive kernel spends most of its time waiting.

Tiling fixes this with workgroup memory. Each 16 × 16 workgroup walks along k in steps of 16. At each step:

  1. Each thread loads one element of a 16 × 16 tile of A and one of B into shared arrays. Barrier.
  2. Each thread does 16 multiply-adds reading only from the shared tiles. Barrier, so no one overwrites a tile still in use.

Same math, same answer, but each value now comes from global memory 16 times less often.

To measure it, the stretch version runs matmul 16 times per frame with #dispatch_count. Watch the fps counter in the corner, then flip TILED to true once your tiled version works. On a fast GPU you may need to raise the dispatch count before the naive version slows down at all.

#storage inputs array<f32>
#storage outputs array<f32>
#workgroup_count fill 16 16 1
#workgroup_count matmul 16 16 1
#dispatch_count matmul 16

const N = 256u;
const TILE = 16u;
const A = 0u;
const B = N * N;
const TILED = false;     // flip once matmul_tiled works, and compare fps

var<workgroup> tileA: array<array<f32, 16>, 16>;
var<workgroup> tileB: array<array<f32, 16>, 16>;

fn pcg(v: u32) -> u32 {
    let state = v * 747796405u + 2891336453u;
    let word = ((state >> ((state >> 28u) + 4u)) ^ state) * 277803737u;
    return (word >> 22u) ^ word;
}

@compute @workgroup_size(16, 16)
fn fill(@builtin(global_invocation_id) gid: vec3u) {
    let index = gid.y * N + gid.x;
    inputs[A + index] = select(0.0, 1.0, gid.x == gid.y);
    inputs[B + index] = f32(pcg(index) >> 8u) / 16777216.0;
}

fn matmul_naive(row: u32, col: u32) -> f32 {
    var acc = 0.0;
    for (var k = 0u; k < N; k++) {
        acc += inputs[A + row * N + k] * inputs[B + k * N + col];
    }
    return acc;
}

// Every thread in the workgroup must call this (it contains barriers)
fn matmul_tiled(row: u32, col: u32, lid: vec2u) -> f32 {
    var acc = 0.0;
    for (var t = 0u; t < N; t += TILE) {
        // TODO: load tileA[lid.y][lid.x] = A[row][t + lid.x]
        //       and   tileB[lid.y][lid.x] = B[t + lid.y][col], then barrier
        // TODO: acc += tileA[lid.y][k] * tileB[k][lid.x] for k in 0..TILE
        // TODO: barrier before the next tile overwrites these
    }
    return acc;
}

@compute @workgroup_size(16, 16)
fn matmul(@builtin(global_invocation_id) gid: vec3u,
          @builtin(local_invocation_id) lid: vec3u) {
    var acc = 0.0;
    if (TILED) { acc = matmul_tiled(gid.y, gid.x, lid.xy); }
    else { acc = matmul_naive(gid.y, gid.x); }
    outputs[gid.y * N + gid.x] = acc;
}

@compute @workgroup_size(16, 16)
fn main_image(@builtin(global_invocation_id) id: vec3u) {
    let size = textureDimensions(screen);
    if (id.x >= size.x || id.y >= size.y) { return; }

    let side = min(size.x / 2u, size.y);
    let panel = id.x / side;
    let local = vec2u(id.x % side, id.y);
    if (panel > 1u || local.y >= side) {
        textureStore(screen, id.xy, vec4f(0.004, 0.004, 0.008, 1.));
        return;
    }
    let ij = min(local * N / side, vec2u(N - 1u));
    let index = ij.y * N + ij.x;
    let c = outputs[index];
    let ok = abs(c - inputs[B + index]) < 1e-5;

    var col = vec3f(c);
    if (panel == 1u) { col = select(vec3f(0.9, 0.2, 0.15), vec3f(0.2, 0.75, 0.4), ok); }
    textureStore(screen, id.xy, vec4f(pow(col, vec3f(2.2)), 1.));
}
#storage inputs array<f32>
#storage outputs array<f32>
#workgroup_count fill 16 16 1
#workgroup_count matmul 16 16 1
#dispatch_count matmul 16

const N = 256u;
const TILE = 16u;
const A = 0u;
const B = N * N;
const TILED = true;

var<workgroup> tileA: array<array<f32, 16>, 16>;
var<workgroup> tileB: array<array<f32, 16>, 16>;

fn pcg(v: u32) -> u32 {
    let state = v * 747796405u + 2891336453u;
    let word = ((state >> ((state >> 28u) + 4u)) ^ state) * 277803737u;
    return (word >> 22u) ^ word;
}

@compute @workgroup_size(16, 16)
fn fill(@builtin(global_invocation_id) gid: vec3u) {
    let index = gid.y * N + gid.x;
    inputs[A + index] = select(0.0, 1.0, gid.x == gid.y);
    inputs[B + index] = f32(pcg(index) >> 8u) / 16777216.0;
}

fn matmul_naive(row: u32, col: u32) -> f32 {
    var acc = 0.0;
    for (var k = 0u; k < N; k++) {
        acc += inputs[A + row * N + k] * inputs[B + k * N + col];
    }
    return acc;
}

fn matmul_tiled(row: u32, col: u32, lid: vec2u) -> f32 {
    var acc = 0.0;
    for (var t = 0u; t < N; t += TILE) {
        // Each thread fetches one element of each tile
        tileA[lid.y][lid.x] = inputs[A + row * N + (t + lid.x)];
        tileB[lid.y][lid.x] = inputs[B + (t + lid.y) * N + col];
        workgroupBarrier();

        // 16 multiply-adds from fast shared memory
        for (var k = 0u; k < TILE; k++) {
            acc += tileA[lid.y][k] * tileB[k][lid.x];
        }
        workgroupBarrier();
    }
    return acc;
}

@compute @workgroup_size(16, 16)
fn matmul(@builtin(global_invocation_id) gid: vec3u,
          @builtin(local_invocation_id) lid: vec3u) {
    var acc = 0.0;
    if (TILED) { acc = matmul_tiled(gid.y, gid.x, lid.xy); }
    else { acc = matmul_naive(gid.y, gid.x); }
    outputs[gid.y * N + gid.x] = acc;
}

@compute @workgroup_size(16, 16)
fn main_image(@builtin(global_invocation_id) id: vec3u) {
    let size = textureDimensions(screen);
    if (id.x >= size.x || id.y >= size.y) { return; }

    let side = min(size.x / 2u, size.y);
    let panel = id.x / side;
    let local = vec2u(id.x % side, id.y);
    if (panel > 1u || local.y >= side) {
        textureStore(screen, id.xy, vec4f(0.004, 0.004, 0.008, 1.));
        return;
    }
    let ij = min(local * N / side, vec2u(N - 1u));
    let index = ij.y * N + ij.x;
    let c = outputs[index];
    let ok = abs(c - inputs[B + index]) < 1e-5;

    var col = vec3f(c);
    if (panel == 1u) { col = select(vec3f(0.9, 0.2, 0.15), vec3f(0.2, 0.75, 0.4), ok); }
    textureStore(screen, id.xy, vec4f(pow(col, vec3f(2.2)), 1.));
}

With matmul, softmax and one more reduction (layer norm), you have the kernels for the forward pass of a small transformer.

After this Where to go next

Debugging without print

WGSL directives

Real WGSL has a few directives at the top of a file. They’re different from this playground’s # lines, which are rewritten before WGSL sees the code.

enable f16;                                      // optional hardware feature
requires readonly_and_readwrite_storage_textures; // newer language feature
diagnostic(off, derivative_uniformity);          // tune a compiler check

Attributes such as @compute, @workgroup_size and @builtin are not directives. They attach to one function or variable.

Your own WebGPU program

This playground hides the host-side setup: requesting a device, creating buffers and bind groups, building pipelines, and reading results back. That setup is where much of the real-world difficulty lives, and it’s the natural next step. A good first program is deliberately boring: add two arrays on the GPU, read the result back, and compare it with the same sum computed in JavaScript. It touches every piece of plumbing with nothing visual to distract you. Then port your softmax and matmul over, and check them against CPU results.

For reference as you go: WebGPU Fundamentals explains the host API step by step, the WGSL specification is the final word on the language, and compute.toys has hundreds of shaders written in the same style as this page.