100 DAYS OF GPU Course preview
Day 34Module 4
in-technical-review

CUDA stencils and the heat equation

Here is the loop almost everyone writes the first time they put a simulation on a GPU.

for (int step = 0; step < 1000; ++step) {
    heatStep<<<grid, block>>>(d_cur, d_next, n, r);
    cudaMemcpy(h_grid, d_next, bytes, cudaMemcpyDeviceToHost);
    std::swap(d_cur, d_next);
}

It is correct. Every step is right and the picture it draws is the picture. It also spends most of its time on line three.

A 2048 by 2048 grid of floats is 16 MiB, so a thousand steps copy 15.6 GiB to the host, one blocking copy at a time, to fetch a grid that was already on the device. Day 9 measured that cost: a vector add whose kernel took 0.786 ms took 45.822 ms once its copies were counted.

This page builds the stencil that loop is calling, gives it a shared-memory tile with a halo, and then measures the part that sets the run time. It is not the kernel.

What a stencil asks memory for

A stencil computes each cell's new value from its own old value and its neighbours' old values. The two-dimensional heat equation, solved explicitly, is the smallest interesting one:

u'[i][j] = u[i][j] + r * (u[i-1][j] + u[i+1][j] + u[i][j-1] +
                          u[i][j+1] - 4 * u[i][j])

Five cells in, one cell out, every cell, every step. That is the memory cost. Each updated cell costs five reads from global memory, and four of those five are cells some neighbouring thread also wants.

The reads themselves are fine: the 32 lanes of a warp hold 32 consecutive columns of one row, so each of the five loads is a contiguous 128-byte span and coalesces into four sectors, exactly as day 11 measured. The access pattern is fine, but the read count is high.

A stencil is the natural shape for shared memory because its reuse is small, local and known before the kernel runs. A block that owns a 32 by 8 patch of output needs exactly one extra cell on each side to compute it, so it loads a 34 by 10 patch, waits at one __syncthreads(), and then every thread reads its five neighbours out of shared memory. That ring of extra cells is the halo, also called the apron, and it is the same object day 20 loaded for a 5 by 5 filter.

The arithmetic is the design. 34 times 10 is 340 loads to produce 32 times 8 = 256 updated cells, which is 1.33 global reads per cell against five. The 84 extra loads are the halo, and a neighbouring block also loads each one.

A tiled matmul reuses a whole row and a whole column; a stencil reuses a ring one cell wide. That is why the tile only has to be two cells bigger in each direction, and also why the win is smaller than a matmul's.

A 32 by 8 output tile uses 256 threads to load a 34 by 10 shared tile, totaling 340 loads or 1.328125 per output.

The thousand steps are the program

You cannot update a grid in place. A thread that writes u[i][j] before its neighbour reads it gives that neighbour a value from the wrong step, and the answer then depends on the block schedule.

So a stencil runs on two buffers and swaps them after every step. CUDA code often calls this ping-pong buffering. The swap changes two host pointers; it does not move data.

Both buffers need the boundary. These kernels never write the boundary ring, because a Dirichlet wall holds its temperature forever, so a run that filled one device buffer and left the other as cudaMalloc returned it reads real wall values on even steps and whatever was in that allocation on odd ones.

cudaMemcpy(d_a, h_init, bytes, cudaMemcpyHostToDevice);
cudaMemcpy(d_b, h_init, bytes, cudaMemcpyHostToDevice);  // both, every time

Now the loop around all that. The CPU habit is that a step is a function which takes data and returns data, so you bring the result back to look at it.

On a GPU, keep the loop on the device and let the host choose when to copy a result. Kernel launches are cheap: day 9 measured 0.016 ms of host time to queue one, so a thousand of that overhead is a rounding error. The copies are not, because they cross PCIe rather than the card's own memory, and a synchronous one also stops the host until it finishes.

Part 2 measures the three copy policies.

When the timestep is the bug

Write the update as a weighted sum instead of a correction. The stability limit is then clear:

u'[i][j] = (1 - 4r) * u[i][j] + r * (four neighbours)

r is alpha * dt / h^2, the diffusion constant times the timestep over the squared cell spacing. The five weights sum to one. While 1 - 4r is not negative, every weight is between zero and one.

The new value is then an average of five old values and cannot be larger than the largest of them. Nothing in the grid can grow. This discrete maximum principle gives the CFL condition r <= 1/4 for this two-dimensional scheme.

