diff --git a/docs/_includes/mesh-band.html b/docs/_includes/mesh-band.html
new file mode 100644
index 0000000..37cfef4
--- /dev/null
+++ b/docs/_includes/mesh-band.html
@@ -0,0 +1,23 @@
+{%- comment -%}
+Mesh-gradient band, drawn by the mesh-band.js canvas runner.
+
+Parameters:
+ colors comma-separated spot colours (palette defaults: green, rust, python blue)
+ weights comma-separated numbers, one per spot; the largest becomes full opacity and
+ the rest are scaled against it. Pass the note's headline ratios so the
+ picture says what the note says.
+ caption one sentence for the figcaption, ideally naming what the weights are.
+ still optional full-size PNG URL: shown instead of the canvas when scripting is off.
+{%- endcomment -%}
+
+
+
+
+ {{ include.caption }}
+ Field: Paper Shaders mesh gradient (Apache-2.0), palette and weights from this note.
+
+
diff --git a/docs/_layouts/default.html b/docs/_layouts/default.html
index c3e0a5b..18305f8 100644
--- a/docs/_layouts/default.html
+++ b/docs/_layouts/default.html
@@ -8,6 +8,7 @@
{% feed_meta %}
{% if page.profiles %}{% endif %}
+ {% if page.mesh_band %}{% endif %}
{% if page.math %}
diff --git a/docs/assets/figures/mage-004/matmul-layouts-mobile.svg b/docs/assets/figures/mage-004/matmul-layouts-mobile.svg
new file mode 100644
index 0000000..b60d011
--- /dev/null
+++ b/docs/assets/figures/mage-004/matmul-layouts-mobile.svg
@@ -0,0 +1,4882 @@
+
+
+
diff --git a/docs/assets/figures/mage-004/matmul-layouts.png b/docs/assets/figures/mage-004/matmul-layouts.png
new file mode 100644
index 0000000..f3b7f6f
Binary files /dev/null and b/docs/assets/figures/mage-004/matmul-layouts.png differ
diff --git a/docs/assets/figures/mage-004/matmul-layouts.svg b/docs/assets/figures/mage-004/matmul-layouts.svg
new file mode 100644
index 0000000..afe44c6
--- /dev/null
+++ b/docs/assets/figures/mage-004/matmul-layouts.svg
@@ -0,0 +1,4882 @@
+
+
+
diff --git a/docs/assets/figures/mage-004/mesh-band.png b/docs/assets/figures/mage-004/mesh-band.png
new file mode 100644
index 0000000..64a234a
Binary files /dev/null and b/docs/assets/figures/mage-004/mesh-band.png differ
diff --git a/docs/assets/figures/mage-007/mesh-band.png b/docs/assets/figures/mage-007/mesh-band.png
new file mode 100644
index 0000000..7298492
Binary files /dev/null and b/docs/assets/figures/mage-007/mesh-band.png differ
diff --git a/docs/assets/mesh-band.js b/docs/assets/mesh-band.js
new file mode 100644
index 0000000..98879a9
--- /dev/null
+++ b/docs/assets/mesh-band.js
@@ -0,0 +1,347 @@
+/* A mesh-gradient band, driven by the numbers in the note it sits under.
+ *
+ * The fragment shader is Paper Shaders' mesh gradient
+ * (https://github.com/paper-design/shaders, Apache-2.0, (c) Paper Design), kept
+ * verbatim; the vertex shader, the palette, the weights and the lifecycle
+ * (visibility, reduced motion, DPR) are local.
+ *
+ * Markup contract, in _includes/mesh-band.html:
+ *
+ *
+ *
+ * weights become each spot's opacity, normalised against the largest: in these
+ * pages a weight is a measured ratio from the note, so the brightest spot is the
+ * largest one. The canvas stays hidden until a frame has been drawn, so a device
+ * without WebGL 2 shows the CSS fallback instead of an empty box.
+ */
+(() => {
+ "use strict";
+
+ const FRAG = `#version 300 es
+precision highp float;
+
+uniform vec2 u_resolution;
+uniform float u_time;
+uniform vec4 u_colors[10];
+uniform float u_colorsCount;
+uniform float u_distortion;
+uniform float u_swirl;
+uniform float u_grainMixer;
+uniform float u_grainOverlay;
+uniform float u_intensity;
+uniform float u_opacityGain;
+
+in vec2 v_objectUV;
+out vec4 fragColor;
+
+vec2 rotate(vec2 uv, float th) {
+ return mat2(cos(th), sin(th), -sin(th), cos(th)) * uv;
+}
+
+float hash21(vec2 p) {
+ p = fract(p * vec2(0.3183099, 0.3678794)) + 0.1;
+ p += dot(p, p + 19.19);
+ return fract(p.x * p.y);
+}
+
+float valueNoise(vec2 st) {
+ vec2 i = floor(st);
+ vec2 f = fract(st);
+ float a = hash21(i);
+ float b = hash21(i + vec2(1.0, 0.0));
+ float c = hash21(i + vec2(0.0, 1.0));
+ float d = hash21(i + vec2(1.0, 1.0));
+ vec2 u = f * f * (3.0 - 2.0 * f);
+ float x1 = mix(a, b, u.x);
+ float x2 = mix(c, d, u.x);
+ return mix(x1, x2, u.y);
+}
+
+float noise(vec2 n, vec2 seedOffset) {
+ return valueNoise(n + seedOffset);
+}
+
+vec2 getPosition(int i, float t) {
+ float a = float(i) * .37;
+ float b = .6 + fract(float(i) / 3.) * .9;
+ float c = .8 + fract(float(i + 1) / 4.);
+
+ float x = sin(t * b + a);
+ float y = cos(t * c + a * 1.5);
+
+ return .5 + .5 * vec2(x, y);
+}
+
+void main() {
+ vec2 uv = v_objectUV;
+ uv += .5;
+ vec2 grainUV = uv * 1000.;
+
+ float mixerGrain = 0.;
+ if (u_grainMixer > 0.) {
+ mixerGrain = .4 * u_grainMixer * (noise(grainUV, vec2(0.)) - .5);
+ }
+
+ const float firstFrameOffset = 41.5;
+ float t = .5 * (u_time + firstFrameOffset);
+
+ float radius = smoothstep(0., 1., length(uv - .5));
+ float center = 1. - radius;
+ for (float i = 1.; i <= 2.; i++) {
+ uv.x += u_distortion * center / i * sin(t + i * .4 * smoothstep(.0, 1., uv.y)) * cos(.2 * t + i * 2.4 * smoothstep(.0, 1., uv.y));
+ uv.y += u_distortion * center / i * cos(t + i * 2. * smoothstep(.0, 1., uv.x));
+ }
+
+ vec2 uvRotated = uv;
+ uvRotated -= vec2(.5);
+ float angle = 3. * u_swirl * radius;
+ uvRotated = rotate(uvRotated, -angle);
+ uvRotated += vec2(.5);
+
+ vec3 color = vec3(0.);
+ float opacity = 0.;
+ float totalWeight = 0.;
+
+ for (int i = 0; i < 10; i++) {
+ if (i >= int(u_colorsCount)) break;
+
+ vec2 pos = getPosition(i, t) + mixerGrain;
+ vec3 colorFraction = u_colors[i].rgb * u_colors[i].a;
+ float opacityFraction = u_colors[i].a;
+
+ float dist = length(uvRotated - pos);
+
+ dist = pow(dist, 3.5);
+ float weight = 1. / (dist + 1e-3);
+ color += colorFraction * weight;
+ opacity += opacityFraction * weight;
+ totalWeight += weight;
+ }
+
+ color /= max(1e-4, totalWeight);
+ opacity /= max(1e-4, totalWeight);
+
+ // Local, so a bright mesh sits on a dark page as a bloom rather than a wash.
+ // Paper's shader averages the spot colours un-premultiplied; on the notebook's
+ // #101217 backdrop that reads as white. These two knobs scale the result.
+ color *= u_intensity;
+ opacity *= u_opacityGain;
+
+ if (u_grainOverlay > 0.) {
+ float grainOverlay = valueNoise(rotate(grainUV, 1.) + vec2(3.));
+ grainOverlay = mix(grainOverlay, valueNoise(rotate(grainUV, 2.) + vec2(-1.)), .5);
+ grainOverlay = pow(grainOverlay, 1.3);
+
+ float grainOverlayV = grainOverlay * 2. - 1.;
+ vec3 grainOverlayColor = vec3(step(0., grainOverlayV));
+ float grainOverlayStrength = u_grainOverlay * abs(grainOverlayV);
+ grainOverlayStrength = pow(grainOverlayStrength, .8);
+ color = mix(color, grainOverlayColor, .35 * grainOverlayStrength);
+
+ opacity += .5 * grainOverlayStrength;
+ }
+ opacity = clamp(opacity, 0., 1.);
+
+ fragColor = vec4(color, opacity);
+}
+`;
+
+ const VERT = `#version 300 es
+precision highp float;
+
+uniform float u_aspect;
+in vec2 a_position;
+out vec2 v_objectUV;
+
+void main() {
+ v_objectUV = vec2(a_position.x * u_aspect, a_position.y);
+ gl_Position = vec4(a_position * 2.0, 0.0, 1.0);
+}
+`;
+
+ const BACKDROP = [0x10 / 255, 0x12 / 255, 0x17 / 255, 1.0];
+
+ function rgb(hex) {
+ const v = hex.trim().replace("#", "");
+ return [
+ parseInt(v.slice(0, 2), 16) / 255,
+ parseInt(v.slice(2, 4), 16) / 255,
+ parseInt(v.slice(4, 6), 16) / 255,
+ ];
+ }
+
+ function compile(gl, type, source) {
+ const shader = gl.createShader(type);
+ gl.shaderSource(shader, source);
+ gl.compileShader(shader);
+ if (!gl.getShaderParameter(shader, gl.COMPILE_STATUS)) {
+ const log = gl.getShaderInfoLog(shader);
+ gl.deleteShader(shader);
+ throw new Error(log || "shader failed to compile");
+ }
+ return shader;
+ }
+
+ function program(gl) {
+ const p = gl.createProgram();
+ gl.attachShader(p, compile(gl, gl.VERTEX_SHADER, VERT));
+ gl.attachShader(p, compile(gl, gl.FRAGMENT_SHADER, FRAG));
+ gl.linkProgram(p);
+ if (!gl.getProgramParameter(p, gl.LINK_STATUS)) {
+ throw new Error(gl.getProgramInfoLog(p) || "program failed to link");
+ }
+ return p;
+ }
+
+ function start(band) {
+ const canvas = band.querySelector("canvas");
+ const gl = canvas.getContext("webgl2", {
+ alpha: true,
+ antialias: false,
+ depth: false,
+ powerPreference: "low-power",
+ preserveDrawingBuffer: true, // lets a still be pulled out of the canvas
+ });
+ if (!gl) return;
+
+ let prog;
+ try {
+ prog = program(gl);
+ } catch (error) {
+ console.warn("mesh-band:", error.message);
+ return;
+ }
+
+ gl.useProgram(prog);
+
+ const buffer = gl.createBuffer();
+ gl.bindBuffer(gl.ARRAY_BUFFER, buffer);
+ gl.bufferData(
+ gl.ARRAY_BUFFER,
+ new Float32Array([-0.5, -0.5, 0.5, -0.5, -0.5, 0.5, 0.5, 0.5]),
+ gl.STATIC_DRAW
+ );
+ const aPosition = gl.getAttribLocation(prog, "a_position");
+ gl.enableVertexAttribArray(aPosition);
+ gl.vertexAttribPointer(aPosition, 2, gl.FLOAT, false, 0, 0);
+
+ const colors = (band.dataset.colors || "#91dbba,#c9b2ff,#93caff")
+ .split(",")
+ .slice(0, 10)
+ .map(rgb);
+ const weights = (band.dataset.weights || "")
+ .split(",")
+ .map(Number)
+ .filter((w) => Number.isFinite(w) && w > 0);
+ const peak = weights.length ? Math.max(...weights) : 1;
+ while (weights.length < colors.length) weights.push(peak);
+
+ const spots = new Float32Array(40);
+ colors.forEach((c, i) => {
+ const alpha = Math.min(1, weights[i] / peak);
+ spots.set([c[0], c[1], c[2], alpha], i * 4);
+ });
+
+ const u = {
+ resolution: gl.getUniformLocation(prog, "u_resolution"),
+ aspect: gl.getUniformLocation(prog, "u_aspect"),
+ time: gl.getUniformLocation(prog, "u_time"),
+ colors: gl.getUniformLocation(prog, "u_colors"),
+ colorsCount: gl.getUniformLocation(prog, "u_colorsCount"),
+ distortion: gl.getUniformLocation(prog, "u_distortion"),
+ swirl: gl.getUniformLocation(prog, "u_swirl"),
+ grainMixer: gl.getUniformLocation(prog, "u_grainMixer"),
+ grainOverlay: gl.getUniformLocation(prog, "u_grainOverlay"),
+ intensity: gl.getUniformLocation(prog, "u_intensity"),
+ opacityGain: gl.getUniformLocation(prog, "u_opacityGain"),
+ };
+
+ gl.uniform4fv(u.colors, spots);
+ gl.uniform1f(u.colorsCount, colors.length);
+ gl.uniform1f(u.distortion, Number(band.dataset.distortion ?? 0.8));
+ gl.uniform1f(u.swirl, Number(band.dataset.swirl ?? 0.55));
+ gl.uniform1f(u.grainMixer, Number(band.dataset.grainMixer ?? 0.05));
+ gl.uniform1f(u.grainOverlay, Number(band.dataset.grainOverlay ?? 0.04));
+ gl.uniform1f(u.intensity, Number(band.dataset.intensity ?? 0.8));
+ gl.uniform1f(u.opacityGain, Number(band.dataset.opacity ?? 0.6));
+ gl.enable(gl.BLEND);
+ gl.blendFunc(gl.SRC_ALPHA, gl.ONE_MINUS_SRC_ALPHA);
+ gl.clearColor(...BACKDROP);
+
+ function resize() {
+ const dpr = Math.min(window.devicePixelRatio || 1, 1.5);
+ const width = Math.max(1, Math.round(canvas.clientWidth * dpr));
+ const height = Math.max(1, Math.round(canvas.clientHeight * dpr));
+ if (canvas.width !== width || canvas.height !== height) {
+ canvas.width = width;
+ canvas.height = height;
+ }
+ gl.viewport(0, 0, width, height);
+ gl.uniform2f(u.resolution, width, height);
+ gl.uniform1f(u.aspect, canvas.clientWidth / Math.max(1, canvas.clientHeight));
+ }
+
+ function frame(seconds) {
+ resize();
+ gl.clear(gl.COLOR_BUFFER_BIT);
+ gl.uniform1f(u.time, seconds);
+ gl.drawArrays(gl.TRIANGLE_STRIP, 0, 4);
+ }
+
+ const reduced = window.matchMedia("(prefers-reduced-motion: reduce)");
+ let visible = true;
+ let raf = 0;
+ let t0 = 0;
+
+ function loop(now) {
+ if (!t0) t0 = now;
+ frame((now - t0) / 1000);
+ raf = requestAnimationFrame(loop);
+ }
+
+ function play() {
+ if (raf || reduced.matches || !visible) return;
+ raf = requestAnimationFrame(loop);
+ }
+
+ function pause() {
+ if (!raf) return;
+ cancelAnimationFrame(raf);
+ raf = 0;
+ }
+
+ band.dataset.ready = "true";
+ frame(0); // a valid first paint, before any rAF tick and regardless of visibility
+ if (!reduced.matches) play();
+
+ if ("IntersectionObserver" in window) {
+ new IntersectionObserver(
+ (entries) => {
+ visible = entries.some((e) => e.isIntersecting);
+ if (visible) play();
+ else pause();
+ },
+ { rootMargin: "100px" }
+ ).observe(band);
+ }
+ reduced.addEventListener("change", () => {
+ if (reduced.matches) {
+ pause();
+ frame(0);
+ } else play();
+ });
+ window.addEventListener("resize", () => {
+ if (reduced.matches) frame(0);
+ });
+ }
+
+ function boot() {
+ document.querySelectorAll(".mesh-band").forEach(start);
+ }
+
+ if (document.readyState === "loading") {
+ document.addEventListener("DOMContentLoaded", boot);
+ } else {
+ boot();
+ }
+})();
diff --git a/docs/assets/notebook.css b/docs/assets/notebook.css
index cec8822..c36696a 100644
--- a/docs/assets/notebook.css
+++ b/docs/assets/notebook.css
@@ -152,6 +152,32 @@ mjx-container[display] { overflow-x: auto; overflow-y: hidden; padding-block: 8p
.profile-plot img { display: block; width: 100%; height: auto; }
.profile-plot figcaption, .profile-limit { font-size: 12px; line-height: 1.7; color: var(--muted); }
.profile-plot figcaption p { margin: 8px 0; }
+
+/* Mesh-gradient band: a canvas hero for notes whose argument is about memory.
+ The CSS backdrop shows until a frame has been drawn, so a device without
+ WebGL 2 (or without scripting) still gets an intentional panel. */
+.mesh-band {
+ position: relative;
+ margin: 26px 0 30px;
+ border: 1px solid var(--rule);
+ border-radius: 10px;
+ overflow: hidden;
+ aspect-ratio: 740 / 220;
+ background:
+ radial-gradient(120% 160% at 18% 28%, rgba(145, 219, 186, .22), transparent 62%),
+ radial-gradient(110% 150% at 78% 62%, rgba(201, 178, 255, .20), transparent 60%),
+ radial-gradient(90% 130% at 52% 92%, rgba(147, 202, 255, .16), transparent 58%),
+ var(--paper);
+}
+.mesh-band canvas { display: block; width: 100%; height: 100%; opacity: 0; transition: opacity .7s ease; }
+.mesh-band[data-ready] canvas { opacity: 1; }
+.mesh-band noscript img { display: block; width: 100%; height: auto; }
+.mesh-band figcaption { padding: 10px 14px 12px; border-top: 1px solid var(--rule); background: var(--surface); font-size: 12px; line-height: 1.7; color: var(--muted); }
+.mesh-band figcaption p { margin: 0 0 6px; }
+.mesh-band-credit { display: block; margin-top: 6px; font: 11px/1.6 var(--mono); color: #6f7b8f; }
+.mesh-band-credit a { color: #6f7b8f; }
+.mesh-band-credit a:hover { color: var(--muted); }
+@media (max-width: 520px) { .mesh-band { aspect-ratio: 320 / 200; } }
.profile-scale { color: var(--ink); }
.profile-links { display: flex; flex-wrap: wrap; gap: 16px; margin-top: 11px; }
.profile-links a { text-decoration: none; }
diff --git a/docs/experiments/mage-004.md b/docs/experiments/mage-004.md
index 9d5c4a6..baa56db 100644
--- a/docs/experiments/mage-004.md
+++ b/docs/experiments/mage-004.md
@@ -2,177 +2,135 @@
title: What the compiler chose
permalink: /experiments/mage-004/
eyebrow: "Field note 004 / Mathematics on a GPU"
-description: "Handing the kernel decisions to a tile compiler — where it beat the hand-written kernels, where it lost, and why most of the apparent difference turned out to be the cost of starting work rather than the cost of doing it."
+description: "A tile compiler makes the memory decisions for you. On five FP32 operations it wins one, draws one and loses three — and the losses are exactly where data reuse is highest, which is where layout decides the answer."
math: true
+mesh_band: true
---
-
-The [earlier notes]({{ '/experiments/mage-003/' | relative_url }}) wrote these kernels by hand: each
-thread owned a fixed block of results, and the width of every load, the layout of shared memory and
-the placement of barriers were deliberate choices with measured effects.
-
-This note gives those decisions to a compiler. A [cuTile Rust](https://github.com/NVlabs/cutile-rs)
-kernel is a single-threaded program over *tiles* — blocks of data — and the compiler decides how many
-warps receive each tile, which values stay in registers, when loads widen to 128 bits, and when a
-tensor-core instruction replaces a multiply. On the same five operations it is competitive: bias +
-GELU 8.19 µs against the hand-written 11.0, layer normalization 10.76 against 10.05. It trails where
-reuse is highest — matrix multiply 131.56 against 80.00, triangle contraction 116.08 against 80.1 —
-and by 3.4× on neighbor aggregation, whose irregular gathers have no safe expression in the tile
-model and run through raw device pointers instead.
-
-Two measurement facts matter more than the ranking. Timing that waits for each call prices the cost
-of *starting* work: 14–23 µs per call against 2–3 µs for the hand-written kernel, which vanishes
-when ten calls share one measurement — bias + GELU goes from 25.76 µs to 8.50 µs against an 8.19 µs
-kernel. And tile shape is a memory decision: the same arithmetic measured 738 µs at the tutorial's
-16×16×8 tile and 131.56 µs at 32×128×32, and the compiler's own autotuner found a shape 27% faster
-than the best of twelve hand-picked ones.
-
-## The five operations, two views
-
-| Operation | PyTorch | Triton | cuTile Rust | cuda-oxide Rust |
-| --- | ---: | ---: | ---: | ---: |
-| Matrix multiplication 1024³ | 56.04 | 83.06 | 131.56 | **80.00** |
-| Bias + GELU 4096×768 | 15.88 | 7.73 | **8.19** | 11.0 |
-| LayerNorm 4096×768 | 11.42 | 8.15 | 10.76 | 10.05 |
-| Triangle contraction 128×32 | 28.64 | 102.19 | 116.08 | 80.1 |
-| Neighbor aggregation 4096×64×65536 | 67.32 | 7.97 | 35.37 | 10.3 |
-
-*Microseconds of GPU kernel time per operation, from one Nsight Systems capture of 100 launches
-each. Lower is better. PyTorch's bias + GELU is two kernels, its triangle contraction three and
-its neighbor aggregation four; every other cell is one.*
+{% include mesh-band.html
+ colors="#93caff,#91dbba,#c9b2ff,#e0a08a"
+ weights="0.74,1,0.62,0.90"
+ still="/assets/figures/mage-004/mesh-band.png"
+ alt="A dark field with four soft spots of colour — blue, green, lavender and warm sand — scaled by how competitive each implementation is."
+ caption="The four spots are the four implementations in the table below, opacity set by the geometric mean of their kernel times relative to the best implementation on each operation: Triton brightest, cuTile Rust dimmest." %}**The claim.** Give a compiler the job of deciding how a GPU kernel places its data, and it will do
+a good job where the reuse is low and a worse one where the reuse is high. On five FP32 operations,
+the [cuTile Rust](https://github.com/NVlabs/cutile-rs) tile kernels beat the hand-written ones on
+bias + GELU ($8.19\ \mu s$ against $11.0$), draw on layer normalization ($10.76$ against $10.05$),
+and lose on matrix multiply ($131.56$ against $80.00$), triangle contraction ($116.08$ against
+$80.1$) and neighbor aggregation ($35.37$ against $10.3$).
+
+Two things support that claim, and both are about memory rather than arithmetic.
+
+## The first argument: reuse is a layout problem
+
+Matrix multiplication is the clearest case, because the arithmetic is free: every implementation
+computes $C = A B$ with the same multiply-adds, and the only question is how often each value has to
+be fetched. A thread that computes one output element reads one value of $A$ and one of $B$ per
+multiply-add. A thread that computes a $4 \times 4$ block of outputs reads four of $A$ and four of
+$B$ and performs sixteen multiply-adds, so shared-memory reads per multiply-add fall from
+$\frac{2}{1}$ to $\frac{8}{16} = 0.5$, and reading those values as 128-bit quads takes it to
+$0.125$ — one instruction fetching the four values a thread needs.
-
-
+
+
-
GPU kernel time per operation, from separate Nsight Systems captures of 100 launches each. Each row has its own scale. The tile kernels win bias + GELU, draw on layer normalization, and trail on the other three.
+
$C = A B$ at $1024^3$ in FP32. On the left the compiler decides how the tile is threaded and how the loads are widened; on the right the kernel author does. The transposed $A$ with a row stride of 68 exists so that each thread's four rows are contiguous — one 128-bit read instead of four scattered 32-bit ones.
-
- Read the plotted values (µs of GPU kernel time)
-
-
GPU kernel time per operation and implementation
-
Operation
PyTorch
Triton
cuTile Rust
cuda-oxide Rust
-
-
Matrix multiplication 1024³
56.04
83.06
131.56
80.00
-
Bias + GELU 4096×768
15.88
7.73
8.19
11.0
-
LayerNorm 4096×768
11.42
8.15
10.76
10.05
-
Triangle contraction 128×32
28.64
102.19
116.08
80.1
-
Neighbor aggregation 4096×64×65536
67.32
7.97
35.37
10.3
-
-
-
-The same run, timed the way most comparisons are timed — one measurement around each call, waiting
-for it to finish:
+The tile compiler is not blind to this. It widens loads too, and it partitions the output into
+$32 \times 128$ tiles with $32 \times 32$ operands per contraction step. What it cannot know is that
+this particular problem wants a different arrangement: the shapes that win are the ones whose
+default layout happens to be close to good, and matrix multiply is not one of them. Its kernel is
+$1.6\times$ slower, and nothing about the arithmetic explains that — only where the operands sit
+when the multiply instruction issues.
+
+## The second argument: the clock measures two things
+
+Time around a call is not the kernel's time. It is submission plus the kernel:
+
+$$\text{span} \;=\; \underbrace{t_{\text{submit}}}_{\text{host builds and queues}} \;+\; \underbrace{t_{\text{kernel}}}_{\text{device executes}} \;+\; \text{idle}$$
+
+The tile runtime submits lazily, and one awaited call costs $14\text{–}23\ \mu s$ of host time
+against $2\text{–}3\ \mu s$ for the hand-written kernel. So a comparison that waits for every call
+is comparing submission paths, and the tile kernels look far worse than they are: bias + GELU
+measures $29.97\ \mu s$ around a kernel that runs in $8.19\ \mu s$.
+
+| Operation | awaited | queued in tens | kernel only |
+| --- | ---: | ---: | ---: |
+| Matrix multiply $1024^3$ | 150.53 | 127.07 | 131.56 |
+| Bias + GELU $4096\times768$ | 25.76 | 8.50 | 8.19 |
+| LayerNorm $4096\times768$ | 30.62 | 10.64 | 10.76 |
+| Triangle contraction $128\times32$ | 136.19 | 121.80 | 116.08 |
+| Neighbor aggregation $4096\times64\times65536$ | 49.28 | 33.28 | 35.37 |
+
+*Microseconds, mean of 100 launches. Queued, every row lands on its kernel time; replaying ten
+launches from a recorded CUDA graph does the same ($7.77\ \mu s$ for bias + GELU).*
+
+## Evidence
+
+Kernel time, one Nsight Systems capture of 100 launches per implementation. Every output was checked
+against PyTorch with TF32 disabled before any timing was believed; the worst error is
+$1.5\times10^{-5}$ on the matrix multiply.
| Operation | PyTorch | Triton | cuTile Rust | cuda-oxide Rust |
| --- | ---: | ---: | ---: | ---: |
-| Matrix multiplication 1024³ | 52.40 | 93.50 | 170.79 | 82.6 |
-| Bias + GELU 4096×768 | 33.02 | 27.93 | 29.97 | 13.3 |
-| LayerNorm 4096×768 | 21.68 | 23.52 | 37.22 | 12.9 |
-| Triangle contraction 128×32 | 61.81 | 99.96 | 157.84 | 83.5 |
-| Neighbor aggregation 4096×64×65536 | 95.07 | 27.00 | 61.22 | 12.9 |
+| Matrix multiplication $1024^3$ | 56.04 | 83.06 | 131.56 | **80.00** |
+| Bias + GELU $4096\times768$ | 15.88 | 7.73 | **8.19** | 11.0 |
+| LayerNorm $4096\times768$ | 11.42 | 8.15 | 10.76 | 10.05 |
+| Triangle contraction $128\times32$ | 28.64 | 102.19 | 116.08 | 80.1 |
+| Neighbor aggregation $4096\times64\times65536$ | 67.32 | 7.97 | 35.37 | 10.3 |
-*Microseconds around each call, mean of three rounds of 100 samples with the implementations
-rotated. Every implementation is measured the same way here, and the ordering is different.*
+*Microseconds of GPU kernel time per operation, summed over however many kernels each operation
+launches. Lower is better.*
-## Where the difference goes
+Layer normalization is worth reading twice. The tile kernel pads a $768$-wide row to $1024$ —
+tile dimensions must be powers of two — so a third of its lanes do nothing, and it still draws with
+the hand-written kernel. Out of padding, it would win.
-The two tables disagree, and the disagreement is measurable rather than mysterious. Taking the
-same launches and queueing ten of them behind one measurement:
+## Tile shape is part of the layout
-| Operation | awaited | queued by ten | kernel |
-| --- | ---: | ---: | ---: |
-| Matrix multiplication 1024³ | 150.53 | 127.07 | 131.56 |
-| Bias + GELU 4096×768 | 25.76 | 8.50 | 8.19 |
-| LayerNorm 4096×768 | 30.62 | 10.64 | 10.76 |
-| Triangle contraction 128×32 | 136.19 | 121.80 | 116.08 |
-| Neighbor aggregation 4096×64×65536 | 49.28 | 33.28 | 35.37 |
-
-Queued, every row lands on its kernel time — bias + GELU and layer norm within 0.5 µs, and matrix
-multiply below its own kernel time because submission overlaps execution. Recording ten launches
-into a replayable graph does the same: bias + GELU comes out at 7.77 µs against an 8.19 µs kernel.
-Timed the same way, Triton improves by 6–7 µs per call and PyTorch by about 1 µs, so with all
-three launch paths matched the wins and losses are the ones in the first table and nothing else.
-
-*Caveat: the cuda-oxide column above is the third note's retained measurement, not taken in the
-same session. The kernel track has since published [mage-006]({{ '/experiments/mage-006/' | relative_url }})
-with lower numbers of its own, so the two Rust columns should not be subtracted from each other.*
-
-## Tile shape is a memory decision
-
-The matrix multiply started from the upstream tutorial's 16 × 16 × 8 tile and was slow. Twelve
-hand-picked configurations took it from 738 µs to 201 µs, and the shape was not monotone in any
-direction: deepening the contraction step cost 42% at a 16 × 16 tile but only 6% at 64 × 64, and
-128 × 128 × 8 was 2.3× *slower* than 128 × 64 × 8 — the signature of a register or occupancy
-limit rather than arithmetic.
-
-| Tile | Span | Tile | Span |
-| --- | ---: | --- | ---: |
-| 16×16×8 | 738.0 | 64×128×8 | 213.9 |
-| 32×32×32 | 488.8 | 128×64×8 | 201.5 |
-| 64×64×32 | 255.4 | 128×128×8 | 469.8 |
-
-Then the library's autotuner searched the same space properly — 36 candidates, each validated
-before timing, kernel time alone, in 17 seconds — and chose **32 × 128 × 32**. Checked in one
-session, that shape measured 137.28 µs against the hand-picked 128 × 64 × 8's 189.44 µs. Every
-number in this note was re-measured with it afterwards, which is where the matrix multiply's
-131.56 µs comes from.
+The same kernel and the same arithmetic, twelve tile shapes: $738\ \mu s$ at the tutorial's
+$16\times16\times8$, $201\ \mu s$ at $128\times64\times8$, and $469\ \mu s$ at $128\times128\times8$
+— $2.3\times$ *slower* than a smaller tile, which is a register or occupancy cliff rather than
+arithmetic. The library's autotuner then searched the space properly: 36 candidates, each validated
+before timing, and its pick measured $137.28\ \mu s$ against the hand-picked tile's $189.44\ \mu s$.
+Every number in the table above was re-measured with it.
## What the compiler would not let us write
-These cost time to find, and each is silent until the compiler or the assembler runs:
-
-- **Tile dimensions must be powers of two.** A row of 768 columns is rejected with
- `failed to compile Tile IR program` and no further detail, so layer norm pads each row to 1024
- and divides by the true width — 33% of its lanes are wasted, and it still draws level.
-- **A partition load indexed by a loop variable does not vary.** The triangle kernel loaded every
- channel through a `for` loop and read the first channel every time; the channel has to come from
- the grid axes instead.
-- **A `Tile<..>` written inside an expression is not rewritten** by the entry macro, so it must
- appear only in a `let` annotation.
-- **`convert_scalar` has no `u32` → `i32`.** The CSR arrays are uploaded as `i32`, which the
- host's bounds check makes safe.
-- **Scalar comparison is not a supported operator**, so the edge walk counts its edges and loops
- over the count.
+- **Tile dimensions must be powers of two.** A $768$-wide row is rejected with
+ `failed to compile Tile IR program` and no further detail.
+- **A partition load indexed by a loop variable does not vary.** The triangle kernel read the first
+ channel every time until the channel came from the grid axes.
+- **A `Tile<..>` inside an expression is not rewritten** by the entry macro; it belongs in a `let`
+ annotation.
+- **`convert_scalar` has no `u32` → `i32`**, so the CSR arrays are uploaded as `i32`.
+- **Scalar comparison is not a supported operator**, so the edge walk counts its edges.
The irregular operation needed raw device pointers (`load_ptr_tko`) to read its row pointers and
edge indices. That escape hatch is the honest boundary of the safe tile model, and it is why that
-kernel is the one that trails furthest.
+kernel trails furthest.
## What the numbers do not establish
-- No hardware counters are available here, so occupancy and bandwidth are inferred from ratios,
- not read.
-- Layer norm's 10.76 µs includes the 768 → 1024 row padding; a native 768-wide tile would do less
- work.
-- One capture per operation: a kernel time here carries no interval of its own. The other column
- comes from three rotating rounds.
-- Measured timings are sensitive to a busy device. A repeat that shared the GPU with another
- capture measured three to twenty times larger spans for the same binary, and it was caught only
- because PyTorch's own kernel moved with it.
-- Compilation, transfers and process startup are excluded throughout.
-
-## Open items
-
-- **Replay for the other three kernels.** Only matrix multiply and bias + GELU can be replayed
- from a recorded graph today; each kernel needs its own capture.
-- **Neighbor aggregation stays 3.4× behind.** It gathers about 17 MB of rows, which is
- L2-bandwidth work on this part, and its 35.37 µs is far above what that traffic accounts for.
-- **Lower precision is a separate question.** FP16, BF16 and TF32 have their own error budgets and
- nothing here speaks to them.
+- No hardware counters are available here, so occupancy and bandwidth are inferred from ratios.
+- One capture per operation: a kernel time carries no interval of its own.
+- The cuda-oxide column is the third note's retained measurement, not taken in the same session.
+- Timings are sensitive to a busy device: a repeat that shared the GPU with another capture showed
+ three to twenty times larger spans for the same binary.
-## Reproduction
+## Reproduce
```bash
source scripts/cutile-env.sh
@@ -182,5 +140,5 @@ cd examples/cutile && cargo build --release && cd ../..
--rounds 3 --iterations 100 --warmup 25
```
-`CUTILE_MATMUL_TILE=BM,BN,BK` picks a different tile; `--mode single|batch|graph` picks the launch
-path; `--features tune` adds the autotuner. Run one device experiment at a time.
+`CUTILE_MATMUL_TILE=BM,BN,BK` picks a tile shape, `--mode single|batch|graph` the launch path, and
+`--features tune` adds the autotuner. Run one device experiment at a time.
diff --git a/docs/experiments/mage-007.md b/docs/experiments/mage-007.md
index 8dd3cac..3530c61 100644
--- a/docs/experiments/mage-007.md
+++ b/docs/experiments/mage-007.md
@@ -2,22 +2,42 @@
title: Wider loads helped two kernels and hurt a third
permalink: /experiments/mage-007/
eyebrow: "Field note 007 / Mathematics on a GPU"
-description: "One gap closed, one halved, and one attempt that measured 2.6x slower and was reverted, each before-and-after pair taken in one session."
+description: "A wider load per thread is a trade: fewer instructions issued, more bytes in flight. Applied to three kernels it paid twice and cost once, and the one it cost is the one that was already using every thread the device has."
math: true
+mesh_band: true
---
-Neighbor aggregation sums, for each output row, the input rows named by that row's edges. Each
-thread read one float per edge — a 32-bit load, one edge at a time — so a row with 16 edges issued
-16 narrow loads in a dependent chain:
+{% include mesh-band.html
+ colors="#91dbba,#93caff,#e08a8a"
+ weights="1.34,1.37,0.38"
+ still="/assets/figures/mage-007/mesh-band.png"
+ alt="A dark field with three soft spots of colour — green, blue and muted red — the red one much dimmer than the other two."
+ caption="The three spots are this note's three changes, opacity set by the speed-up each measured: neighbour aggregation (1.34×, green), LayerNorm (1.37×, blue) and the rejected GELU variant (0.38×, red)." %}**The claim.** Making each thread load four values at once instead of one is a memory decision, not
+an optimisation: it cuts the number of load instructions by four and puts $4\times$ the bytes in
+flight per instruction. It pays when a kernel is short of instructions to issue. It costs when the
+grid already fills the machine and the only thing hiding memory latency is how many threads are
+running. The same change, applied to three kernels, gave $1.34\times$ and $1.37\times$ speedups,
+and then a $2.6\times$ slowdown.
+
+## Argument one: a gather wants fewer, wider loads
+
+Neighbor aggregation computes one output row per input row:
+
+$$y_{i,f} \;=\; \sum_{e \in \mathrm{row}(i)} w_e \, x_{i_e,\, f}$$
+
+The edge list is data — row $i$ owns `rowptr[i]..rowptr[i+1]` entries of `indices` — so each thread
+walks its own row's edges and asks global memory for one value of $x$ per edge. At
+$4096 \times 64$ with $65536$ edges each thread issues 16 loads of 32 bits, one at a time, and each
+one's address depends on the previous edge's index. There is almost nothing to overlap.
+ alt="Two schematics of one thread's memory traffic. Before, one feature per thread: a 32-bit load per edge, one edge in flight. After, four features per thread in a single 128-bit load, two edges unrolled and in flight.">
-
What one thread reads in the neighbor kernel, before and after. Drawing, not a measurement: the load widths and the edges in flight are read from the two kernels' code.
+
What one thread reads before and after. Drawing, not a measurement — the load widths and the edges in flight come from the two kernels' code. Four features per thread makes one 128-bit load; unrolling the edge walk by two keeps two of those loads in flight at once.
Download SVGPNG
@@ -27,33 +47,57 @@ thread read one float per edge — a 32-bit load, one edge at a time — so a ro
-Four features per thread makes that one 128-bit load, and unrolling the edge walk by two puts two
-gathers in flight. Measured 9.45 → 7.06 µs across 100 launches at 4096×64×65536. Triton's kernel,
-which loads 32 edges at once into a tile, is still ahead at 5.97 µs; the gathered rows are about
-17 MB, which is L2 traffic on this part, so more edges in flight is what is left.
-
-Layer normalization at 4096-wide rows needed no rewrite. The kernel for these rows splits each row
-across two warps, and the two partial sums meet in 64 bytes of shared memory behind one barrier.
-It was gated at width 2048, so 4096 fell through to a kernel that walks the row in a single warp.
-The gate was stale: the only structural requirement is that each warp's span be a whole number of
-32-lane steps, and at 4096 each half is 2048 elements, which is. Removing the gate measured
-255.41 → 186.80 µs, against Triton's 181.59.
-
-Bias + GELU took the same treatment and lost 2.6×. The shape is 4096×768: 3.1M elements, or 12,288
-blocks of 256 threads with one element each, which already fills the device. Four elements per
-thread cuts that to 786,432 threads and gives each thread four serial `tanh` evaluations. Measured
-27.00 µs against 10.29 µs, and reverted. Vector width pays when a thread is short of work to
-issue; it costs when the grid is already saturating the machine.
-
-| Operation | Before | After | Reference (other session) |
+Four features per thread turns each of those loads into a 128-bit quad, and unrolling the walk by
+two gives the memory system two independent requests instead of a dependency chain. Measured across
+100 launches: $9.45 \to 7.06\ \mu s$. Triton's kernel, which loads $32$ edges at once into a tile
+and reduces over them, is still ahead at $5.97\ \mu s$: the gathered rows are about
+$65536 \times 256\ \mathrm{B} = 16.8\ \mathrm{MB}$, which is L2 traffic on this part, so the only
+thing left to win is more requests in flight.
+
+## Argument two: sometimes nothing needs to change
+
+Layer normalization computes, per row:
+
+$$y_f \;=\; \gamma_f \, \frac{x_f - \mu}{\sqrt{\sigma^2 + \epsilon}} \;+\; \beta_f$$
+
+and the kernel for it holds the row in registers, so the split of the row across warps is fixed by
+size. At 4096-wide rows the kernel that handles them splits each row across two warps and the two
+partial sums meet in 64 bytes of shared memory behind one barrier; it was gated at width 2048, so
+$4096$ fell through to a kernel that walks the whole row in a single warp. The gate was stale. The
+only structural requirement is that each warp's span be a whole number of 32-lane steps, and at
+$4096$ each half is $2048$ elements, which is. Removing the gate measured $255.41 \to
+186.80\ \mu s$, against Triton's $181.59\ \mu s$ — the largest single win of the three, from
+deleting one condition.
+
+## Argument three: the counterexample
+
+Bias + GELU has no reuse and no gathering. It applies
+
+$$\mathrm{gelu}(z) \;=\; \tfrac{1}{2} z \left(1 + \tanh\!\left(\sqrt{\tfrac{2}{\pi}}\left(z + 0.044715\,z^3\right)\right)\right)$$
+
+to $z = x + \text{bias}$, elementwise. At $4096 \times 768$ that is 3.1M elements; with one element
+per thread it is 12,288 blocks of 256 threads, which already fills the device, and every thread's
+work is independent. Four elements per thread cuts that to 786,432 threads and gives each thread
+four `tanh` evaluations in sequence. Measured $27.00\ \mu s$ against $10.29\ \mu s$ — the change was
+reverted, and the scalar kernel stayed.
+
+That is the whole lesson of this note: **the same edit is a $1.4\times$ win on two kernels and a
+$2.6\times$ loss on a third**, and what separates them is whether the kernel was short of
+instructions or short of threads.
+
+## Evidence
+
+| Operation | before | after | reference (other session) |
| --- | ---: | ---: | --- |
-| Neighbor aggregation 4096×64×65536 | 9.45 | **7.06** | Triton 5.97, PyTorch 120.93 |
-| LayerNorm 4096×4096 | 255.41 | **186.80** | Triton 181.59 |
-| Bias + GELU 4096×768 | 10.29 (kept) | 27.00 (reverted) | Triton 7.76, cuTile 7.97 |
+| Neighbor aggregation $4096\times64\times65536$ | 9.45 | **7.06** | Triton 5.97, PyTorch 120.93 |
+| LayerNorm $4096\times4096$ | 255.41 | **186.80** | Triton 181.59 |
+| Bias + GELU $4096\times768$ | 10.29 (kept) | 27.00 (reverted) | Triton 7.76, cuTile 7.97 |
*Microseconds of GPU kernel time, mean of 100 launches after 25 warm-up launches. Each before and
-after pair was measured in one session on an idle device; the references are quoted from mage-004
-and mage-006 and are drawn hatched in the figure below because they are not pairings.*
+after pair was measured in one session on an idle device; the references come from mage-004 and
+mage-006 and are hatched in the figure below because they are not pairings. Worst full-output error
+against PyTorch: $4.8\times10^{-7}$ (LayerNorm), $3.6\times10^{-7}$ (neighbor), $5.3\times10^{-6}$
+(GELU).*
@@ -63,7 +107,7 @@ and mage-006 and are drawn hatched in the figure below because they are not pair
alt="Three rows of horizontal bars of GPU kernel time in microseconds, each row before and after one change. Neighbor aggregation falls from 9.45 to 7.06 with Triton at 5.97 for reference. LayerNorm at 4096 by 4096 falls from 255.41 to 186.80 with Triton at 181.59. Bias plus GELU keeps its scalar kernel at 10.29 while the feature-quad variant that measured 27.00 was reverted, with Triton at 7.76 and cuTile Rust at 7.97.">
-
The same three changes as bars, with the other-session references hatched. Lower is better.
+
The three changes as bars, with the other-session references hatched. Lower is better.
Download SVGPNG
@@ -85,15 +129,15 @@ and mage-006 and are drawn hatched in the figure below because they are not pair
-## Limits of these numbers
+## What the numbers do not establish
-- No hardware counters are available, so "L2 traffic" is arithmetic on the bytes moved, not a
- measurement.
+- No hardware counters are available, so "$16.8$ MB of L2 traffic" is arithmetic on the bytes moved,
+ not a measurement.
- One capture per point: a kernel time carries no interval of its own.
-- The references come from other sessions (`scripts/evolve_capture.py --all` runs every
- implementation in one session, and has not been run on these shapes).
-- The LayerNorm change is only validated where the harness points it: 4096×4096, 4096×512,
- 3072×1024, 2048×2048. Widths above 4096 still fall through.
+- The references come from other sessions; `scripts/evolve_capture.py --all` runs every
+ implementation in one pass and has not been run on these shapes.
+- The LayerNorm gate is only validated where the harness points it: $4096\times4096$,
+ $4096\times512$, $3072\times1024$, $2048\times2048$. Widths above 4096 still fall through.
- GELU's rejection is one pair, not a sweep: two features per thread, and quads with other block
sizes, were not tried.
diff --git a/scripts/plot-matmul-layouts.py b/scripts/plot-matmul-layouts.py
new file mode 100644
index 0000000..95ba6b3
--- /dev/null
+++ b/scripts/plot-matmul-layouts.py
@@ -0,0 +1,134 @@
+"""Two libraries, the same matrices, different memory.
+
+Draws how cuTile Rust and cuda-oxide place the same C = A B problem: one tile per
+program on the left, one register tile per thread on the right, with the shared
+memory layouts and load widths annotated. Schematic, taken from both kernels'
+code, not a measurement.
+
+Run: uv run --script scripts/plot-matmul-layouts.py
+"""
+from pathlib import Path
+
+import matplotlib
+
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt # noqa: E402
+from matplotlib.patches import FancyArrowPatch, Rectangle # noqa: E402
+
+ROOT = Path(__file__).resolve().parents[1]
+OUT = ROOT / "docs/assets/figures/mage-004"
+
+BG, INK, MUTED, RULE = "#101217", "#edf0f5", "#a0a9b9", "#303641"
+CUTILE, OXIDE, SHARED = "#c9b2ff", "#91dbba", "#93caff"
+
+
+def cell_grid(ax, x, y, cols, rows, size, face, edge=RULE, lw=.6, z=3):
+ for c in range(cols):
+ for r in range(rows):
+ ax.add_patch(Rectangle((x + c * size, y - r * size), size - .6, size - .6,
+ facecolor=face, edgecolor=edge, linewidth=lw, zorder=z))
+
+
+def label(ax, x, y, text, color=MUTED, size=8.2, ha="left", weight="normal"):
+ ax.text(x, y, text, color=color, fontsize=size, ha=ha, va="top", weight=weight)
+
+
+def panel_title(ax, title, subtitle, color):
+ ax.text(1, 98, title, color=color, fontsize=12, weight="bold", va="top")
+ ax.text(1, 90.5, subtitle, color=MUTED, fontsize=8.6, va="top")
+
+
+def cutile_panel(ax):
+ ax.set_xlim(0, 120)
+ ax.set_ylim(0, 100)
+ ax.axis("off")
+ panel_title(ax, "cuTile Rust: one tile per program", "the compiler owns the threads and the layout", CUTILE)
+
+ # the output, partitioned into 32 x 128 tiles
+ ax.add_patch(Rectangle((8, 22), 62, 62, facecolor="none", edgecolor=RULE, linewidth=1, zorder=2))
+ for c in range(4):
+ for r in range(4):
+ lit = (c, r) == (1, 1)
+ ax.add_patch(Rectangle((8 + c * 15.5, 22 + r * 15.5), 15.5 - .7, 15.5 - .7,
+ facecolor=CUTILE if lit else "none", edgecolor=RULE,
+ linewidth=.8, alpha=.95 if lit else 1, zorder=3))
+ label(ax, 8, 88, r"$C = A B$", color=INK, size=11)
+ label(ax, 8, 84, r"$C: 1024 \times 1024$ in $32 \times 128$ tiles", size=8)
+ label(ax, 8, 19, r"one program per lit tile", size=8, color=CUTILE)
+
+ # what it reads per K step
+ ax.add_patch(FancyArrowPatch((74, 59), (86, 59), arrowstyle="-|>", mutation_scale=10,
+ color=CUTILE, linewidth=1.3, zorder=4))
+ label(ax, 86, 66, r"$A$ tile $32\times32$", size=8, color=INK)
+ label(ax, 86, 60, r"$B$ tile $32\times32$", size=8, color=INK)
+ label(ax, 86, 54, r"$32$ steps of $K$", size=8)
+ label(ax, 86, 47, "widened to 128-bit loads", size=8, color=MUTED)
+ label(ax, 86, 41, "tile shape is fixed at", size=8, color=MUTED)
+ label(ax, 86, 35.5, "compile time, so it is part", size=8, color=MUTED)
+ label(ax, 86, 30, "of what gets compiled", size=8, color=MUTED)
+
+
+def oxide_panel(ax):
+ ax.set_xlim(0, 120)
+ ax.set_ylim(0, 100)
+ ax.axis("off")
+ panel_title(ax, "cuda-oxide: one register tile per thread", "the kernel author owns both", OXIDE)
+
+ # the block tile and one thread's 4 x 4 registers
+ ax.add_patch(Rectangle((8, 22), 62, 62, facecolor="none", edgecolor=RULE, linewidth=1, zorder=2))
+ for c in range(16):
+ for r in range(16):
+ lit = 4 <= c < 8 and 4 <= r < 8
+ ax.add_patch(Rectangle((8 + c * 3.875, 22 + r * 3.875), 3.875 - .35, 3.875 - .35,
+ facecolor=OXIDE if lit else "none", edgecolor=RULE,
+ linewidth=.4, zorder=3))
+ label(ax, 8, 88, r"$C = A B$", color=INK, size=11)
+ label(ax, 8, 84, r"$C: 64 \times 64$ per block, $16 \times 16$ threads", size=8)
+ label(ax, 8, 19, "lit = one thread's 4 x 4 registers", size=8, color=OXIDE)
+
+ # shared memory: A transposed so a thread's four rows are contiguous
+ ax.add_patch(FancyArrowPatch((74, 59), (86, 59), arrowstyle="-|>", mutation_scale=10,
+ color=OXIDE, linewidth=1.3, zorder=4))
+ label(ax, 86, 66, r"shared $A^T$: $64 \times 68$", size=8, color=INK)
+ label(ax, 86, 60, r"row stride 68, not 64", size=8, color=SHARED)
+ label(ax, 86, 54, "one 128-bit read per", size=8, color=MUTED)
+ label(ax, 86, 48, "row, four rows at once", size=8, color=MUTED)
+ label(ax, 86, 41, r"shared $B$: $64 \times 64$", size=8, color=INK)
+ label(ax, 86, 35.5, r"$0.125$ shared reads per", size=8, color=SHARED)
+ label(ax, 86, 29.5, r"multiply-add", size=8, color=SHARED)
+ label(ax, 86, 22, "stride 68 avoids bank conflicts", size=8, color=MUTED)
+
+
+def main():
+ fig, axes = plt.subplots(1, 2, figsize=(10.4, 3.9))
+ fig.subplots_adjust(left=.01, right=.99, top=.80, bottom=.02, wspace=.05)
+ fig.text(.01, .965, "The same product, two ways of placing it in memory", color=INK,
+ fontsize=13, weight="bold", va="top")
+ fig.text(.01, .915, r"$C = A B$ at $1024^3$, FP32. The compiler's kernel and the hand-written one "
+ "move different amounts of data to reach the same answer.", color=MUTED,
+ fontsize=9, va="top")
+ cutile_panel(axes[0])
+ oxide_panel(axes[1])
+ for ax in axes:
+ ax.set_facecolor(BG)
+ fig.patch.set_facecolor(BG)
+ OUT.mkdir(parents=True, exist_ok=True)
+ metadata = {"Date": None, "Description":
+ "Schematic of the matmul tile layouts in cuTile Rust and cuda-oxide; not a "
+ "measurement. See docs/experiments/mage-004.md."}
+ for name in ("matmul-layouts", "matmul-layouts-mobile"):
+ if name.endswith("mobile"):
+ fig.set_size_inches(4.2, 6.4)
+ axes[0].set_position([.02, .52, .96, .36])
+ axes[1].set_position([.02, .06, .96, .36])
+ fig.savefig(OUT / f"{name}.svg", metadata=metadata)
+ svg = OUT / f"{name}.svg"
+ svg.write_bytes(b"\n".join(line.rstrip() for line in svg.read_bytes().splitlines()) + b"\n")
+ if not name.endswith("mobile"):
+ fig.savefig(OUT / f"{name}.png", dpi=200, metadata=metadata)
+ plt.close(fig)
+ print("wrote", OUT / "matmul-layouts.svg")
+
+
+if __name__ == "__main__":
+ main()
+ 
+ alt="Two schematics of one thread's memory traffic. Before, one feature per thread: a 32-bit load per edge, one edge in flight. After, four features per thread in a single 128-bit load, two edges unrolled and in flight.">
