code/day34-stencils/heat.cuThis is the source used by the lesson and its recorded evidence. Compile commands and expected output live in the directory README.
// SPDX-License-Identifier: MIT
//
// Day 34: the 2D heat equation, a haloed stencil tile, and a thousand steps
// that never leave the device.
//
// Three kernels over the same 2048 x 2048 grid, and one of them is a floor:
//
// copyGrid one read and one write per cell, no neighbours. The
// floor: no stencil over this grid can beat it.
// heatStepGlobal five reads per updated cell, straight from global
// memory.
// heatStepTiled one block loads a 34 x 10 apron into shared memory,
// then every one of its 256 threads reads its five
// neighbours from there. 340 global reads for 256 cells.
//
// Then the part the day is actually about. The same kernel runs 1000 steps
// three times, and the only thing that changes is where the grid lives
// between steps: never copied back, copied back every hundredth step, or
// copied back every step, which is what a first version usually does so its
// author can watch the simulation. The arithmetic is identical in all three.
//
// Then stability. The explicit scheme this program uses is
//
// 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])
//
// with r = alpha * dt / h^2. Rewrite it as a weighted sum of the five cells
// and the centre weight is 1 - 4r. At r = 0.25 that weight is zero and every
// weight in the sum is still non-negative, so the new value is an average of
// old values and cannot leave their range. Past r = 0.25 the centre weight
// is negative, the sum stops being an average, and the alternating
// checkerboard pattern grows by |1 - 8r| every step. That is the CFL
// condition for this scheme in two dimensions, and this program measures
// both sides of it.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -o heat heat.cu
// Run: ./heat
//
// Verified 2026-08-30 on a Tesla T4 (compute capability 7.5), driver
// 595.84, CUDA 12.6 (V12.6.85). Transcript: evidence/run-2026-08-30.txt
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <vector>
#include <cuda_runtime.h>
// The one error macro. This file is standalone, the way a Compiler Explorer
// embed is, so it carries its own verbatim copy. `err_` carries a trailing
// underscore so it cannot collide with a variable at the call site, and the
// do/while makes the macro one statement so it survives a braceless `if`.
#define CUDA_CHECK(call) \
do { \
cudaError_t err_ = (call); \
if (err_ != cudaSuccess) { \
std::fprintf(stderr, "CUDA error %s:%d: %s: %s\n", __FILE__, \
__LINE__, #call, cudaGetErrorString(err_)); \
std::exit(EXIT_FAILURE); \
} \
} while (0)
// The grid is kN x kN cells including the boundary ring, so the interior
// that actually changes is (kN - 2) x (kN - 2). 2048 squared floats is 16
// MiB per buffer, and the program holds two of them on the device, which is
// what ping-pong costs.
constexpr int kN = 2048;
constexpr size_t kCells = static_cast<size_t>(kN) * static_cast<size_t>(kN);
constexpr size_t kInterior =
static_cast<size_t>(kN - 2) * static_cast<size_t>(kN - 2);
constexpr int kSteps = 1000;
constexpr int kSnapshotEvery = 100;
// The output tile one block owns, and the apron it has to load to compute it.
// 32 wide means one warp covers one tile row, so every load in the fill loop
// below is 32 consecutive floats.
constexpr int kTileX = 32;
constexpr int kTileY = 8;
constexpr int kHalo = 1;
constexpr int kSharedX = kTileX + 2 * kHalo;
constexpr int kSharedY = kTileY + 2 * kHalo;
constexpr int kSharedCells = kSharedX * kSharedY;
// The 1D launch used by copyGrid and by nothing else.
constexpr int kThreadsPerBlock = 256;
// The five-point stencil: the cell itself plus its four edge neighbours.
constexpr int kTaps = 5;
// r = alpha * dt / h^2, the only number that decides whether the simulation
// is a simulation or a pile of floating-point noise. Two values at or under
// the bound, one just past it.
constexpr float kStableR = 0.24f;
constexpr float kBoundaryR = 0.25f;
constexpr float kUnstableR = 0.26f;
constexpr int kNumR = 3;
constexpr int kStabilitySteps = 400;
constexpr int kCheckSteps = 4;
constexpr float kRelTolerance = 1e-4f;
constexpr int kWarmupRuns = 3;
constexpr int kTimedRuns = 10;
constexpr float kHot = 1.0f;
constexpr float kCold = 0.0f;
// The maximum principle is checked in float against a double reference, so it
// gets a little room for rounding and nothing else. A run that breaks the
// bound breaks it by orders of magnitude, not by an ulp.
constexpr float kMaxPrincipleSlack = 1e-5f;
static_assert(kTileX * kTileY % 32 == 0,
"the block is a whole number of warps");
static_assert(kTileX * kTileY <= 1024,
"1024 threads is the per-block maximum on every compute "
"capability this course targets");
static_assert(static_cast<size_t>(kSharedCells) * sizeof(float) <= 48u * 1024u,
"the apron must fit the 48 KiB of shared memory a Turing block "
"gets without an opt-in call to cudaFuncSetAttribute");
static_assert(kStableR <= 0.25f && kBoundaryR <= 0.25f,
"both of these are meant to be at or under the 2D CFL bound of "
"one quarter, where 1 - 4r is still non-negative");
static_assert(kUnstableR > 0.25f,
"this one is meant to be past the bound, or part 3 has nothing "
"to show");
static_assert(kSteps % 2 == 0,
"with an even step count the ping-pong ends with the newest "
"grid back in the buffer it started in, which is also what a "
"kernel doing two steps per launch needs");
// Fills the initial condition: a cold grid with the left wall held hot and a
// hot square in the middle.
//
// The square's side is odd and its edges are sharp on purpose. A sharp edge
// carries every spatial frequency the grid can represent, including the
// alternating one that decides stability, so part 3 does not have to inject
// a perturbation to have something to amplify.
static void fillInitial(float* out, int n) {
for (size_t i = 0; i < static_cast<size_t>(n) * static_cast<size_t>(n);
++i) {
out[i] = kCold;
}
for (int row = 0; row < n; ++row) {
out[static_cast<size_t>(row) * n] = kHot;
}
const int half = n / 8;
const int lo = n / 2 - half;
const int hi = lo + 2 * half + 1;
for (int row = lo; row < hi; ++row) {
for (int col = lo; col < hi; ++col) {
out[static_cast<size_t>(row) * n + col] = kHot;
}
}
}
// One step of the reference, host only, in double. Written to be obviously
// right rather than fast; it never allocates and the caller owns both
// buffers.
//
// It copies the boundary rather than skipping it, because the reference has
// one buffer pair like the GPU does and a boundary cell that is not written
// here would be whatever the other buffer held two steps ago. The kernels
// take the other route and never write the boundary at all, which works only
// because both device buffers start out holding it.
static void heatStepCpu(const double* in, double* out, int n, double r) {
const size_t stride = static_cast<size_t>(n);
for (int row = 0; row < n; ++row) {
for (int col = 0; col < n; ++col) {
const size_t c = static_cast<size_t>(row) * stride + col;
if (row == 0 || col == 0 || row == n - 1 || col == n - 1) {
out[c] = in[c];
} else {
out[c] = in[c] + r * (in[c - 1] + in[c + 1] + in[c - stride] +
in[c + stride] - 4.0 * in[c]);
}
}
}
}
// Returns the first index where got and want differ by more than the relative
// tolerance, or n if they agree everywhere.
static size_t firstMismatch(const float* got, const double* want, size_t n,
float relTolerance) {
for (size_t i = 0; i < n; ++i) {
const double scale = (want[i] == 0.0) ? 1.0 : std::fabs(want[i]);
if (std::fabs(static_cast<double>(got[i]) - want[i]) >
relTolerance * scale) {
return i;
}
}
return n;
}
static float maxAbs(const float* x, size_t n) {
float biggest = 0.0f;
for (size_t i = 0; i < n; ++i) {
const float v = std::fabs(x[i]);
if (v > biggest) {
biggest = v;
}
}
return biggest;
}
static bool allFinite(const float* x, size_t n) {
for (size_t i = 0; i < n; ++i) {
if (!std::isfinite(x[i])) {
return false;
}
}
return true;
}
// One read and one write per cell, no neighbours and no arithmetic.
//
// Memory: consecutive lanes take consecutive cells, so one warp's 32
// addresses cover 128 contiguous bytes and cost four 32-byte sectors.
//
// This is the floor row. A stencil over this grid has to read every cell and
// write every interior cell whatever else it does, so it cannot beat this,
// and measuring it here rather than borrowing day 11's number keeps this
// day's block shape, buffer size and clock state out of the comparison.
//
// Launch assumption: any 1D grid that covers n. The guard is the bounds
// check.
__global__ void copyGrid(const float* __restrict__ in, float* __restrict__ out,
size_t n) {
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
if (i < n) {
out[i] = in[i];
}
}
// One thread updates one interior cell, reading all five of its inputs from
// global memory.
//
// Memory: 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. The pattern is not the problem here, the count is: five
// reads per updated cell, four of which some neighbouring thread also wants.
//
// The kernel never writes the boundary ring. That is correct only because
// both device buffers are initialised with the boundary values, which is the
// bug this file's README warns about twice.
//
// Launch assumption: a 2D grid whose x times y covers the (n - 2) squared
// interior. Both guards are bounds checks.
// snippet: global-kernel
__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);
}
}
// end snippet: global-kernel
// The same step, with the block's neighbourhood loaded once into shared
// memory and read from there five times.
//
// Memory: the fill loop below is the only global traffic, and it walks the
// apron row by row, so one warp's 32 addresses are again 32 consecutive
// floats. 340 apron cells for 256 output cells is 1.33 global reads per
// updated cell against five, and the 84 extra loads are the halo, which is
// what a tile with a halo costs.
//
// No padding column. Day 15 added one to a transpose because a warp there
// walked a column of the tile; here every one of the five shared reads is a
// row walk with consecutive lanes on consecutive words, so the 32 lanes
// already land in 32 different banks and a 35th column would buy nothing.
//
// Launch assumption: blockDim must be exactly (kTileX, kTileY), because the
// output mapping below takes threadIdx.x and threadIdx.y as the position
// inside the tile. The fill loop is written against blockDim rather than
// against the constants so that it stays correct while somebody is changing
// the tile shape.
// snippet: tiled-kernel
__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);
}
}
// end snippet: tiled-kernel
// Times a launch with CUDA events and returns the mean milliseconds per run.
// A host-side clock around a launch measures the launch, not the kernel,
// because launches are asynchronous. Day 9 takes that apart.
//
// This is the one template and the one lambda allowed in module 1 to 3 code.
// Copy it verbatim.
template <typename LaunchFn>
static float timeKernel(LaunchFn launch) {
cudaEvent_t start, stop;
CUDA_CHECK(cudaEventCreate(&start));
CUDA_CHECK(cudaEventCreate(&stop));
// Warm up this kernel, not just the first kernel in the program. Lazy
// module loading has been the default since CUDA 12.2 on Linux, so the
// first launch of each kernel pays its own load.
for (int i = 0; i < kWarmupRuns; ++i) {
launch();
}
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaEventRecord(start));
for (int i = 0; i < kTimedRuns; ++i) {
launch();
}
CUDA_CHECK(cudaEventRecord(stop));
CUDA_CHECK(cudaEventSynchronize(stop));
CUDA_CHECK(cudaGetLastError());
float ms = 0.0f;
CUDA_CHECK(cudaEventElapsedTime(&ms, start, stop));
CUDA_CHECK(cudaEventDestroy(start));
CUDA_CHECK(cudaEventDestroy(stop));
return ms / kTimedRuns;
}
// Every kernel in part 1 reads the grid once and writes it once, so this is
// the traffic all three rows are divided by. The two stencils leave the
// boundary ring unwritten, which is 0.40 percent of the grid on a 2048
// square; the figure ignores that so the three rows share one denominator.
static double bandwidthGBs(float ms) {
const double bytes = 2.0 * static_cast<double>(kCells) * sizeof(float);
return bytes / (static_cast<double>(ms) * 1.0e-3) / 1.0e9;
}
// Copies the initial condition into both device buffers.
//
// Both, every time. A ping-pong kernel that never writes the boundary reads
// the boundary out of whichever buffer it was handed, so a run that
// initialised only one of them alternates between the real wall temperature
// and whatever the other allocation happened to contain.
static void resetDevice(float* d_a, float* d_b, const float* h_init,
size_t bytes) {
CUDA_CHECK(cudaMemcpy(d_a, h_init, bytes, cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(d_b, h_init, bytes, cudaMemcpyHostToDevice));
}
// Runs kSteps steps of heatStepTiled with the buffers swapped after each one,
// and returns the elapsed milliseconds for the whole loop.
//
// `copyEvery` is the only variable between the three calls in main. Zero
// means the grid never comes back until the run is over. One means every
// step ends with 16 MiB crossing PCIe and the host blocked until it lands.
//
// The events bracket a host loop rather than a single launch, which is
// deliberate: with a synchronous copy in the loop the GPU sits idle while
// the host waits, and that idle time is exactly what this measurement is
// for. timeKernel is not used here because it runs thirteen iterations, and
// a thousand steps is not a microbenchmark.
//
// kSteps is even and the static_assert at the top of the file keeps it that
// way, so the newest grid is back in d_a when the loop ends.
static float runSteps(float* d_a, float* d_b, dim3 grid, dim3 block,
int copyEvery, float* h_snapshot) {
cudaEvent_t start, stop;
CUDA_CHECK(cudaEventCreate(&start));
CUDA_CHECK(cudaEventCreate(&stop));
const size_t bytes = kCells * sizeof(float);
float* cur = d_a;
float* next = d_b;
CUDA_CHECK(cudaEventRecord(start));
for (int step = 0; step < kSteps; ++step) {
heatStepTiled<<<grid, block>>>(cur, next, kN, kStableR);
if (copyEvery > 0 && (step + 1) % copyEvery == 0) {
CUDA_CHECK(
cudaMemcpy(h_snapshot, next, bytes, cudaMemcpyDeviceToHost));
}
float* swap = cur;
cur = next;
next = swap;
}
CUDA_CHECK(cudaEventRecord(stop));
CUDA_CHECK(cudaEventSynchronize(stop));
CUDA_CHECK(cudaGetLastError());
float ms = 0.0f;
CUDA_CHECK(cudaEventElapsedTime(&ms, start, stop));
CUDA_CHECK(cudaEventDestroy(start));
CUDA_CHECK(cudaEventDestroy(stop));
return ms;
}
// Runs kStabilitySteps steps at one value of r and brings the grid back.
// Nothing here is timed; the question is only what the numbers do.
static void runAtR(float* d_a, float* d_b, dim3 grid, dim3 block, float r,
const float* h_init, float* h_out) {
const size_t bytes = kCells * sizeof(float);
resetDevice(d_a, d_b, h_init, bytes);
float* cur = d_a;
float* next = d_b;
for (int step = 0; step < kStabilitySteps; ++step) {
heatStepTiled<<<grid, block>>>(cur, next, kN, r);
float* swap = cur;
cur = next;
next = swap;
}
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_out, cur, bytes, cudaMemcpyDeviceToHost));
}
// Runs kCheckSteps steps of one kernel and compares the result against the
// double-precision reference. Returns EXIT_SUCCESS only when every cell
// agrees and the boundary ring came through untouched.
static int checkAgainstCpu(const char* name, bool tiled, float* d_a, float* d_b,
dim3 grid, dim3 block, const float* h_init,
const double* h_want, float* h_got) {
const size_t bytes = kCells * sizeof(float);
resetDevice(d_a, d_b, h_init, bytes);
float* cur = d_a;
float* next = d_b;
for (int step = 0; step < kCheckSteps; ++step) {
if (tiled) {
heatStepTiled<<<grid, block>>>(cur, next, kN, kStableR);
} else {
heatStepGlobal<<<grid, block>>>(cur, next, kN, kStableR);
}
float* swap = cur;
cur = next;
next = swap;
}
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_got, cur, bytes, cudaMemcpyDeviceToHost));
const size_t bad = firstMismatch(h_got, h_want, kCells, kRelTolerance);
if (bad != kCells) {
std::fprintf(stderr,
"%s wrong after %d steps at row %zu col %zu: got %.9g, "
"want %.9g\n",
name, kCheckSteps, bad / static_cast<size_t>(kN),
bad % static_cast<size_t>(kN),
static_cast<double>(h_got[bad]), h_want[bad]);
return EXIT_FAILURE;
}
// The boundary is Dirichlet, so every cell of the ring has to hold the
// value it started with, exactly. This is the check that catches a run
// which initialised one device buffer and not the other: the wall
// temperature then alternates between the real value and whatever the
// second allocation contained, and it does it in a way that still looks
// like a plausible picture.
for (int k = 0; k < kN; ++k) {
const size_t top = static_cast<size_t>(k);
const size_t bottom = static_cast<size_t>(kN - 1) * kN + k;
const size_t left = static_cast<size_t>(k) * kN;
const size_t right = static_cast<size_t>(k) * kN + (kN - 1);
const size_t ring[4] = {top, bottom, left, right};
for (int e = 0; e < 4; ++e) {
if (h_got[ring[e]] != h_init[ring[e]]) {
std::fprintf(stderr,
"%s moved a boundary cell at index %zu: got "
"%.9g, want %.9g\n",
name, ring[e], static_cast<double>(h_got[ring[e]]),
static_cast<double>(h_init[ring[e]]));
return EXIT_FAILURE;
}
}
}
return EXIT_SUCCESS;
}
int main() {
const int device = 0;
CUDA_CHECK(cudaSetDevice(device));
cudaDeviceProp prop;
CUDA_CHECK(cudaGetDeviceProperties(&prop, device));
const size_t bytes = kCells * sizeof(float);
const dim3 block(kTileX, kTileY);
const dim3 grid(static_cast<unsigned int>((kN - 2 + kTileX - 1) / kTileX),
static_cast<unsigned int>((kN - 2 + kTileY - 1) / kTileY));
const int copyBlocks =
static_cast<int>((kCells + kThreadsPerBlock - 1) / kThreadsPerBlock);
std::vector<float> h_init(kCells);
std::vector<float> h_got(kCells);
std::vector<float> h_snapshot(kCells);
std::vector<double> h_refA(kCells);
std::vector<double> h_refB(kCells);
fillInitial(h_init.data(), kN);
// The reference, kCheckSteps steps in double from the same start.
for (size_t i = 0; i < kCells; ++i) {
h_refA[i] = static_cast<double>(h_init[i]);
}
for (int step = 0; step < kCheckSteps; ++step) {
heatStepCpu(h_refA.data(), h_refB.data(), kN,
static_cast<double>(kStableR));
h_refA.swap(h_refB);
}
float* d_a = nullptr;
float* d_b = nullptr;
CUDA_CHECK(cudaMalloc(&d_a, bytes));
CUDA_CHECK(cudaMalloc(&d_b, bytes));
int status = EXIT_SUCCESS;
float copyMs = 0.0f;
float globalMs = 0.0f;
float tiledMs = 0.0f;
const int copyEvery[3] = {0, kSnapshotEvery, 1};
float loopMs[3] = {0.0f, 0.0f, 0.0f};
std::vector<float> h_final[3];
const float rValues[kNumR] = {kStableR, kBoundaryR, kUnstableR};
float rMax[kNumR] = {0.0f, 0.0f, 0.0f};
bool rFinite[kNumR] = {true, true, true};
if (checkAgainstCpu("heatStepGlobal", false, d_a, d_b, grid, block,
h_init.data(), h_refA.data(),
h_got.data()) != EXIT_SUCCESS) {
status = EXIT_FAILURE;
}
if (status == EXIT_SUCCESS &&
checkAgainstCpu("heatStepTiled", true, d_a, d_b, grid, block,
h_init.data(), h_refA.data(),
h_got.data()) != EXIT_SUCCESS) {
status = EXIT_FAILURE;
}
if (status == EXIT_SUCCESS) {
// Part 1. Every timed launch reads d_a and writes d_b and the two are
// never swapped, so each kernel sees the same input on every run and
// no reset has to sit inside the measurement.
resetDevice(d_a, d_b, h_init.data(), bytes);
copyMs = timeKernel([&] {
copyGrid<<<copyBlocks, kThreadsPerBlock>>>(d_a, d_b, kCells);
});
globalMs = timeKernel(
[&] { heatStepGlobal<<<grid, block>>>(d_a, d_b, kN, kStableR); });
tiledMs = timeKernel(
[&] { heatStepTiled<<<grid, block>>>(d_a, d_b, kN, kStableR); });
// Part 2. Same kernel, same step count, same arithmetic. The only
// thing that moves is how often the grid is copied back to the host.
//
// One discarded run first. Each row below is a single measurement,
// not a timeKernel mean, so without this the first row pays whatever
// clock ramp part 1's short kernels left behind and reads slower
// than rows that do strictly more work.
resetDevice(d_a, d_b, h_init.data(), bytes);
(void)runSteps(d_a, d_b, grid, block, 0, h_snapshot.data());
for (int which = 0; which < 3; ++which) {
resetDevice(d_a, d_b, h_init.data(), bytes);
loopMs[which] = runSteps(d_a, d_b, grid, block, copyEvery[which],
h_snapshot.data());
h_final[which].resize(kCells);
CUDA_CHECK(cudaMemcpy(h_final[which].data(), d_a, bytes,
cudaMemcpyDeviceToHost));
}
// Part 3. Three values of r, same kernel, same initial condition.
for (int which = 0; which < kNumR; ++which) {
runAtR(d_a, d_b, grid, block, rValues[which], h_init.data(),
h_got.data());
rMax[which] = maxAbs(h_got.data(), kCells);
rFinite[which] = allFinite(h_got.data(), kCells);
}
}
// Real branches returning EXIT_FAILURE, not assert(), because CI builds
// Release and NDEBUG deletes assert().
//
// Where the grid lives between steps cannot change the arithmetic, so the
// three runs of part 2 have to produce bit-identical grids. If they do
// not, one of them is reading a buffer somebody else is writing and no
// timing on this page means anything.
if (status == EXIT_SUCCESS) {
for (int which = 1; which < 3; ++which) {
for (size_t i = 0; i < kCells; ++i) {
if (h_final[which][i] != h_final[0][i]) {
std::fprintf(stderr,
"copy policy changed the answer at index "
"%zu: %.9g with copyEvery %d, %.9g with 0\n",
i, static_cast<double>(h_final[which][i]),
copyEvery[which],
static_cast<double>(h_final[0][i]));
status = EXIT_FAILURE;
break;
}
}
if (status != EXIT_SUCCESS) {
break;
}
}
}
// The discrete maximum principle. At r at or under one quarter every
// weight in the update is non-negative and they sum to one, so a new
// value is an average of old ones and the largest value in the grid can
// never grow. This is the same inequality the CFL condition comes from,
// and it is checkable rather than quotable.
if (status == EXIT_SUCCESS) {
const float initialMax = maxAbs(h_init.data(), kCells);
for (int which = 0; which < kNumR; ++which) {
if (rValues[which] > 0.25f) {
continue;
}
if (!rFinite[which] ||
rMax[which] > initialMax * (1.0f + kMaxPrincipleSlack)) {
std::fprintf(stderr,
"maximum principle broken at r = %.3f: max |u| "
"is %.9g after %d steps, and it started at "
"%.9g\n",
static_cast<double>(rValues[which]),
static_cast<double>(rMax[which]), kStabilitySteps,
static_cast<double>(initialMax));
status = EXIT_FAILURE;
}
}
}
// Arithmetic, not a claim about the card: a zero or negative elapsed time
// means the event pair never separated and every figure derived from it
// would be infinite or nonsense.
if (status == EXIT_SUCCESS) {
if (copyMs <= 0.0f || globalMs <= 0.0f || tiledMs <= 0.0f) {
std::fprintf(
stderr, "timing broken in part 1: %.6f, %.6f and %.6f ms\n",
static_cast<double>(copyMs), static_cast<double>(globalMs),
static_cast<double>(tiledMs));
status = EXIT_FAILURE;
}
for (int which = 0; which < 3; ++which) {
if (loopMs[which] <= 0.0f) {
std::fprintf(stderr, "timing broken in part 2 row %d: %.6f\n",
which, static_cast<double>(loopMs[which]));
status = EXIT_FAILURE;
}
}
}
if (status == EXIT_SUCCESS) {
const double mib = static_cast<double>(bytes) / (1024.0 * 1024.0);
const double tiledReads = static_cast<double>(kSharedCells) /
(static_cast<double>(kTileX) * kTileY);
std::printf("GPU: %s (compute capability %d.%d)\n", prop.name,
prop.major, prop.minor);
std::printf(
"Shared memory per block: %zu bytes. Max threads per block: %d\n",
prop.sharedMemPerBlock, prop.maxThreadsPerBlock);
std::printf(
"Grid: %d x %d cells, %.1f MiB per buffer, two buffers on the "
"device\n",
kN, kN, mib);
std::printf(
"Tile: %d x %d output from a %d x %d apron, %zu bytes of shared "
"memory\n",
kTileX, kTileY, kSharedX, kSharedY,
static_cast<size_t>(kSharedCells) * sizeof(float));
std::printf("Interior cells updated per step: %zu\n\n", kInterior);
std::printf("Global reads per updated cell, from the constants\n");
std::printf(" heatStepGlobal %8.3f\n", static_cast<double>(kTaps));
std::printf(" heatStepTiled %8.3f (%d apron cells for %d)\n\n",
tiledReads, kSharedCells, kTileX * kTileY);
std::printf("Part 1: one step, mean of %d runs after %d warm-ups\n",
kTimedRuns, kWarmupRuns);
std::printf("%-22s %11s %9s %12s\n", "kernel", "time (ms)", "GB/s",
"% of copy");
std::printf("%-22s %11s %9s %12s\n", "---------------------",
"----------", "--------", "-----------");
std::printf("%-22s %11.3f %9.1f %12.1f\n", "copyGrid, no stencil",
static_cast<double>(copyMs), bandwidthGBs(copyMs), 100.0);
std::printf("%-22s %11.3f %9.1f %12.1f\n", "heatStepGlobal",
static_cast<double>(globalMs), bandwidthGBs(globalMs),
100.0 * static_cast<double>(copyMs / globalMs));
std::printf("%-22s %11.3f %9.1f %12.1f\n\n", "heatStepTiled",
static_cast<double>(tiledMs), bandwidthGBs(tiledMs),
100.0 * static_cast<double>(copyMs / tiledMs));
std::printf(
"Part 2: %d steps of heatStepTiled. Only the copy policy "
"changes.\n",
kSteps);
std::printf("%-28s %11s %13s %13s %11s\n", "policy", "total (ms)",
"per step (ms)", "to host (GiB)", "vs row 1");
std::printf("%-28s %11s %13s %13s %11s\n",
"---------------------------", "----------", "------------",
"------------", "----------");
for (int which = 0; which < 3; ++which) {
const int every = copyEvery[which];
const double copies =
(every == 0) ? 0.0 : static_cast<double>(kSteps / every);
char label[64];
if (every == 0) {
std::snprintf(label, sizeof(label), "never copied back");
} else if (every == 1) {
std::snprintf(label, sizeof(label), "copied back every step");
} else {
std::snprintf(label, sizeof(label),
"copied back every %d steps", every);
}
std::printf("%-26s %11.3f %13.3f %13.1f %11.1f\n", label,
static_cast<double>(loopMs[which]),
static_cast<double>(loopMs[which]) / kSteps,
copies * mib / 1024.0,
static_cast<double>(loopMs[which] / loopMs[0]));
}
std::printf(
"\nAll three rows above produced a bit-identical final grid. The "
"kernel,\nthe step count and the arithmetic were the same; only "
"the copies moved.\n\n");
std::printf("Part 3: %d steps at three timesteps\n", kStabilitySteps);
std::printf("%8s %10s %12s %22s %9s\n", "r", "1 - 4r", "|1 - 8r|",
"max |u| at the end", "finite");
std::printf("%8s %10s %12s %22s %9s\n", "-------", "---------",
"-----------", "---------------------", "--------");
for (int which = 0; which < kNumR; ++which) {
const double r = static_cast<double>(rValues[which]);
std::printf("%8.3f %10.3f %12.3f %22.6g %9s\n", r, 1.0 - 4.0 * r,
std::fabs(1.0 - 8.0 * r),
static_cast<double>(rMax[which]),
rFinite[which] ? "yes" : "no");
}
std::printf(
"\nThe grid started with max |u| = %.3f. A row whose centre "
"weight 1 - 4r is\nnegative is no longer averaging its "
"neighbours, and |1 - 8r| is what the\ncheckerboard pattern is "
"multiplied by on every one of the %d steps.\n",
static_cast<double>(maxAbs(h_init.data(), kCells)),
kStabilitySteps);
}
CUDA_CHECK(cudaFree(d_a));
CUDA_CHECK(cudaFree(d_b));
return status;
}