Past a quarter the centre weight goes negative, the sum stops being an average, and the scheme amplifies instead of smoothing. The pattern it amplifies fastest is the checkerboard, the one that alternates sign between neighbouring cells, and it is multiplied by |1 - 8r| on every step. At r = 0.26 that is 1.08, which sounds harmless until you compound it: a factor of 2,200 over a hundred steps and more than a trillion over four hundred.

What a learner sees is a picture that goes chequered, then to numbers with no physical meaning, growing without bound toward inf and nan if the run is long enough.

The first guess is often a kernel bug, so the next check focuses on indices and barriers. But the kernel can be correct. Halve the timestep and the same binary produces a stable simulation.

One dimension allows r <= 1/2, two allows a quarter, three allows a sixth: the bound is 1 / (2d), because each extra dimension puts two more neighbours' weight on the same centre.

Measuring it without measuring the copies

Full program in code/day34-stencils/heat.cu. The harness checks its numbers four ways.

The baseline is measured here, not borrowed. copyGrid reads every cell and writes it back with no neighbours and no arithmetic, from the same grid in the same process, so no correct stencil can beat it. Taking a copy bandwidth from day 11 instead would fold that day's block shape, buffer size and clock state into this answer.

Only one thing changes between the rows of part 2. Same kernel, same thousand steps, same arithmetic, and the only variable is how often the grid is copied to the host. The program then requires the three final grids to be bit-identical, so "the copy policy changed the answer" is a failure rather than a footnote.

Correctness is checked before anything is timed. Both kernels run four steps and are compared cell by cell against a double-precision CPU reference, and every cell of the boundary ring has to hold exactly the value it started with.

Events, and a warm-up per kernel. Day 9 covers why a host clock around a launch measures the launch. Part 2 is timed once rather than ten times, because a thousand steps is not a microbenchmark.

The global-memory version is the stencil written the obvious way:

__global__ void heatStepGlobal(const float* __restrict__ in,
                               float* __restrict__ out, int n, float r) {
    const int col = static_cast<int>(blockIdx.x * blockDim.x + threadIdx.x) + 1;
    const int row = static_cast<int>(blockIdx.y * blockDim.y + threadIdx.y) + 1;
    if (col < n - 1 && row < n - 1) {
        const size_t stride = static_cast<size_t>(n);
        const size_t c = static_cast<size_t>(row) * stride + col;
        const float centre = in[c];
        out[c] = centre + r * (in[c - 1] + in[c + 1] + in[c - stride] +
                               in[c + stride] - 4.0f * centre);
    }
}

The tiled version differs only in where the five reads come from. Inspect the fill loop: 256 threads have 340 apron cells to load, so it is strided by the block size rather than written as one load per thread, and 84 threads take a second cell.

__global__ void heatStepTiled(const float* __restrict__ in,
                              float* __restrict__ out, int n, float r) {
    __shared__ float tile[kSharedY][kSharedX];

    const int tileCol = static_cast<int>(blockIdx.x) * kTileX + 1;
    const int tileRow = static_cast<int>(blockIdx.y) * kTileY + 1;

    const int threads = static_cast<int>(blockDim.x * blockDim.y);
    const int flat = static_cast<int>(threadIdx.y * blockDim.x + threadIdx.x);
    for (int s = flat; s < kSharedCells; s += threads) {
        const int sy = s / kSharedX;
        const int sx = s - sy * kSharedX;
        const int gy = tileRow - kHalo + sy;
        const int gx = tileCol - kHalo + sx;
        const bool inside = (gx >= 0 && gx < n && gy >= 0 && gy < n);
        tile[sy][sx] =
            inside ? in[static_cast<size_t>(gy) * static_cast<size_t>(n) + gx]
                   : kCold;
    }
    __syncthreads();

    const int sx = static_cast<int>(threadIdx.x) + kHalo;
    const int sy = static_cast<int>(threadIdx.y) + kHalo;
    const int col = tileCol + static_cast<int>(threadIdx.x);
    const int row = tileRow + static_cast<int>(threadIdx.y);
    if (col < n - 1 && row < n - 1) {
        const float centre = tile[sy][sx];
        out[static_cast<size_t>(row) * static_cast<size_t>(n) + col] =
            centre + r * (tile[sy][sx - 1] + tile[sy][sx + 1] +
                          tile[sy - 1][sx] + tile[sy + 1][sx] - 4.0f * centre);
    }
}

