GPU stencil
research
project
Synchronous vs block-async
Jacobi on one RTX 5050
512² → 18432² grids
Separating synchronization
from temporal locality.
1.3×
Speedup from exact temporal
blocking once the grid leaves L2
About the project
Block-asynchronous iterative methods relax global synchronization, trading cheaper iterations for extra iterations caused by stale data. I examined this trade-off for a five-point Jacobi solver of an implicit diffusion step on a single consumer GPU, an NVIDIA GeForce RTX 5050 Laptop GPU, on grids from 512² to 18432².
The block-asynchronous kernel does two things at once: it tolerates stale data, and it reads each grid point from memory only once per D sweeps. The second effect is temporal locality, which an exact, fully synchronous method can also exploit. This project separates the two.
Learn more→What I
did
The problem01 / 07
Jacobi update with uⁿ held fixed
The transient equation ∂u/∂t = α∇²u + f is discretized on an N×N unit square. Each backward-Euler step solves a linear system; Jacobi sweeps it to convergence. The top boundary is held at 100, the others at 0, and the source is either a point at the centre or a Gaussian with σ = N/16. Each run solves ten consecutive time steps in float32.
The iteration matrix has spectral radius below 4r/(1+4r), so the number of sweeps grows with r but not with N: about 35, 226 and 1417 sweeps per step for r = 1, 10 and 100.
A step converges when the relative true residual ρ = ‖b − Au‖∞ / ‖b‖∞ drops below 10⁻⁵, checked every M = 8 sweeps. Because A is strictly diagonally dominant, the residual directly bounds the solution error. It is computed in a difference form that avoids float32 cancellation, which the naive form underestimated by about 2.5× at r = 100.
Can relaxing global synchronization reduce time to convergence despite the extra iterations it causes, and under which conditions?
Can the best synchronization strategy be predicted from features known before execution?
A mistake worth sharing. The first solver used the current iterate in place of uⁿ. It still converged, but to the steady-state Poisson problem, needing 3.7×10⁵ sweeps at N = 512 and 1.25×10⁶ at N = 1024, a count that grows as O(N²). Fixing it changed iteration counts by three orders of magnitude and made r a meaningful parameter.
Synchronous
Each sweep is its own kernel launch. No block may start sweep k+1 until every block has finished sweep k.
Block-asynchronous
A block sweeps its tile D times in shared memory, using frozen halo values from its neighbours, before anyone exchanges data.
S = Tsync / Tx = 1 / (c · s) — c: cost per sweep · s: sweep ratio
Try it · staleness demo
The same equations, re-implemented in JavaScript on a 96² grid with 32×8 tiles. Pick r and D; the page runs ten time steps and compares the last one for each strategy.
CPU float64, illustrative only. The ratio s matches the GPU closely because it is set by the numerics, not by N; whether fewer launches pay for extra sweeps is set by the hardware, which the GPU results below answer.
Method · 02 / 07
All kernels share one arithmetic definition and are validated bit for bit against a NumPy reference, down to identical sweep counts.
One thread per point in 32×8 blocks, ping-pong buffers. A reduction kernel computes both residual norms every M sweeps with warp shuffles and copies them to the host, which decides whether to continue.
Strengthened later with a fused right-hand side (uⁿ + Δt·f computed once per step), cutting traffic from 16 to 12 bytes per point per sweep.
After Anzt et al.: each block loads its tile plus a one-cell halo into shared memory and sweeps it D ∈ {1, 2, 4, 8} times with __syncthreads() between sweeps. Values read across a tile edge can be up to D−1 sweeps old.
With D = 1 it reproduces the synchronous iterates exactly.
D is just a kernel argument, so it can change between launches for free. G-asc and G-desc step D up or down as ρ crosses 10⁻², 10⁻³, 10⁻⁴ (latched, to stop oscillation).
G-fb is a hill-climbing controller: start at D = 1, measure −ln(ρₖ/ρₖ₋₁) per ms, double D while the rate holds, halve it on a >10% drop. A settling window skips the residual bump right after each change.
The ghost-zone technique of Meng & Skadron. Each block loads a halo D cells wide; in sweep s it updates only points at least s+1 cells from the edge. After D sweeps the tile holds exactly what D synchronous sweeps would produce.
Same data reuse as block-async, zero staleness. It pays with redundant work: a Tx×Ty tile loads (Tx+2D)(Ty+2D) points.
Round-to-nearest intrinsics (__fadd_rn, __fmul_rn, __fdiv_rn) fix the operation order and prevent fused multiply-adds, confirmed in the generated PTX. The GPU performs exactly the same float32 operations as the NumPy reference, so any difference means a bug, not rounding.
Grids like N = 70 leave partial tiles at the edges to catch boundary bugs.
Setup03 / 07
| GPU | RTX 5050 Laptop · CC 12.0 · 20 SMs · 8 GB |
| L2 cache | 32 MiB |
| DRAM bandwidth | ≈ 384 GB/s nominal (128-bit bus) |
| Toolkit | CUDA 13.3 · driver 610.57.04 |
| Flags | -O3 -arch=sm_120 --fmad=false |
| Host | Fedora 44 · Python 3.14 · GCC 16.2 |
| Grid sizes N | 512, 1024 · extensions 1024–18432 |
| Coupling r | 1, 10, 100 |
| Time steps | 10 per run · ρ < 10⁻⁵ · M = 8 |
| Repeats | 3 (main) · 5 (N ≤ 2048) · 3 (N ≥ 4096) |
Main matrix. N ∈ {512, 1024}, both sources, three values of r and eight strategies with three repeats: 288 runs, all of which met the tolerance.
Extension experiments. Designed afterwards to explain the main matrix: a check-interval diagnostic, tile-shape and overlapped-check tests, the fused right-hand side and temporal blocking, and grid sizes up to 18432², the largest that fits in device memory. SM clock locked at 3090 MHz, 300 ms warm-up.
All ten time steps are timed with CUDA events, including convergence checks. All reported values are medians over repeats.
No hardware counters. Nsight Compute failed with a permission error on this driver setup, so cache residency is inferred from effective bandwidth (bytes per sweep × sweeps ÷ time) against the nominal DRAM peak.
Results · in cache04 / 07
Configurations where plain synchronous Jacobi is fastest. The two exceptions (N = 1024, r = 100) are within timing noise.
Extra sweeps D8 needs because of stale halos. D4 needs 1.35–1.59×; D1 exactly as many.
Best asynchronous config relative to best synchronous config, after trying three tile shapes and overlapped checks.
Table IV. Median of 3 runs, 32×8 tiles. Blue cells are faster than sync, red cells slower. Hover a cell for details.
Kernel launches are queued asynchronously, and global sweeps are served largely from L2: the fused sync kernel sustains 580–640 GB/s at N = 1024, well above the 384 GB/s DRAM peak, consistent with a 12–16 MiB working set in the 32 MiB L2.
The block-async kernel also carries a fixed cost per launch: at N = 512, a D1 launch took ≈ 9.4 µs and a D8 launch ≈ 19.4 µs against ≈ 4.4 µs for a sync sweep. Each extra in-tile sweep costs ≈ 1.4 µs, but loading and storing the tile costs about two synchronous sweeps.
The synchronization that costs time, the host round trip that decides convergence, is something inter-block asynchrony can't remove.
Table V. Gaussian source, 5 repeats, 300 ms warm-up. B = blocking checks, O = overlapped checks.
Results · beyond L205 / 07
Once the grid outgrows L2, the synchronous kernel becomes bandwidth-bound. Kernels that sweep several times per pass through device memory now have a real advantage, and both temporal blocking and block-asynchronous relaxation overtake the synchronous kernel.
The speedups did not keep growing with N. From 2048² to 18432² the working set grows 81×, yet speedups stay inside a flat band. Beyond L2 there is no further cache level to fall out of.
The two factors of S = 1/(c·s) explain why: s is fixed by the numerics (D8@32×16 needs exactly 1.738× the sync sweeps at r = 10 for every N from 1024 to 16384), and c is fixed by the memory regime (0.39–0.42 for D8, 0.71–0.76 for T4).
Exact temporal blocking over the fused sync baseline, out of cache, with no staleness at all.
Block-asynchronous over the same baseline at r = 10 (1.44–1.48× at r = 100, 1.07–1.44× at r = 1).
What relaxing synchronization adds over temporal blocking at r ≥ 10. At r = 1 it is 16–19% slower.
Fig. 6 / Table VIII at r = 10. A strategy wins exactly when c·s < 1. Block-async: cheap sweeps, but more of them. Temporal blocking: no extra sweeps, but nearly the full cost of a sync sweep because of redundant halo work.
How optimistic was the first result?
The first out-of-cache measurement looked like a clear win for asynchrony. It didn't survive a fair comparison, and how it shrank is one of the main findings.
Takeaway. When an asynchronous kernel also changes the memory access pattern, the comparison needs a synchronous baseline that changes it the same way. Fusing the right-hand side alone removed ≈ 15% of the apparent advantage; most of what remained came from temporal locality, not stale data.
Table VI. Tunfused/Tfused, Gaussian source. Sync gets 1.17–1.24× faster out of cache; D8 and T8 barely move, because they already read those arrays once per launch.
Adaptive & learned selection06 / 07
G-fb came within 10% of the best static strategy per configuration without being told N or r, and beat every static strategy in one configuration. At r = 100 it changed D between 186 and 212 times per run, each change reacting to the residual disturbance of the previous one.
Out of cache it was 1.17–1.54× faster than the unfused sync kernel at N = 2048 and 4096, within 10–15% of the best static choice. In cache it lost (0.69–0.84×) because its search pushes towards larger D, which doesn't pay there. A version that can fall back to synchronous behaviour would recover much of that.
Depth-2 decision tree, 16 configs
Split values are midpoints between measured points. The data establish only that the switch lies between N = 1024 (ratio 0.5) and 2048 (ratio 2.0), an interval containing the physical L2 capacity, and between r = 1 and 10.
X = pd.DataFrame({"ws_over_L2": 16 * N**2 / 2**25, "r": r}) y = fam[["sync", "temporal", "staleness"]].idxmin(axis=1) for held_out in set(groups): # config or grid size test = groups == held_out tree = DecisionTreeClassifier(max_depth=2, random_state=52) tree.fit(X[~test], y[~test]) pred[test] = tree.predict(X[test])
On the main matrix, a tree on (N, source, r) picked the fastest static strategy every time, but only two configurations had a non-sync label, so a second selector was trained on the 16-config extension dataset: classes sync (3), temporal blocking (5) and block-async (8).
Conclusion07 / 07
While the working set fits in L2, plain synchronous Jacobi is fastest. Barriers between blocks are cheap; the costly synchronization is the host round trip for convergence checks, which inter-block asynchrony can't touch.
Sweeping several times per pass through memory gives a step change of ≈ 1.3–1.5× that stays flat up to the largest grid that fits in memory. About 1.3× is available to an exact method with no staleness.
Relaxing synchronization adds 10–13% over temporal blocking when time steps need many iterations (r ≥ 10), because a one-cell halo is cheaper than redundant halo work. When steps are short (r = 1) it is 16–19% slower.
Working set relative to L2, and r. A depth-2 tree recovers this rule and cuts the best fixed strategy's mean slowdown by almost two-thirds (13.6% → 4.9%), but cannot extrapolate across the cache boundary.