There is no padding column on that tile. Day 15 added one to a transpose because a warp there walked a column of the tile; here all five shared reads walk a row with consecutive lanes on consecutive words, so the 32 lanes already land in 32 different banks and a 35th column would buy nothing but shared memory. Check that before adding the + 1.

Bank conflicts are the day that teaches it.

Note. The tiled kernel may not win. Five reads per cell is a small amount of reuse, and every row is read three times within a few microseconds by threads on the same SM, so L1 and L2 have a good chance of serving most of it already. Shared memory would then be doing by hand what the cache was doing for free. That is what part 1 is for, and on this card it settled the question: the tiled kernel lost, 0.283 ms against 0.173. The Results section takes it apart.

Results

Measured. Tesla T4, driver 595.84, CUDA 12.6 (V12.6.85), built with nvcc -O3 -arch=sm_75. Captured 2026-09-01 on the project's verification node; full transcript in the page's evidence file.

GPU: Tesla T4 (compute capability 7.5)
Shared memory per block: 49152 bytes. Max threads per block: 1024
Grid: 2048 x 2048 cells, 16.0 MiB per buffer, two buffers on the device
Tile: 32 x 8 output from a 34 x 10 apron, 1360 bytes of shared memory
Interior cells updated per step: 4186116

Global reads per updated cell, from the constants
  heatStepGlobal        5.000
  heatStepTiled         1.328   (340 apron cells for 256)

Part 1: one step, mean of 10 runs after 3 warm-ups
kernel                   time (ms)      GB/s    % of copy
---------------------   ----------  --------  -----------
copyGrid, no stencil         0.138     244.0        100.0
heatStepGlobal               0.173     194.2         79.6
heatStepTiled                0.283     118.6         48.6

Part 2: 1000 steps of heatStepTiled. Only the copy policy changes.
policy                        total (ms) per step (ms) to host (GiB)    vs row 1
---------------------------   ----------  ------------  ------------  ----------
never copied back              200.894         0.201           0.0         1.0
copied back every 100 steps     230.566         0.231           0.2         1.1
copied back every step        4117.288         4.117          15.6        20.5

All three rows above produced a bit-identical final grid. The kernel,
the step count and the arithmetic were the same; only the copies moved.

Part 3: 400 steps at three timesteps
       r     1 - 4r     |1 - 8r|     max |u| at the end    finite
 -------  ---------  -----------  ---------------------  --------
   0.240      0.040        0.920                      1       yes
   0.250      0.000        1.000                      1       yes
   0.260     -0.040        1.080            4.83157e+09       yes

The grid started with max |u| = 1.000. A row whose centre weight 1 - 4r is
negative is no longer averaging its neighbours, and |1 - 8r| is what the
checkerboard pattern is multiplied by on every one of the 400 steps.

Tiling made it slower, and that is the result worth keeping. The global version reads 5.000 values per updated cell and runs in 0.173 ms. The tiled version reads 1.328 and runs in 0.283 ms.

Fewer global reads took more time.

The reason is that the reads the tile eliminates were not expensive. A five-point stencil re-reads its neighbours, and those neighbours were just read by an adjacent thread, so they are sitting in L1 and L2. The tile replaces a cheap cache hit with shared memory traffic, a barrier, and an apron that is 33 percent larger than the output it serves: 340 cells loaded for 256 updated.

This is the counterweight to day 13, where tiling nearly doubled a transpose. There, the reads being eliminated were uncoalesced DRAM traffic.

Here they are cache hits. Shared memory is worth it when it replaces slow reads, not merely repeated ones, and the only way to know which you have is to measure both.

The copies cost 20.5x, the page's main result. A thousand steps with the grid left on the device took 200.894 ms. The same thousand steps copied back after every one took 4117.288 ms and moved 15.6 GiB to the host.

All three policies produced a bit-identical final grid.

The row in between, copied back every 100 steps, reads 230.566 ms, 1.1x. An earlier capture of this table had the middle row beating the first, which made no sense for a row doing strictly more work; a discarded warm-up pass before part 2 removed it, and the rows now sit in the only order the work supports.

The blow-up is smaller than the compounding suggests. The model says the checkerboard component is multiplied by 1.080 per step at r = 0.260, which compounds to about 2e13 over 400 steps. The run ends at 4.83157e+09, still finite, because the initial grid's checkerboard component is far smaller than 1: the multiplier is exact, but the starting amplitude is not 1.

The sign of 1 - 4r correctly predicts growth or decay, as the table shows, and the exact magnitude at step 400 depends on where the pattern started.

Run it yourself

Use a CUDA GPU that supports the lesson's minimum compute capability. The README gives this build line:

nvcc -std=c++17 -O3 -arch=sm_75 -o heat heat.cu

An embed is out: two 16 MiB grids on the device, about three thousand timed launches and a thousand 16 MiB copies to the host will not finish in the 20 seconds Compiler Explorer allows.

Shrinking the grid until it fits would change the copy cost this page measures. If you have no GPU, read /setup/learn-cuda-without-a-gpu.

Exercise

Add a kernel that takes two timesteps per launch. Load a two-cell halo, run the stencil once inside shared memory, run it again, and write the tile out, so a thousand steps cost five hundred launches and five hundred passes over global memory. Then say what those two steps cost you on the way in.

Time: 35 to 50 minutes. Submit: heat.cu with the new kernel, your two new rows, and one sentence on when you would ship it.

Check: the harness runs your kernel for the same thousand steps and compares every cell against the shipped single-step run, printing the row, the column and both values for the first cell that disagrees. It counts your global reads per updated cell and holds you to 1.5 per timestep, which the naive five-read version fails and a correct two-step tile passes with room. It reports your time per timestep against heatStepTiled either way, because the ratio is the answer and the pass is only the gate.

Hint 1

The block already has everything it needs to compute its own tile once. What would it need to have loaded to compute that tile a second time without asking global memory again?

Hint 2

Work backwards. The second step needs a one-cell ring around the output tile, and that ring has to be correct after the first step, so the first step has to compute it, so the load has to cover a ring around the ring. Count the apron, then count how many cells each of the two steps updates.

Solution

A 36 by 12 apron. The first step updates the inner 34 by 10 into a second shared array, the second updates the 32 by 8 tile from that. Any cell whose global coordinates put it on the domain boundary keeps its value in both steps, which is the same rule the single-step kernel follows and the reason the edge tiles need no special case.

That is 432 loads for 256 cells and two timesteps: 1.69 per launch, which looks worse than 1.33, and 0.84 per cell-step, which is the useful ratio. You pay for it in redundant loads, in shared memory (two arrays now), and in a first step where 256 threads have 340 cells to update, so some of them go round twice before the barrier.

The general trade is this: a stencil tile can cover more time steps as well as more grid cells. That needs a larger halo and more repeated loads, but it reduces trips to memory. Whether that trade is worth taking is a measurement, not a rule, and on a kernel already running at the copy rate there is little left to gain.

Pitfalls

Your simulation goes chequered, then to huge numbers, then to nan. The timestep is past the CFL bound, so the centre weight 1 - 4r is negative and the checkerboard pattern grows by |1 - 8r| every step.

Halve dt and run the same binary. Halving the cell spacing to get a finer picture means quartering dt, which is why refining a grid gets expensive faster than it looks.

You updated the grid in place. One buffer, and a thread reads a neighbour that has already advanced a step.

It does not crash, and the picture still looks like heat, but the error changes with the block schedule. Use two buffers and swap them after each step.

You filled one device buffer and not the other. The kernels never write the boundary ring, so the wall temperature alternates between the set value and whatever cudaMalloc handed you. Copy the initial condition into both, and check the ring exactly rather than within a tolerance.

You guarded the barrier instead of the load. On a grid that is not a whole number of tiles, some threads own no output cell, and the tempting move is to return early or to put the bounds check around the __syncthreads().

It may not hang. The threads that stay then read a tile that the others did not finish filling, so the error changes with the schedule.

Guard the load and the store, never the barrier. Day 14 takes apart what a barrier does and does not promise.

You copied the grid back so you could watch it. That is part 2 of this page. Keep the loop on the device and snapshot every hundredth step, or move the copy to a second stream with pinned memory so it overlaps the next step.

Day 51 is where the copies start hiding.

You added a padding column out of habit. Padding costs shared memory and buys nothing when every warp access walks a row of the tile, which is the case for all five reads here. Work out which bank each lane lands in before you spend the memory.

Go deeper

Next

Day 35 builds a radix sort on the scan from day 32, which uses data in a different way. This topic returns twice: day 51 moves the measured copy onto a second stream so it overlaps the next step, and day 43 asks the same question about matmul, where the answer comes out the other way because the kernel, not the loop, is what is left to fix.