CUDA and GPU Parallel Computing Engineering

7.4 Managing Producer Consumer Pipelines with Streams and Events

A producer-consumer pipeline on the GPU is about keeping work moving without forcing everyone to wait for everyone else. The producer stage generates inputs and launches kernels; the consumer stage transforms results and possibly launches follow-up kernels. The key engineering goal is to overlap independent work while preserving correct ordering where dependencies exist.

Core Concepts That Make Pipelines Work

A CUDA stream is an ordered queue of operations. Operations in the same stream execute in issue order; operations in different streams may overlap. Events are lightweight markers that can be recorded into a stream and later used to gate work in another stream.

A practical pipeline uses three rules:

1. Use separate streams for independent stages so the GPU can schedule them concurrently.
2. Record an event at the exact point data becomes ready in the producer stream.
3. Make the consumer wait on that event before it touches the produced data.

This avoids the common mistake of using a global device synchronize, which kills overlap by waiting for everything.

A Systematic Pipeline Layout

Consider processing chunks of a large array. Each chunk goes through:

• Producer: copy or preprocess chunk, then run a “stage A” kernel.
• Consumer: wait for stage A completion, then run “stage B” kernel.

To keep the GPU busy, you process multiple chunks in flight. A typical pattern is double buffering: chunk 0 and chunk 1 alternate roles across streams.

Example: Double-Buffered Chunk Processing

Assume pinned host memory for transfers and device buffers sized for one chunk. Two sets of device buffers allow the producer to work on chunk i+1 while the consumer processes chunk i.

cudaStream_t prod, cons;
cudaEvent_t ready[2];

cudaStreamCreate(&prod);
cudaStreamCreate(&cons);
for (int b = 0; b < 2; b++) cudaEventCreateWithFlags(&ready[b], cudaEventDisableTiming);

for (int i = 0; i < numChunks; i++) {
  int b = i & 1;
  // Producer: H2D + stage A
  cudaMemcpyAsync(d_in[b], h_in[i], bytes, cudaMemcpyHostToDevice, prod);
  stageA<<<grid, block, 0, prod>>>(d_in[b], d_mid[b]);
  cudaEventRecord(ready[b], prod);

  // Consumer: wait + stage B
  cudaStreamWaitEvent(cons, ready[b], 0);
  stageB<<<grid2, block2, 0, cons>>>(d_mid[b], d_out[i]);
}

cudaStreamSynchronize(cons);

This works because the consumer stream is gated per buffer index. The event is recorded after stage A writes d_mid[b], and the consumer waits before reading it.

Choosing Where to Record Events

Record events at boundaries where the consumer truly needs completion. Recording too early can cause the consumer to read partially written data; recording too late reduces overlap because the consumer waits longer than necessary.

A good heuristic is: record immediately after the last operation in the producer stream that writes the shared data region for that chunk.

Handling Transfers and Compute Together

If you also overlap host-device transfers, use a dedicated transfer stream. The producer stream can depend on transfer completion via an event, then run stage A. The consumer stream can depend on stage A completion via another event.

cudaStream_t xfer, prod, cons;
cudaEvent_t xferDone[2], stageADone[2];

for (int i = 0; i < numChunks; i++) {
  int b = i & 1;
  cudaMemcpyAsync(d_in[b], h_in[i], bytes, cudaMemcpyHostToDevice, xfer);
  cudaEventRecord(xferDone[b], xfer);

  cudaStreamWaitEvent(prod, xferDone[b], 0);
  stageA<<<grid, block, 0, prod>>>(d_in[b], d_mid[b]);
  cudaEventRecord(stageADone[b], prod);

  cudaStreamWaitEvent(cons, stageADone[b], 0);
  stageB<<<grid2, block2, 0, cons>>>(d_mid[b], d_out[i]);
}

The transfer stream isolates copy scheduling, while events preserve correctness.

Common Failure Modes and Fixes

• Missing cudaStreamWaitEvent: the consumer may start before data is ready. Fix by waiting on the correct event for the correct buffer index.
• Reusing buffers too soon: if you only have one buffer but multiple chunks are in flight, you’ll overwrite data. Fix by using double buffering or more.
• Recording the event before the write completes: if stage A is asynchronous and you record too early, the event won’t protect the consumer. Fix by recording after the kernel launch in the same stream.

A Quick Mental Model

Think of each event as a “data ready ticket” tied to a specific buffer slot. Streams decide when operations are allowed to run; events decide which operations are allowed to run based on data readiness. When you keep those responsibilities separate, the pipeline stays both fast and correct.

7.5 Overlapping Compute and Transfers with Asynchronous Execution

Overlapping compute with data transfers is mostly about one thing: keeping the GPU busy while the CPU prepares the next chunk. In CUDA, the practical lever is asynchronous copies issued into streams, paired with kernels that run in those same streams. If the transfers and kernels are independent, the hardware can execute them concurrently.

Core Idea and Preconditions

1. Use multiple streams so the GPU has independent work queues.
2. Use asynchronous operations like cudaMemcpyAsync so the host thread does not block.
3. Ensure the data path supports overlap: host memory should be pinned (page-locked), otherwise transfers may serialize.
4. Avoid implicit synchronization: calls that force device-wide ordering (for example, certain blocking syncs) will collapse overlap.

A simple mental model: stream A handles “copy chunk i” then “compute chunk i”, while stream B handles “copy chunk i+1” then “compute chunk i+1”. When chunk i is computing, chunk i+1 can be copying.

Double Buffering with Two Streams

Chunk your workload so each chunk can be transferred and processed independently. Allocate two device buffers for inputs and outputs, and alternate which buffer each stream uses.

Key rule: operations in the same stream execute in order, but operations in different streams may overlap. That means you can safely enqueue copy-then-kernel pairs per stream without extra synchronization.

Example: Enqueue Copy and Kernel per Stream

// Assume pinned host buffers h_in, h_out and device buffers d_in[2], d_out[2]
cudaStream_t s[2];
cudaStreamCreate(&s[0]);
cudaStreamCreate(&s[1]);

for (int i = 0; i < numChunks; ++i) {
  int b = i & 1;                 // double buffer index
  size_t off = (size_t)i * chunkElems;
  size_t n = min(chunkElems, totalElems - off);

  cudaMemcpyAsync(d_in[b], h_in + off, n*sizeof(float),
                  cudaMemcpyHostToDevice, s[b]);

  myKernel<<<gridFor(n), block, 0, s[b]>>>(d_in[b], d_out[b], n);

  cudaMemcpyAsync(h_out + off, d_out[b], n*sizeof(float),
                  cudaMemcpyDeviceToHost, s[b]);
}

cudaDeviceSynchronize();

This pattern overlaps transfers for chunk i+1 with computation for chunk i, as long as the GPU and PCIe/NVLink path can run copy engines concurrently with kernels.

When You Need Dependencies: Events

Sometimes chunk i computation depends on results produced by chunk i-1 on the GPU. In that case, you cannot rely on stream ordering across streams. Use events to express a dependency without forcing a global sync.

cudaEvent_t done[2];
for (int b = 0; b < 2; ++b) cudaEventCreateWithFlags(&done[b], cudaEventDisableTiming);

// After enqueuing kernel for chunk i in stream s[b], record an event.
// Before starting dependent work for chunk i+1 in stream s[b2], make it wait.
// cudaEventRecord(done[b], s[b]);
// cudaStreamWaitEvent(s[b2], done[b], 0);

Record the event in the stream where the producing kernel runs, then make the consuming stream wait on that event. This keeps unrelated work overlapped while preserving correctness.

Avoiding Overlap Killers

• Synchronous copies (cudaMemcpy) block the host and often serialize the timeline.
• Unpinned host memory can prevent true overlap because the driver may stage through temporary buffers.
• Device-wide synchronization inside the loop collapses concurrency.
• Reusing the same staging buffers across streams without ordering guarantees can cause data races. Double buffering exists to prevent that.

Measuring Overlap Correctly

Use a timeline view to confirm concurrency. Look for regions where copy operations and kernels overlap in time. If you see a strict “copy then compute then copy then compute” pattern, check pinned memory, stream usage, and whether any sync calls occur between enqueues.

Practical Checklist for Engineering Use

1. Split work into chunks that fit staging buffers.
2. Allocate pinned host buffers for inputs and outputs.
3. Create at least two streams and two sets of device staging buffers.
4. Enqueue MemcpyAsync → kernel → MemcpyAsync in the same stream.
5. Use events only when cross-stream dependencies exist.
6. Validate correctness first with small sizes, then scale chunk counts.

Done carefully, overlap turns “waiting for data” into “useful work”, and the GPU spends more time computing than waiting for the CPU to catch up.

8.1 Pinned Memory Allocation and Transfer Optimization

Pinned (page-locked) host memory reduces transfer overhead by preventing the operating system from paging host buffers out while a DMA engine is moving data. The result is more predictable throughput and better overlap with GPU work when you use streams correctly. The trick is to pin only what you need, for as long as you need it, and to structure transfers so the GPU is never waiting on a slow host.

Foundational Concepts That Drive the Design

Pinned memory changes two things that matter for engineering:

1. Transfer path stability: pageable memory may require staging through an internal buffer, adding extra copies and latency.
2. Asynchronous behavior: pinned buffers can participate in truly asynchronous cudaMemcpyAsync calls, which lets you overlap copies with kernel execution.

A practical mental model: pageable memory is like sending a package that might be temporarily rerouted; pinned memory is like handing the courier a fixed address and a ready-to-ship crate.

Allocating Pinned Memory Without Creating New Problems

Pinning is not free. Pinned pages reduce the flexibility of the OS memory manager, so pinning too much can hurt overall system performance. A good default is to pin buffers that:

• are reused across many iterations (e.g., per-timestep fields in a simulation),
• are large enough that copy overhead dominates,
• are contiguous and naturally aligned.

For small one-time transfers, pageable memory is often fine because the overhead of pinning and unpinning can outweigh the gains.

Example: Reusable Pinned Buffers with Streamed Transfers

// Allocate pinned host buffers once
cudaHostAlloc(&h_in,  bytes, cudaHostAllocDefault);
cudaHostAlloc(&h_out, bytes, cudaHostAllocDefault);

cudaStream_t s;
cudaStreamCreate(&s);

for (int iter = 0; iter < iters; ++iter) {
  fill_host_input(h_in, iter);

  cudaMemcpyAsync(d_in,  h_in,  bytes, cudaMemcpyHostToDevice, s);
  myKernel<<<grid, block, 0, s>>>(d_in, d_out);
  cudaMemcpyAsync(h_out, d_out, bytes, cudaMemcpyDeviceToHost, s);
}

cudaStreamSynchronize(s);

cudaFreeHost(h_in);
cudaFreeHost(h_out);
cudaStreamDestroy(s);

This pattern works because the host buffers are pinned and the operations are issued on a stream. If your kernel and copies are independent, the runtime can overlap them; if not, the stream ordering still guarantees correctness.

Transfer Geometry That Prevents Hidden Slowdowns

Pinned memory helps, but it cannot fix inefficient transfer geometry. Keep these rules in mind:

• Prefer contiguous layouts: copying strided slices forces extra work and may reduce effective bandwidth.
• Transfer in chunks that match your pipeline: if you stream through a large array, split it into chunks that are big enough to amortize overhead but small enough to overlap with compute.
• Use the right size types: transferring float arrays as bytes is fine, but ensure the byte count matches the element count times sizeof(T).

Example: Chunked Transfers for Overlap

size_t chunkElems = 1<<20; // tune for your workload
size_t chunkBytes = chunkElems * sizeof(float);

for (int k = 0; k < numChunks; ++k) {
  size_t offset = k * chunkElems;

  cudaMemcpyAsync(d_in + offset, h_in + offset,
                  chunkBytes, cudaMemcpyHostToDevice, s[k%2]);
  myKernel<<<grid, block, 0, s[k%2]>>>(d_in + offset, d_out + offset);
  cudaMemcpyAsync(h_out + offset, d_out + offset,
                  chunkBytes, cudaMemcpyDeviceToHost, s[k%2]);
}

Using two streams (s[0] and s[1]) implements a simple double-buffering scheme: while one chunk is copying, another can be computing.

Measuring Whether Overlap Actually Happens

Pinned memory is only useful if it improves end-to-end time. Measure three durations separately:

• host-to-device copy time,
• kernel time,
• device-to-host copy time.

Then check whether the timeline shows overlap. If copies and kernels are serialized, the issue is usually one of these:

• you reused the same stream for dependent stages,
• the kernel depends on data that is still transferring,
• the chunk sizes are too small to hide latency,
• the GPU is memory-bound and cannot run concurrently with copy.

Lifecycle Management That Keeps Systems Stable

Pinning should be treated like reserving a resource, not like a casual convenience. Allocate pinned buffers once per buffer lifetime, reuse them across iterations, and free them when the owning scope ends. If you repeatedly allocate and free pinned memory inside a loop, you pay overhead and risk memory pressure from pinned pages.

A clean engineering rule: pin at initialization, reuse during the hot loop, unpin at shutdown. Your transfer code stays simple, and your system stays predictable.

8.2 Stream Design for Concurrency Across Independent Work

Streams let you overlap work that does not depend on each other: independent kernel launches, independent memory transfers, and even independent stages of a pipeline. The trick is to design a schedule that respects dependencies while keeping the GPU busy and the PCIe bus useful.

Foundational Concepts for Stream Concurrency

A CUDA stream is an ordered queue of operations. Operations in the same stream execute in issue order; operations in different streams may overlap if the hardware has resources and if the operations are compatible with overlap. “Compatible” usually means: no cross-stream dependency that forces serialization, and enough copy engines or compute capacity to run concurrently.

Think in three layers:

1. Dependency graph: what must happen before what.
2. Placement: which stream each operation goes into.
3. Resource fit: whether the GPU can run those operations at the same time.

A practical rule: start with a dependency graph, then assign each independent chain to its own stream, and finally verify overlap with profiling.

Designing a Dependency Graph That Actually Helps

Suppose you have three independent tasks per iteration: A computes, B transfers input, and C computes on transferred data. A common mistake is to put everything into one stream, which guarantees correctness but prevents overlap.

Instead, define dependencies explicitly:

• B must finish before C starts.
• A has no dependency on B or C.

This yields two streams:

• Stream 0: A.
• Stream 1: B then C.

If you also have output transfers, you can add a third stream for independent output copies, but only if those copies do not depend on the same buffers being written by kernels in conflicting ways.

Stream Assignment Patterns That Scale Beyond Toy Examples

Pattern 1: One Stream per Independent Chain

Use separate streams for independent chains of work. This is the simplest approach and often enough for engineering workloads.

Pattern 2: Double Buffering for Iterative Pipelines

For iterative solvers or time-stepping, you can overlap “iteration i transfer” with “iteration i compute” by alternating between two buffer sets. Stream assignment becomes systematic:

• Stream for H2D copy of buffer set 0 while compute uses buffer set 1.
• Swap roles next iteration.

Pattern 3: Staged Pipelines with Events

When you need cross-stream coordination, use events to avoid global synchronization. Record an event at the end of a producer operation, then make the consumer stream wait on that event.

Concrete Example with Two Streams and Double Buffering

Assume each iteration needs:

• Copy input chunk in[k] to device.
• Run kernel on that chunk.
• Copy results back.

You can use two device buffers d_in[2] and d_out[2] and two streams s0 and s1.

cudaStream_t s0, s1;
cudaStreamCreate(&s0);
cudaStreamCreate(&s1);

for (int i = 0; i < iters; ++i) {
  int buf = i & 1;
  cudaStream_t s = (buf == 0) ? s0 : s1;

  cudaMemcpyAsync(d_in[buf], h_in[i], bytes, cudaMemcpyHostToDevice, s);
  myKernel<<<grid, block, 0, s>>>(d_in[buf], d_out[buf]);
  cudaMemcpyAsync(h_out[i], d_out[buf], bytes, cudaMemcpyDeviceToHost, s);
}

cudaStreamSynchronize(s0);
cudaStreamSynchronize(s1);

This works because each stream uses its own buffer set, so operations do not fight over the same memory. The GPU can overlap copy and compute when hardware supports it.

Avoiding the Usual Concurrency Traps

1. Accidental serialization via shared buffers: if two streams write the same device pointer, you reintroduce ordering constraints. Use separate buffers per stream or per pipeline stage.
2. Host memory not pinned: cudaMemcpyAsync is only truly asynchronous for page-locked host memory. If host buffers are pageable, overlap shrinks dramatically.
3. Implicit sync from default stream usage: mixing the legacy default stream with other streams can cause surprising ordering. Use explicit streams consistently.
4. Too many tiny operations: concurrency helps, but launch overhead still matters. Batch work when it is safe.

A Systematic Checklist Before You Trust the Timeline

• Each independent chain has its own stream.
• Any cross-stream dependency is represented with events or buffer separation.
• Host buffers for async copies are pinned.
• You never reuse the same device buffer in two streams without a defined ordering.
• You validate overlap using a timeline view, not just total runtime.

When these conditions hold, stream design becomes less about “more streams” and more about a clean schedule that matches the hardware’s ability to run copies and kernels at the same time.

8.3 Unified Memory Usage with Explicit Prefetching and Guidance

Unified Memory (UM) lets you allocate data once and access it from both the CPU and GPU. The runtime migrates pages on demand, which is convenient but can be slow when the first touch happens inside a kernel. Explicit prefetching and guidance moves that cost to predictable points and nudges the runtime toward the access pattern you actually intend.

Core Concepts You Need Before Touching Prefetch

UM is page-based. Your arrays are split into pages, and each page has a current “residency” location. When a CPU or GPU first touches a page, the runtime may migrate it. If that first touch occurs during a kernel, you pay migration latency while the GPU is waiting.

Two practical rules follow.

1. Make first touch happen on purpose. Prefetch pages to the device before launching kernels that will read or write them.
2. Tell the runtime what you plan to do. Guidance hints help the runtime choose whether pages should stay on the GPU, be treated as read-mostly, or be migrated back.

A Systematic Workflow for Explicit Prefetching

Step 1: Allocate with UM

Use cudaMallocManaged for arrays shared across CPU and GPU. Keep allocations large enough to amortize overhead, but not so large that you thrash memory.

Step 2: Decide Residency per Phase

Most engineering workloads have phases: initialization on CPU, compute on GPU, and optional post-processing on CPU. Prefetching should match those phases.

• Initialization phase: prefetch to CPU (or just let CPU touch first).
• Compute phase: prefetch to the active GPU.
• Post-processing: prefetch back to CPU if you will read results there.

Step 3: Prefetch Before Kernel Launch

Prefetching is a host-side operation. Call it after you finish CPU writes and before you launch the kernel that will consume the data.

Step 4: Use Guidance for Access Patterns

Guidance is not magic; it’s a hint. Still, it can reduce unnecessary migrations and improve cache behavior.

Common guidance choices:

• Read-mostly: if the GPU reads a dataset repeatedly while the CPU rarely writes it.
• Preferred location: if a region should stay on the GPU for the duration of compute.

Step 5: Synchronize at the Right Boundaries

Prefetching and kernel execution are asynchronous with respect to the host unless you synchronize. Use cudaDeviceSynchronize() or stream synchronization when you must ensure residency before accessing results.

Example: Prefetching a Managed Array for a Compute Phase

// Allocate managed memory
float* a = nullptr;
size_t bytes = N * sizeof(float);
cudaMallocManaged(&a, bytes);

// CPU initialization
for (size_t i = 0; i < N; i++) a[i] = 1.0f;

// Prefetch to GPU before kernel touches it
int device = 0;
cudaGetDevice(&device);
cudaMemPrefetchAsync(a, bytes, device, 0);

// Kernel launch uses a on the GPU
myKernel<<<grid, block>>>(a, N);

// Ensure kernel completion before CPU reads
cudaDeviceSynchronize();

// Prefetch back if CPU will read results
cudaMemPrefetchAsync(a, bytes, cudaCpuDeviceId, 0);

This pattern avoids the “GPU first touch” penalty. Without prefetching, the first kernel access can trigger page migration while the GPU is already scheduled, which often shows up as a stall.

Example: Guidance for Read-Mostly Data

If you have a constant-ish lookup table used by many kernels, guidance can help keep it on the GPU.

float* table = nullptr;
cudaMallocManaged(&table, tableBytes);

// Initialize on CPU
initTable(table, M);

// Hint that GPU will read it frequently
cudaMemAdvise(table, tableBytes, cudaMemAdviseSetReadMostly, 0);

// Prefetch to GPU for the compute phase
cudaMemPrefetchAsync(table, tableBytes, 0, 0);

// Launch kernels that read table
for (int k = 0; k < K; k++)
  myLookupKernel<<<grid, block>>>(table, M);

cudaDeviceSynchronize();

The key is that guidance should match reality: if the CPU writes to table often, read-mostly hints can backfire by encouraging the runtime to keep stale pages around.

Practical Checks That Prevent Common Mistakes

• Bounds and indexing still matter. UM doesn’t protect you from out-of-range accesses; it just makes memory movement more forgiving.
• Measure first-touch timing. If you see stalls, confirm that prefetch happens after CPU writes and before the kernel.
• Avoid mixing conflicting advice. Preferred location and read-mostly hints should be consistent with your phase plan.

Integrated Summary for Engineering Use

Treat UM as a page-migration system you can steer. Prefetching sets the stage so kernels don’t pay migration latency mid-flight. Guidance refines the runtime’s decisions by reflecting your actual access pattern. When you combine both with phase-based residency planning and correct synchronization, UM becomes predictable enough for high-performance engineering rather than just convenient for demos.

8.4 Minimizing Launch Overhead with Kernel Fusion Techniques

Kernel launch overhead is the tax you pay every time the host asks the GPU to start work. For large kernels, the tax is small compared to execution time; for small or highly segmented workloads, it can dominate. Kernel fusion reduces the number of launches by combining multiple operations into one kernel, but it also changes resource usage and can make debugging harder. The engineering goal is to fuse safely and only when it improves end-to-end time.

Foundational Model of Launch Overhead

A typical pipeline looks like this: host prepares arguments, launches a kernel, GPU schedules blocks, and the host later synchronizes or enqueues dependent work in a stream. Launch overhead includes CPU-side work (building launch parameters, driver calls) and GPU-side setup (queueing and scheduling). If your workload is split into many kernels—say, elementwise transform, then reduction, then another transform—you may spend more time launching than computing.

A quick sanity check is to time a single fused candidate against the original sequence using the same input sizes and stream behavior. If the fused kernel is not faster, you likely increased register pressure, reduced occupancy, or introduced extra memory traffic.

Fusion Decision Rules That Prevent Accidental Slowdowns

Start with operations that:

• Are elementwise or have compatible iteration spaces (same indexing scheme).
• Can share intermediate values without writing them to global memory.
• Have simple dependency ordering (producer-to-consumer within the same thread or block).
• Do not require different synchronization scopes beyond what a single kernel can provide.

Avoid fusing when:

• The fused kernel would require large shared memory buffers for multiple stages.
• The stages have very different memory access patterns that would fight each other.
• You need global synchronization between stages; a single kernel cannot synchronize across all blocks.

Practical Fusion Patterns with Examples

Example: Fuse Two Elementwise Kernels

Suppose you currently do:

1. y = a*x + b
2. z = y*y + c

Two launches write y to global memory and read it back. Fusion computes z directly:

__global__ void fused_ax2(const float* x, float* z, float a, float b, float c, int n){
  int i = blockIdx.x * blockDim.x + threadIdx.x;
  if(i >= n) return;
  float y = a * x[i] + b;
  z[i] = y * y + c;
}

This reduces global memory traffic and eliminates one launch. The fused kernel uses a couple extra registers for y, but the memory savings often outweigh it.

Example: Fuse Stencil Compute with Local Post-Processing

For a 2D stencil, you might compute u_next and then apply a clamp or boundary adjustment in a second kernel. If the post-processing depends only on the computed value for each output element, fuse it into the same kernel. The key is to keep halo loads and shared memory tiling unchanged, then apply the post-step before writing the final output.

Example: Fuse Reduction into a Block-Local Pipeline

Reductions are trickier because they often require multiple passes. A safe fusion target is block-local reduction steps that can be completed within one kernel stage. For instance, you can fuse “load + partial reduction + write partial sums” into one kernel, then keep the final global reduction as a separate kernel.

This is not “one kernel for everything”; it’s “fuse what fits the synchronization boundary.”

Systematic Workflow for Choosing and Validating Fusion

1. Measure the baseline: time the original sequence with the same stream and synchronization points.
2. Identify intermediate writes: if a kernel writes an intermediate that the next kernel immediately reads, fusion is a prime candidate.
3. Prototype one fusion step: fuse two adjacent kernels first, not the whole pipeline.
4. Check resource pressure: if the fused kernel uses many more registers or shared memory, occupancy may drop and erase the benefit.
5. Validate correctness: test boundaries, odd sizes, and alignment-sensitive cases.
6. Re-measure end-to-end: compare total time, not just kernel time, because fewer launches can change scheduling and overlap.

A Small Engineering Example with a Realistic Outcome

If your pipeline is transform -> reduce -> transform, you might fuse the two transforms but keep the reduction separate. The fused transform kernel removes one launch and avoids writing an intermediate array. The reduction remains separate because it typically needs a different access pattern and may require multiple stages. This hybrid approach often beats both extremes: “no fusion” and “everything fused.”

Practical Debugging Without Losing Your Mind

When fused kernels misbehave, you need a way to isolate which stage is wrong. A practical technique is to temporarily add a compile-time switch that writes intermediate results to global memory for a small test size. Once correctness is confirmed, remove the intermediate writes and keep the fused fast path.

// Pseudocode sketch
// if(DEBUG) write y to global; else keep y in registers

This keeps the production kernel lean while giving you a controlled view of the intermediate values during development.

8.5 Building End-to-End Pipelines for Iterative Solvers

An iterative solver on the GPU is rarely “one kernel.” It’s a pipeline: prepare inputs, run compute kernels, reduce or communicate partial results, check convergence, and repeat. The engineering goal is to keep the GPU busy while keeping data movement predictable and correctness easy to verify.

Pipeline Mental Model

Think in stages with explicit data ownership. A typical Krylov-style loop has:

1. Precompute: build or update vectors and any operator parameters.
2. Apply Operator: compute (y = A x) (often sparse or stencil-like).
3. Vector Updates: combine vectors using axpy-like operations.
4. Reductions: compute dot products and norms for step sizes and stopping criteria.
5. Convergence Check: decide whether to continue.
6. Swap and Repeat: update (x), residuals, and any cached intermediates.

Each stage should either be fully on-device or have a clearly bounded host interaction. If the host needs a scalar for convergence, design the pipeline so the scalar arrives without stalling the entire GPU.

Data Flow and Memory Strategy

Start by deciding where each vector lives. A practical rule: keep all iteration vectors in device memory and only move scalars (like a single norm) to the host.

Use a consistent layout for vectors (flat 1D arrays) and for operators (CSR for sparse, structured grids for stencils). That consistency makes kernel indexing and correctness checks less error-prone.

For reductions, prefer a two-phase approach: kernel produces partial sums per block, then a second kernel (or a small reduction kernel) finalizes the result. This avoids forcing every thread to contend for a single global location.

Streams and Overlap Without Guesswork

Use multiple CUDA streams only when you can name what overlaps. A common pattern:

• Stream A runs the operator and vector update kernels.
• Stream B runs reductions that depend on completed updates.
• Host-side convergence check happens after the reduction scalar is ready.

To avoid accidental race conditions, enforce dependencies with events: record an event after the last kernel that produces inputs for a reduction, then make the reduction stream wait on that event.

Example Pipeline for a Simple Iterative Method

Below is a minimal structure for a loop that applies an operator, updates vectors, computes a norm, and stops.

for (int iter = 0; iter < maxIters; ++iter) {
  // 1) Apply operator: y = A * x
  applyOperator<<<grid, block, 0, streamA>>>(A, x, y);

  // 2) Vector update: r = b - y
  axpy<<<grid, block, 0, streamA>>>(b, y, r);

  // 3) Reduction: norm = ||r||
  cudaEventRecord(evReady, streamA);
  cudaStreamWaitEvent(streamB, evReady, 0);

  partialReduce<<<gridR, blockR, 0, streamB>>>(r, partial);
  finalizeReduce<<<1, 256, 0, streamB>>>(partial, normDev);

  // 4) Copy scalar for convergence
  cudaMemcpyAsync(&normHost, normDev, sizeof(double),
                  cudaMemcpyDeviceToHost, streamB);
  cudaStreamSynchronize(streamB);

  if (normHost < tol) break;

  // 5) Update x using r and step size
  updateX<<<grid, block, 0, streamA>>>(x, r, step);
}

The key detail is that only normHost is synchronized. Everything else stays asynchronous until it must be ordered.

Convergence Check Without Pipeline Stalls

If you synchronize every iteration on the host, you’ll often throttle throughput. A better approach is to keep the convergence scalar on the host only as needed and ensure the GPU work for the next iteration doesn’t depend on the host decision.

For methods where the next iteration’s kernels require the convergence decision, you can still reduce stall time by:

• making the reduction scalar computation fast and predictable,
• using a single scalar copy per iteration,
• keeping the rest of the pipeline in device memory.

Practical Checklist for Integrated Correctness

Before optimizing, verify the pipeline as a whole:

• Confirm the residual definition matches the math you intend (signs and scaling matter).
• Validate that reductions compute the same scalar as a CPU reference for a small problem.
• Check that vector swaps are consistent with kernel expectations (no stale pointers).
• Ensure every kernel uses the same indexing convention and bounds logic.

Once correctness is stable, optimization becomes a matter of moving the bottleneck: memory bandwidth in operator kernels, reduction efficiency in dot products, or launch ordering in the pipeline. The pipeline structure above keeps those bottlenecks visible instead of hiding them behind accidental synchronization.

9.1 Selecting Multi-GPU Work Partitioning Strategies

Multi-GPU scaling starts with a boring question: how do you split the work so each GPU stays busy without drowning in communication? The best partitioning strategy depends on (1) where the data lives, (2) how much boundary data must move, and (3) whether the workload is naturally uniform or varies across the domain.

Foundational Partitioning Choices

Domain Decomposition

Domain decomposition splits the problem space into subdomains, one per GPU. Each GPU computes its interior cells and exchanges boundary (halo) regions with neighbors. This fits stencils, PDE solvers, and grid-based simulations.

A practical rule: if your computation needs nearby neighbors, domain decomposition is usually the cleanest mapping. If your computation is mostly independent per item, other strategies may reduce halo traffic.

Data Decomposition

Data decomposition splits arrays or batches across GPUs while keeping the same overall index space. Each GPU processes its chunk and later combines results. This fits embarrassingly parallel workloads like per-particle kernels, per-sample inference, or batched transforms.

A practical rule: if you can aggregate results with a small number of reductions, data decomposition can be efficient.

Task Decomposition

Task decomposition assigns different kernels or stages to different GPUs. This can work when stages have different costs, but it often complicates synchronization and data movement because intermediate results must move between devices.

A practical rule: use task decomposition when stages are clearly separable and the intermediate data is small or already replicated.

Strategy Selection Criteria

Communication Volume Versus Compute

For domain decomposition, communication scales with halo surface area, while compute scales with subdomain volume. If you double GPUs by slicing thinner slabs, halo-to-compute ratio rises and performance can flatten.

For data decomposition, communication scales with the size of the final reduction or gather. If you can keep reductions small, scaling is often smoother.

Load Balance

Even with perfect partitioning, real workloads can be uneven. For example, adaptive meshes or irregular sparsity produce hotspots. If one GPU finishes early, others wait at synchronization points.

A practical rule: measure per-GPU time for a few representative inputs, then adjust partition sizes or use finer-grained chunking.

Data Locality and Layout

Partitioning should match memory access patterns. If each GPU repeatedly touches the same neighbor regions, keep those regions contiguous in memory and minimize strided access.

A practical rule: when splitting multidimensional arrays, partition along the fastest-changing dimension only if it preserves coalesced access; otherwise partition along a dimension that keeps contiguous tiles per GPU.

Example: 3D Stencil Domain Decomposition

Assume a 3D grid of size Nx * Ny * Nz and a 7-point stencil. With G GPUs, a common choice is to split along Nz into slabs.

• Each GPU owns a slab of Nz/G planes.
• Each iteration computes interior planes.
• Before computing boundary planes, GPUs exchange one halo plane with the neighbor GPUs.

If you instead split into 2D blocks (e.g., Nx and Nz), each GPU has more neighbors, and halo exchange grows. Sometimes that’s still worth it because each GPU’s subdomain becomes smaller and fits better in cache/shared memory, but you should expect more communication endpoints.

A practical rule: start with 1D slabs for simplicity, then move to 2D/3D blocks only if load balance or memory constraints demand it.

Example: Batched Independent Work with Reduction

Suppose you have B independent samples, each producing a scalar loss. With G GPUs, split the batch into chunks of size B/G.

• GPU i computes losses for samples [i*chunk, (i+1)*chunk).
• Each GPU produces a partial sum.
• A final reduction combines partial sums.

This avoids halo exchange entirely. The only communication is the reduction of G scalars (or small arrays if you compute multiple metrics).

A practical rule: if your final aggregation is small, data decomposition is hard to beat.

Example: Handling Uneven Workloads

Consider a simulation where some regions are more active (e.g., higher refinement or more nonzeros). If you use fixed geometric partitions, some GPUs will spend more time.

A practical approach is to partition by work units rather than pure geometry. For instance, you can map tiles to GPUs using a simple greedy assignment based on estimated tile cost, then still exchange halos based on tile adjacency.

This keeps the halo logic intact while improving load balance.

Implementation Checklist

• Choose decomposition type based on dependency structure.
• Estimate communication volume: halo surface area for domain decomposition, reduction size for data decomposition.
• Ensure partitions preserve coalesced access and contiguous tiles.
• Validate load balance with timing on representative inputs.
• Keep synchronization minimal: overlap halo exchange with interior computation when possible.

When these pieces align, scaling stops being a guessing game and becomes an engineering exercise: measure, adjust partition geometry, and confirm that the communication pattern matches the workload’s dependency graph.

9.2 Peer To Peer Access And Topology Aware Data Paths

Peer to peer (P2P) access lets one GPU read or write another GPU’s memory without staging through host memory. That matters because host staging adds latency and burns PCIe bandwidth. Topology aware data paths go one step further: they choose the communication route that matches how your GPUs are physically connected, so you don’t accidentally route traffic through a slower hop.

Foundations: What P2P Actually Changes

When P2P is enabled, a kernel on GPU A can issue global memory loads/stores to pointers allocated on GPU B. The CUDA runtime maps those pointers to the correct device address space and uses the hardware path between devices. If P2P is not enabled, the same pointer usage is invalid or forces you to copy through host memory (depending on your API choices).

Two practical constraints drive the design:

• Reachability: Not every GPU pair can support direct access. Some systems require going through the host.
• Path quality: Even when P2P works, the effective bandwidth and latency depend on whether the GPUs share a PCIe switch, are connected via NVLink, or must traverse multiple links.

Enabling P2P Correctly

Start by checking whether each GPU pair supports P2P. In engineering terms, you want a matrix that tells you which pairs are safe and which are likely to be slow.

int canAccess = 0;
cudaDeviceCanAccessPeer(&canAccess, devA, devB);
if (canAccess) {
  cudaSetDevice(devA);
  cudaDeviceEnablePeerAccess(devB, 0);
}

A common mistake is enabling access only in one direction. Many applications need both directions for halo exchange or all-reduce style patterns, so treat the matrix as directional and enable what you need.

Topology Aware Routing: Choosing the Right Pairings

Topology aware routing begins with the observation that “GPU 0 to GPU 3” is not a single thing. It can be:

• Direct fast link (e.g., NVLink-like connectivity)
• Same PCIe switch (often good)
• Different switches (often slower)
• Host-mediated (worst case)

Your goal is to map communication-heavy neighbors in your domain decomposition onto the best-connected GPU pairs.

Practical Example: Halo Exchange Neighbor Selection

Suppose you decompose a 3D grid across 4 GPUs arranged in a logical ring for simplicity. If you blindly connect neighbors as (0↔1, 1↔2, 2↔3, 3↔0), you may place the heaviest halo exchange on a slow path.

Instead:

1. Query topology and build a “link quality” score per GPU pair.
2. Assign logical neighbors so that each GPU’s two halo partners use the best available links.
3. Keep the domain mapping consistent with your indexing so you don’t add extra copies.

Even if your ring stays a ring, the GPU IDs behind each logical position can change.

Measuring and Verifying the Path

Topology queries tell you what the system can do, not what it will do under your workload. Verification should be part of the engineering loop.

A simple bandwidth test uses a large buffer on GPU B and measures copy or load/store throughput from GPU A. You don’t need a full benchmark suite to catch obvious problems like host-mediated transfers.

When you measure, keep these details consistent:

• Use pinned host memory only if you’re testing host staging; otherwise focus on device-to-device.
• Use the same transfer size and alignment.
• Run multiple iterations and take a stable median.

Example: Two-GPU P2P Halo Exchange Skeleton

The key idea is to allocate device buffers on each GPU, enable P2P, then exchange boundary slabs using device-to-device copies or direct kernel reads.

// Setup
cudaSetDevice(devA);
float* bufA = allocOnDevice(devA, bytes);

cudaSetDevice(devB);
float* bufB = allocOnDevice(devB, bytes);

// Enable P2P both ways if needed
cudaDeviceCanAccessPeer(&okAB, devA, devB);
cudaDeviceCanAccessPeer(&okBA, devB, devA);
if (okAB) { cudaSetDevice(devA); cudaDeviceEnablePeerAccess(devB, 0); }
if (okBA) { cudaSetDevice(devB); cudaDeviceEnablePeerAccess(devA, 0); }

// Exchange using device-to-device copy
cudaSetDevice(devA);
cudaMemcpyPeerAsync(dstOnA, devA, srcOnB, devB, bytes, streamA);

To avoid synchronization stalls, pair this with events: record completion on the sender stream, then wait on the receiver stream before launching the compute kernel that consumes the halo.

Engineering Checklist for Topology Aware P2P

• Build a pairwise capability matrix using device reachability checks.
• Enable peer access only for pairs you actually use.
• Choose domain neighbor mappings based on link quality, not GPU IDs.
• Use stream-based copies and event synchronization to overlap exchange with computation.
• Validate with a targeted bandwidth/latency test so you catch slow paths early.

With these steps, P2P becomes a controlled tool rather than a “works on my machine” feature. Your communication pattern stays correct, and your performance stops being a mystery box.

9.3 Using NCCL for Collective Operations Across Devices

When you scale a CUDA workload across multiple GPUs, you quickly run into the same question: how do you move and combine data efficiently without turning every iteration into a host-driven traffic jam? NCCL (NVIDIA Collective Communications Library) answers that by providing collective operations—broadcast, reduce, all-reduce, gather, and more—optimized for GPU-to-GPU communication.

Foundations of NCCL Collectives

A collective operation involves a group of ranks (one per GPU process) that participate in the same call sequence. The key engineering rule is simple: every rank must call the same collective with matching parameters (same communicator, same operation type, same tensor sizes). If one rank diverges, you get hangs that are painful to debug.

NCCL also assumes you are working with device pointers. That means you typically stage data on the GPU, then call NCCL so it can move bytes directly between devices.

A practical mental model is “one collective call, many coordinated transfers.” For example, all-reduce sums a tensor across all ranks and returns the result to every rank. That pattern shows up in distributed gradient aggregation and in global norm computations.

Communicator Setup and Rank Discipline

Before any collective, you create an NCCL communicator that defines the participating ranks and their mapping to devices. In engineering terms, this is where you prevent accidental cross-wiring.

A typical flow is:

1. Choose a rank per process and set the active CUDA device.
2. Create an NCCL unique ID on one process and broadcast it to all ranks.
3. Initialize an ncclComm_t per rank using that ID and the rank index.
4. Run collectives in a consistent order across ranks.

The “consistent order” part matters because NCCL uses internal state. If your code conditionally calls an all-reduce only on some ranks, you will eventually pay for it.

Choosing the Right Collective

Common collectives map cleanly to scientific and engineering tasks:

• All-Reduce: sum or max across ranks, result on all ranks. Use for global reductions and gradient-like aggregations.
• Reduce: combine to a single root rank. Use when only one GPU needs the final value.
• Broadcast: send one rank’s data to all ranks. Use for distributing parameters or initial conditions.
• All-Gather: collect distinct tensors from all ranks and distribute the concatenation to all ranks. Use for assembling distributed vectors.

For performance, prefer collectives that match your dataflow. If you need a global sum everywhere, all-reduce beats reduce+extra broadcast.

Stream Integration and Overlap

NCCL operations are asynchronous with respect to the host when launched on a CUDA stream. The engineering trick is to place NCCL calls on the same stream that produces the input data, or to use explicit stream synchronization so NCCL never reads incomplete buffers.

A common pattern is:

• Launch a kernel to compute a tensor on each GPU.
• Call NCCL on that tensor in the same stream.
• Launch the next kernel that consumes the reduced result.

This keeps dependencies on the GPU timeline rather than bouncing through the CPU.

Example All-Reduce for a Global Sum

Below is a minimal sketch of an all-reduce on device memory. The code assumes you already created comm, set cudaSetDevice, and have sendbuf and recvbuf allocated on the GPU.

// Assume: comm is initialized, sendbuf/recvbuf are device pointers.
// Assume: stream is a valid CUDA stream.

size_t count = N; // number of elements
ncclDataType_t dtype = ncclFloat32;

// Sum across ranks, result in recvbuf on every rank
ncclResult_t r = ncclAllReduce(
    sendbuf,
    recvbuf,
    count,
    dtype,
    ncclSum,
    comm,
    stream
);

if (r != ncclSuccess) {
  // Handle error
}

// Later: launch kernels that use recvbuf

If you want a global maximum instead, switch ncclSum to ncclMax. If your tensor is half precision, choose the matching NCCL datatype and ensure your computation semantics are acceptable.

Engineering Checklist for Collective Calls

Before you trust a collective in a tight iteration loop, verify these points:

• Every rank executes the same collective calls in the same order.
• The element count and datatype match exactly across ranks.
• The input buffer is fully produced on the GPU before NCCL reads it.
• The output buffer is consumed only after NCCL completes on the relevant stream.

Once those are in place, NCCL becomes a reliable building block: you write the math once, and the communication pattern stays consistent across GPUs without you manually orchestrating each transfer.

9.4 Implementing Domain Decomposition with Halo Exchange

Domain decomposition splits a large 3D or 2D grid across multiple GPUs. Each GPU owns a subdomain and also keeps a thin “halo” region around its boundary so it can compute stencil-like updates without waiting for every neighbor cell. The key engineering goal is to overlap communication with computation while keeping indexing and correctness airtight.

Core Concepts and Data Ownership

Start by defining a global grid of size \(N_x \times N_y \times N_z\). Choose a partition direction, commonly \(x\) for simplicity: GPU \(g\) owns \([x_{g0}, x_{g1})\). The halo thickness equals the stencil radius \(r\). For a 7-point stencil, \(r=1\); for a 27-point stencil with larger reach, \(r\) may be 2.

Each GPU stores its local array with extra halo layers. If the owned region is \(n_x\) in x, the allocated storage is \(n_x + 2r\). The interior region maps to indices \([r, r+n_x)\) in the local array. This convention makes halo updates a pure copy into the halo slots.

Halo Exchange Plan

For each time step (or iteration), you need neighbor boundary values before computing updates that depend on them.

1. Pack or directly copy halo faces from the interior boundary of each GPU into send buffers.
2. Exchange halos with neighbors using peer-to-peer copies or explicit device-to-device copies.
3. Compute interior first, using data that does not depend on newly received halos.
4. Compute boundary next, after halo data arrives.

For an \(x\)-partition, each GPU exchanges two faces: left neighbor and right neighbor. If a GPU is at the global boundary, it either uses fixed boundary conditions or mirrors values depending on your physics.

Example: 1D Halo Exchange for a 3D Stencil

Assume you partition along \(x\) and update a field \(u\) using neighbors in \(x\), \(y\), and \(z\). The halo exchange only concerns \(x\) faces; \(y\) and \(z\) are local to each GPU.

Let \(r=1\). GPU \(g\) owns x indices \([x_{g0}, x_{g1})\) globally. In local storage, interior x runs from local index 1 to local index \(n_x\). The halo slots are local x=0 (left halo) and local x=\(n_x+1\) (right halo).

• Send left face: local x=1 values to neighbor’s right halo.
• Send right face: local x=\(n_x\) values to neighbor’s left halo.

A practical way to keep copies simple is to store the array in a layout where the x-face is contiguous. For example, use a structure like \(u[x][y][z]\) with x as the fastest-changing index in memory, or flatten indices so that the face stride is 1.

Example: Streamed Overlap with Interior and Boundary Kernels

Use two CUDA streams per GPU: one for halo exchange copies and one for computation. The pattern is:

• Launch halo exchange copies in the comm stream.
• Launch an interior kernel in the compute stream that only reads interior neighbors. For \(r=1\), the interior kernel can skip x indices 1 and \(n_x\) because those depend on halos.
• Synchronize on the comm stream completion.
• Launch a boundary kernel that computes x indices 1 and \(n_x\) using the updated halo.

This ordering avoids a full device-wide sync and keeps the GPU busy while data moves.

Indexing Rules That Prevent Off-by-One Bugs

Use a single mapping function conceptually:

• Global x index \(x\) maps to local \(x_{loc} = x - x_{g0} + r\).
• Interior region satisfies \(r \le x_{loc} < r+n_x\).
• Left halo corresponds to \(x_{loc}=0\), right halo to \(x_{loc}=r+n_x\).

When packing faces, always pack from the interior boundary indices that match what the neighbor expects. For \(r=1\), that means pack from \(x_{loc}=r\) and \(x_{loc}=r+n_x-1\) if you use zero-based interior indexing in your kernel.

Synchronization and Correctness Checks

Correctness hinges on two things: halos must be updated before boundary computations, and boundary conditions must be applied consistently.

• For internal GPUs, halos come from neighbors.
• For edge GPUs, halos come from boundary conditions. For example, Dirichlet boundaries set halo values to a constant; Neumann boundaries mirror the nearest interior value.

A simple sanity check is to run with two GPUs and compare against a single-GPU run on the same global grid. If the halo exchange and indexing are correct, the results match bitwise for deterministic operations.

Practical Engineering Notes

Minimize packing overhead by copying directly from the array if the face is contiguous. If it is not contiguous, consider a lightweight pack kernel that writes into a contiguous send buffer; then the exchange becomes a single contiguous copy per face. Either way, keep the halo thickness equal to the stencil radius so you never “almost” have the data you need.

9.5 Coordinating Streams and Events Across Multiple GPUs

When you scale to multiple GPUs, the hard part is rarely “making kernels run.” It’s making the whole timeline line up: transfers, halo exchange, and compute must overlap without reading data that hasn’t arrived yet. Streams and events give you that control, but only if you treat them as a dependency graph rather than a pile of queues.

Core Idea: Treat Streams as Lanes

Each GPU has its own default stream and additional streams. Operations issued to the same stream execute in order; operations in different streams can overlap. To coordinate across streams, you record an event in one stream and make another stream wait on it. On multi-GPU systems, you also need to decide where the dependency lives: on the producing GPU, on the consuming GPU, or both.

A practical rule: use one stream for each “stage” on each GPU (compute, H2D/D2H, halo send, halo receive). Then connect stages with events so the consumer stream starts only after the producer has completed the relevant copy.

Step 1: Choose a Communication Pattern

For domain decomposition with halo exchange, each GPU typically needs to:

1. compute interior cells that do not depend on halo data,
2. exchange boundary data with neighbors,
3. compute boundary-adjacent cells once halos arrive.

That maps cleanly to two compute streams per GPU: one for interior and one for boundary. Transfers for halo send and receive can use dedicated copy streams.

Step 2: Use Events to Encode Dependencies

On each GPU, create events for:

• “halo receive complete”
• “halo send complete” (often optional if you only need the receive side to proceed)
• “interior compute complete” if boundary compute depends on it

Record events after the relevant operation in the producing stream. Then call cudaStreamWaitEvent in the consuming stream.

Step 3: Keep Copies and Compute Overlapped

A common timeline looks like this:

• GPU k starts interior compute in compute_interior[k].
• In parallel, it starts halo receive in copy_recv[k] and halo send in copy_send[k].
• When halo receive finishes, compute_boundary[k] waits on the “halo receive complete” event.
• If boundary compute also needs interior results (for example, shared buffers), it waits on “interior compute complete” too.

This avoids a global barrier that would stall every GPU.

Example: Two GPUs with Halo Exchange and Overlapped Compute

Assume GPU 0 sends its right halo to GPU 1, and GPU 1 sends its left halo to GPU 0. Each GPU has two compute streams and two copy streams.

// Pseudocode sketch for GPU 0 and GPU 1
// Streams: ci[k], cb[k], cs[k], cr[k]
// Events: eHaloRecv[k], eInteriorDone[k]

// GPU 0
cudaMemcpyAsync(sendBuf0, hostOrPeerSrc0, bytes, cudaMemcpyDeviceToDevice, cs[0]);
cudaEventRecord(eHaloRecv[0], cr[0]); // record after halo receive copy completes
launchInteriorKernel<<<... , ci[0]>>>(...);
cudaEventRecord(eInteriorDone[0], ci[0]);

cudaStreamWaitEvent(cb[0], eHaloRecv[0], 0);
cudaStreamWaitEvent(cb[0], eInteriorDone[0], 0);
launchBoundaryKernel<<<... , cb[0]>>>(...);

// GPU 1 mirrors the pattern with its own neighbor direction

The key detail is that boundary compute does not start when the program reaches it; it starts when the events it waits on become true.

Step 4: Decide Where to Synchronize

You only need a device-wide synchronization when you need final results on the host or when a later phase truly depends on everything. For phase boundaries, prefer:

• stream-level synchronization for the streams that produce the required outputs,
• event-based synchronization for internal dependencies.

A clean pattern is: after launching all work for a phase, synchronize only the streams that write the final buffers for that phase, then proceed.

Step 5: Validate with Timeline Reasoning

Before you chase micro-optimizations, confirm that overlap actually happens. If boundary compute starts late, check whether:

• the event was recorded in the correct stream,
• the wait was issued to the correct consuming stream,
• the halo copy direction matches the neighbor mapping,
• the buffers used by the copy are not being overwritten early.

A good sanity check is to temporarily insert a deliberate delay in the copy stream and verify that only the boundary compute shifts, not the interior compute.

10.1 Scaling Stencil Solvers with Overlap of Compute and Exchange

Stencil solvers spend time in two places: updating grid points (compute) and moving boundary data between partitions (exchange). Overlap means you start exchanging halo regions early, then keep the GPU busy computing interior points while communication is in flight. The trick is to structure both the kernel work and the host-side scheduling so neither side waits unnecessarily.

Core Idea and Work Decomposition

Assume a 3D 7-point stencil on a domain split across GPUs by subdomains with one-cell halos. Each iteration needs:

• Interior update: points whose stencil footprint stays entirely inside the local subdomain.
• Halo update: points adjacent to partition boundaries that depend on neighbor halo values.

To overlap effectively, you pipeline these steps:

1. Launch a kernel for interior points.
2. Pack and start halo exchange for the next iteration’s boundary values.
3. Launch a kernel for halo-adjacent points after halos arrive.

This requires double buffering for grid arrays (e.g., u and u_next) so that interior computation for iteration k can proceed while halos for iteration k+1 are being transferred.

Step-by-Step Execution Strategy

1) Partition Geometry and Halo Packing

Define your local grid dimensions as nx, ny, nz excluding halos, and allocate arrays with halo padding like (nx+2) x (ny+2) x (nz+2). Halo packing extracts contiguous planes when possible. For example, for the x-direction neighbors:

• Send left halo: u[1, :, :] (first interior plane)
• Send right halo: u[nx, :, :] (last interior plane)
• Receive into: u[0, :, :] and u[nx+1, :, :]

If your data layout makes these planes non-contiguous, consider reorganizing to keep halo regions contiguous so packing can be a single async copy rather than a scatter kernel.

2) Stream and Event Layout

Use at least two CUDA streams per GPU:

• Compute stream: runs interior and halo kernels.
• Transfer stream: runs async device-to-host (or device-to-device) copies and any packing kernels.

Then use events to connect them. The compute stream should not wait for halo exchange before interior work starts. The halo kernel should wait on a “halo received” event.

3) Interior Kernel Launch

Launch an interior kernel that updates indices that do not depend on halos. For a one-cell halo stencil, interior indices are:

• i = 2..nx-1, j = 2..ny-1, k = 2..nz-1 (adjust for your exact indexing)

This kernel reads only local u values and writes u_next. It can run immediately after the previous iteration’s buffer swap.

4) Start Halo Exchange Early

While the interior kernel runs, pack halo planes for the next iteration and start communication.

A common pattern is:

• Record an event after interior kernel begins (or after any writes to halo send buffers are complete).
• On the transfer stream, pack halo send buffers from u into contiguous buffers.
• Start non-blocking sends/receives.

If you use GPU-aware MPI, you can often send directly from device buffers. If not, use pinned host buffers and async copies.

5) Halo Kernel After Receive Completion

When halo receives complete, signal an event that the halo data is ready. Then launch the halo-adjacent kernel that updates points whose stencil footprint touches halos.

For the same one-cell halo stencil, halo-adjacent indices include planes at i=1 and i=nx, and similarly for j and k.

6) Swap Buffers and Repeat

After halo kernel finishes, swap u and u_next and proceed to the next iteration. The swap is where correctness is enforced: interior and halo kernels must write to the same destination buffer for that iteration.

Example: Minimal Pipelined Loop Skeleton

for (int iter = 0; iter < iters; ++iter) {
  // 1) Interior compute
  launch_interior_kernel<<<grid, block, 0, computeStream>>>(u, u_next);

  // 2) Pack halos and start exchange
  pack_halos<<<packGrid, packBlock, 0, transferStream>>>(u, sendBuf);
  mpi_irecv(recvBuf, ...);
  mpi_isend(sendBuf, ...);

  // 3) Wait for halos then compute boundary-adjacent
  wait_for_mpi_receives_and_record_event(haloReadyEvent);
  cudaStreamWaitEvent(computeStream, haloReadyEvent, 0);
  launch_halo_kernel<<<grid, block, 0, computeStream>>>(u, u_next);

  // 4) Swap
  std::swap(u, u_next);
}

The key is that interior compute and halo exchange overlap in time. If your halo packing is slow or forces synchronization, the overlap shrinks and scaling suffers.

Correctness Checklist That Prevents “It Works on My Machine”

• No halo reads in interior kernel: interior indices must exclude any point whose stencil touches halo cells.
• Send buffers are stable: don’t overwrite u or reuse send buffers until the send completes.
• Halo kernel waits: boundary-adjacent kernel must not start until receives are complete and halo data is written.
• Consistent iteration semantics: halos exchanged must correspond to the same source buffer used by the next iteration’s kernels.

When these rules hold, overlap becomes a predictable performance lever rather than a lucky accident.

10.2 Scaling Reductions and Global Norm Computations

Reductions are where performance goes to be measured, not guessed. A “global norm” is a reduction over many elements, often followed by a scalar operation and sometimes a collective communication across GPUs. Scaling it means keeping three things under control: (1) how much work each GPU does, (2) how efficiently each GPU reduces locally, and (3) how cheaply GPUs combine their partial results.

Local Reduction Foundations

Start with the simplest decomposition: each GPU computes a partial sum (or partial sum of squares) for its assigned data, then the host or a device collective combines those partials into the final scalar.

A typical global L2 norm looks like this conceptually:

1. Each GPU computes partial = Σ x_i^2 over its chunk.
2. GPUs combine partials into total = Σ partial.
3. The norm is sqrt(total).

The local kernel should reduce in two stages: intra-block reduction to produce one value per block, then a second kernel (or the same kernel with grid-stride and atomics, depending on constraints) to reduce block results. The key engineering choice is to avoid a long dependency chain where every thread waits on every other thread.

Choosing a Reduction Strategy That Scales

For dense reductions, a common pattern is:

• Use a block-level reduction in shared memory.
• Reduce within a warp using warp-level primitives to reduce synchronization overhead.
• Write one partial per block to global memory.
• Reduce the partials until a single value remains.

This structure scales because the expensive synchronization is limited to within a block, and the number of blocks shrinks quickly across stages.

Example Local Kernel Structure

Below is a compact sketch of a two-stage approach. It assumes the input is already on the device and uses a grid-stride loop to cover arbitrary sizes.

__global__ void partialSquares(const float* x, float* blockOut, int n) {
  extern __shared__ float s[];
  int tid = threadIdx.x;
  int gid = blockIdx.x * blockDim.x + tid;

  float sum = 0.0f;
  for (int i = gid; i < n; i += gridDim.x * blockDim.x) {
    float v = x[i];
    sum += v * v;
  }

  s[tid] = sum;
  __syncthreads();

  for (int stride = blockDim.x / 2; stride > 32; stride >>= 1) {
    if (tid < stride) s[tid] += s[tid + stride];
    __syncthreads();
  }

  // Final warp reduction without __syncthreads
  if (tid < 32) {
    // assume warpReduceSum is available
    float wsum = warpReduceSum(s[tid]);
    if (tid == 0) blockOut[blockIdx.x] = wsum;
  }
}

A second kernel reduces blockOut to a single scalar. If you keep the partial array small, the second stage is cheap and predictable.

Global Combination Across GPUs

Once each GPU has its partial sum, the global norm is just a collective reduction of scalars. The engineering goal is to minimize overhead and avoid unnecessary host-device round trips.

A practical workflow is:

• Partition the vector across GPUs with equal-sized chunks when possible.
• Launch the local partial reduction on each GPU in its own stream.
• Synchronize only where needed to ensure partial scalars are ready.
• Perform an all-reduce (sum) across devices for the partial scalars.
• Compute sqrt on the resulting scalar.

If you use NCCL for the all-reduce, the communication cost is mostly a function of the number of GPUs and the scalar size (one float or double). That means the local reduction quality matters more than the communication payload.

Numerical and Validation Details

Global norms are sensitive to accumulation order. Even if the math is “the same,” floating-point associativity means results can differ slightly across GPU counts. For engineering sanity, validate with a tolerance and compare against a CPU reference using the same precision.

If you need tighter agreement, consider accumulating in double for the sum of squares even if the input is float. That increases arithmetic cost, but it often pays off in predictable error bounds.

Case-Driven Example: L2 Norm Across Multiple GPUs

Suppose you have n = 10^8 floats and g = 4 GPUs. Split the input into four contiguous chunks of size about n/4. Each GPU runs partialSquares to produce blockOut values, then reduces those to one partialSum[g]. Finally, all-reduce sums partialSum[0..3] into totalSum, and compute sqrt(totalSum).

The scaling behavior is straightforward: if local reduction is efficient, the all-reduce of a single scalar is rarely the bottleneck. If local reduction is inefficient (too many stages, excessive synchronization, or poor memory access), then scaling stalls because each GPU spends more time producing its partial.

The engineering takeaway is simple: treat the global norm as two problems—local reduction quality and global scalar combination—and optimize them separately without mixing concerns.

10.3 Scaling Sparse Matrix Operations with Partitioned Storage

Sparse matrix work scales well only when the data layout matches the computation. Partitioned storage is the practical bridge: you split the matrix into blocks that fit your execution and communication patterns, then store each block in a sparse format that minimizes wasted work.

Core Idea of Partitioned Storage

Start with a sparse matrix A and a vector x. The baseline operation is y = A x. The scaling problem appears when you distribute rows across devices: each device needs the x entries referenced by its local nonzeros. Partitioned storage tackles this by aligning three things: (1) row ownership, (2) where nonzeros live, and (3) which x segments must be available locally.

A common approach is row partitioning. Device k owns a contiguous range of rows. For each owned row, you store its nonzeros locally. Then you build an index set of required x elements, often called a gather list. This list drives communication and also guides how you pack x for faster device access.

Choosing a Partition Strategy

Row partitioning is usually the simplest because it keeps y ownership local. Two practical variants matter for engineering:

1. Contiguous row blocks: easy to implement, stable mapping, good for predictable memory access.
2. Graph-aware row blocks: reduce the size of gather lists by grouping rows that reference similar x indices.

Even if you start with contiguous blocks, you can still apply the same engineering steps: compute gather lists, pack x, and run a local sparse kernel.

Partitioned Storage Layout

Within each device, store the local submatrix A_k using a sparse format that matches your sparsity structure.

• CSR for irregular row lengths: each row has its own nonzero count, and CSR stores row pointers plus column indices.
• ELL or HYB for structured sparsity: if many rows share similar nonzero counts, ELL reduces index overhead.

For scaling, the key is that the local format should support fast traversal of nonzeros while keeping indices small and predictable.

A practical engineering trick is to remap global column indices to local x offsets. After you build the gather list, you create a mapping col_global -> col_local. Then your local sparse kernel reads x_local[col_local] instead of chasing global indices.

Data Flow and Execution Steps

A systematic pipeline for device k looks like this:

1. Partition rows: determine row range [r0, r1) for device k.
2. Extract local nonzeros: build A_k from A restricted to those rows.
3. Build gather list: collect all global column indices referenced by A_k.
4. Pack x: create x_local by gathering x_global[gather[i]].
5. Remap indices: convert each nonzero’s column index to its position in gather.
6. Compute: run y_local = A_k x_local.
7. Scatter y: place y_local into the global y range owned by device k.

This pipeline keeps the sparse kernel focused on arithmetic and memory reads, while communication and indexing happen in preparation steps.

Example: Two-Device Row Partition with CSR

Assume A is stored in CSR globally. You split rows into two ranges: device 0 owns rows 0..R-1, device 1 owns rows R..N-1.

For device 0:

• Extract CSR arrays for rows 0..R-1 to form CSR\( A_0 \).
• While scanning A_0’s column indices, build gather_0 = unique(sorted(columns referenced)).
• Pack x_local_0 by reading x_global at gather_0.
• Remap each column index in A_0 to its position in gather_0.

Now the local SpMV kernel can do:

• For each local row i: sum over local nonzeros j in that row
• Read x_local_0[col_local(j)]
• Write y_local_0[i]

Device 1 repeats the same steps independently. After both devices finish, each writes its y_local into the global y range it owns. No device needs to write the other’s rows, so synchronization is straightforward.

Engineering Details That Prevent Slowdowns

• Unique gather lists: duplicates inflate x_local size and waste bandwidth. Use a uniqueness step when building gather.
• Remapping cost vs. kernel speed: remapping adds overhead, but it often pays off because the kernel avoids global index lookups.
• Load balance by nonzeros: equal row counts can still produce uneven work. Partition by nonzero counts when possible.
• Index width: if gather lists are small, store local column indices in narrower integer types to reduce memory traffic.

Case Study: When Partitioning Helps and When It Doesn’t

If your matrix is nearly dense in the sense that each row references most columns, gather lists become large. Partitioning still works, but communication and x packing dominate, and the sparse kernel may spend more time waiting on data than computing. In contrast, if each row references a limited neighborhood of columns, gather lists stay compact, remapping reduces index overhead, and scaling improves because each device can compute with mostly local data.

The engineering takeaway is simple: partitioning is effective when it reduces the size and fragmentation of the x data each device must fetch, not just when it splits rows evenly.

10.4 Scaling Batched Workloads with Device Local Scheduling

Batched workloads show up when you have many independent problems that share the same kernel logic: many right-hand sides in a linear solve, many small simulations, many short time windows, or many parameter sweeps. Scaling across multiple GPUs is easiest when each batch item can run mostly independently. The trick is to schedule those items so each GPU stays busy while minimizing cross-device traffic.

Core Idea: Keep Work Local

Device local scheduling means each GPU pulls work from its own queue and only exchanges data when it must. You partition the batch into chunks per GPU, then within each GPU you schedule the chunk using streams so transfers and compute overlap.

A practical mental model: each GPU has three lanes—H2D transfers, kernel compute, and D2D or D2H transfers for results. If you feed the compute lane with a steady stream of ready batch items, the GPU stops waiting. If you also overlap transfers for the next items, you hide latency behind useful work.

Work Partitioning Strategies

Start with a partitioning rule that matches your data layout.

1. Contiguous batch slices: GPU 0 gets items [0..B0), GPU 1 gets [B0..B1), etc. This is simplest when your batch dimension is stored contiguously.
2. Strided assignment: GPU g gets items g, g+N, g+2N.... This can help when some items have different sizes and you want a rough load balance.
3. Two-level partitioning: first split by size class (small vs large items), then split each class across GPUs. This avoids one GPU getting stuck with all the heavy items.

A useful engineering check: measure per-item runtime variance on one GPU. If variance is high, prefer strided or two-level partitioning.

Device Local Scheduling with Streams

Within a GPU, you typically use multiple streams to pipeline batch items. Each stream handles one “stage” of the pipeline for a sequence of items.

A clean pattern is double buffering per stream: while stream 0 computes item k, stream 1 transfers item k+1 and stream 2 can transfer results for item k-1.

Example: Pipelined Batched Kernel Launch

Assume each batch item needs an input vector and produces an output vector. You allocate per-GPU pinned host buffers and per-GPU device buffers for two buffers per stream.

// Pseudocode sketch for one GPU
for (int s = 0; s < numStreams; ++s) {
  cudaStreamCreate(&streams[s]);
}
for (int i = start; i < end; ++i) {
  int buf = i % 2;
  int s = i % numStreams;

  // 1) Copy input for item i
  cudaMemcpyAsync(d_in[buf], h_in[i], bytes,
                  cudaMemcpyHostToDevice, streams[s]);

  // 2) Compute
  myKernel<<<grid, block, 0, streams[s]>>>(d_in[buf], d_out[buf]);

  // 3) Copy output back
  cudaMemcpyAsync(h_out[i], d_out[buf], bytes,
                  cudaMemcpyDeviceToHost, streams[s]);
}
for (int s = 0; s < numStreams; ++s) cudaStreamSynchronize(streams[s]);

The key is that buf prevents overwriting buffers still in flight, and s spreads work across streams so the GPU can overlap operations.

Handling Dependencies Without Breaking Locality

Sometimes batch items share read-only data, like a common matrix or geometry. If that data is identical across GPUs, you can replicate it once per GPU and keep the rest local. If batch items depend on each other, you must introduce ordering, but you can still keep most data local by exchanging only the minimal boundary state.

A common compromise: replicate read-only inputs, and only transfer small “control” results that determine the next steps. This keeps bandwidth usage predictable.

Example: Load Balancing with Strided Assignment

Suppose you have 10,000 batch items and each item runtime varies because some inputs trigger extra iterations. If you split contiguously, GPU 1 might get a cluster of “hard” items.

Use strided assignment: GPU g runs items g, g+N, g+2N... where N is the number of GPUs. This spreads hard items across devices without any runtime prediction. Then apply the same per-GPU stream pipeline to its assigned items.

Practical Tuning Checklist

• Choose buffer count: start with double buffering; increase only if you see stalls from overwriting or insufficient overlap.
• Tune stream count: too few streams underutilize overlap; too many can increase overhead and reduce cache locality.
• Align batch layout: ensure each batch item’s input and output are contiguous so transfers are large and efficient.
• Validate correctness per stream: if kernels write to shared device buffers, you must isolate buffers per in-flight item.

Device local scheduling is not magic; it’s disciplined bookkeeping. Partition the batch so each GPU gets a fair share, pipeline transfers and compute with streams, and keep data movement intentional. The result is scaling that looks boring—in the best way.

10.5 Measuring Scaling Efficiency with Controlled Experiments

Scaling efficiency answers a simple question: when you add more GPUs, how much of the extra work actually turns into extra throughput? The trick is to measure it under controlled conditions so you don’t accidentally reward faster hardware, luckier scheduling, or a different batch composition.

Core Definitions and What You Measure

Start by defining the quantities you will compute every time.

• Strong scaling: fixed total problem size, increasing number of GPUs. Efficiency should stay near 1 if the added GPUs help.
• Weak scaling: problem size grows with GPU count so each GPU gets the same amount of work. Efficiency is about whether overhead grows slower than the work.

Use a consistent time window. For iterative solvers, measure from the start of the first iteration to the end of the last iteration, excluding one-time setup like memory allocation. For multi-GPU runs, include communication that occurs during the measured window, not just the kernel time.

Compute:

• Speedup: \(S(N)=T(1)/T(N)\)
• Parallel efficiency: \(E(N)=S(N)/N\)

If \(E(N)\) drops, it’s not automatically “bad scaling.” It means overheads are consuming more of the total time.

Controlled Experiment Design

A controlled experiment keeps everything constant except the scaling variable.

1. Fix the workload: same input data distribution, same stopping criteria, same number of iterations (or same residual threshold).
2. Fix the batch composition: if your solver uses batches of right-hand sides, keep batch size and ordering identical.
3. Fix the GPU mapping: pin processes to specific GPUs and keep the same device order across runs.
4. Fix the communication pattern: same halo width, same partitioning scheme, same collective operations.
5. Warm up: run a few iterations before timing to stabilize clocks, caches, and memory behavior.

A practical rule: if you can’t explain why two runs differ besides GPU count, you don’t have a controlled experiment.

Timing Strategy That Doesn’t Lie

Use two layers of timing.

• Wall-clock time: what the user experiences. Measure with a host timer around the synchronized region.
• Component timing: what the system spends time on. Use GPU events for kernel and device-to-device transfers, and host timers for orchestration overhead.

Synchronize carefully. For example, record a start event after the last dependency is enqueued, and record an end event after the last dependency is enqueued, then synchronize on the end event before reading elapsed time.

Example: Strong Scaling for a Halo Exchange Stencil

Suppose you solve a 3D stencil with domain decomposition. Each GPU updates its subdomain and exchanges halo planes.

Controlled setup:

• Same global grid size for all runs.
• Same halo width and same number of iterations.
• Same partitioning rule so each GPU gets the same number of interior cells.

Measurement plan:

• Time the full iteration loop.
• Record component times per iteration: compute kernels, halo pack/unpack, and halo exchange.

Interpretation:

• If compute time per iteration drops roughly like \(1/N\) but halo exchange stays flat, efficiency will fall at larger \(N\).
• If compute time doesn’t drop much, you likely have load imbalance or synchronization stalls.

Example: Weak Scaling for Batched Reductions

Now consider a batched reduction where each GPU handles a fixed number of vectors. Weak scaling keeps per-GPU work constant.

Controlled setup:

• Increase global batch size proportionally with GPU count.
• Keep the reduction tree structure consistent.

Measurement plan:

• Time only the reduction phase and any required inter-GPU collectives.
• Track both total time and the fraction spent in collectives.

Interpretation:

• If efficiency stays high, overhead grows slowly compared to work.
• If efficiency drops, the collective portion is growing faster than the local compute.

Repeatability and Statistical Sanity

Run multiple trials per \(N\) and report the median time. If the spread is large, you’re measuring scheduling noise or contention rather than scaling behavior. Also compare component timings: a consistent shift from compute to communication across trials is meaningful; random swings usually aren’t.

A Simple Efficiency Checklist

Before trusting \(E(N)\):

• The workload and stopping criteria match.
• The timing window includes communication and excludes setup.
• GPU mapping and partitioning are consistent.
• Warm-up is done.
• Trials are repeated and summarized consistently.

When these conditions hold, efficiency becomes a reliable diagnostic tool rather than a mystery score.

11.1 Optimizing a 3D Stencil Kernel With Tiling and Shared Memory

A 3D stencil computes each grid point from a neighborhood, commonly a 7-point or 27-point pattern. The performance challenge is usually memory traffic: each output point needs several input values, and naive kernels repeatedly fetch overlapping regions from global memory. Tiling plus shared memory reduces those redundant fetches by loading a block’s “working slab” once, then reusing it for many output points.

Baseline Kernel: What You Start With

Assume a 3D grid nx * ny * nz stored in a flat array with idx = x + nx*(y + ny*z). A simple interior 7-point stencil might compute:

out[x,y,z] = c0*in[x,y,z] + c1*(in[x-1,y,z] + in[x+1,y,z] + in[x,y-1,z] + in[x,y+1,z] + in[x,y,z-1] + in[x,y,z+1])

In a naive kernel, each thread loads its center and all neighbors from global memory. Threads in the same block share neighbors, but the hardware cache may not capture enough reuse, especially when the working set is larger than cache capacity. The result is high global memory load count and low arithmetic intensity.

Tiling with Halo: The Core Idea

For a 7-point stencil, each output cell needs neighbors at ±1 in x, y, and z. If a thread block computes a tile of size (Bx, By, Bz) outputs, it also needs a halo of width 1 around that tile in all directions. That means shared memory must hold (Bx+2) * (By+2) * (Bz+2) input values.

A practical mapping is:

• Each thread corresponds to one output cell in the interior of the tile.
• Threads cooperatively load shared memory for both interior and halo cells.
• After a __syncthreads(), each thread reads only from shared memory to compute its output.

This turns many global loads into one-time shared loads per block.

Shared Memory Staging: Cooperative Loading

Shared memory is limited, so tile sizes must fit. For float values, shared usage is 4 * (Bx+2)*(By+2)*(Bz+2) bytes. For example, Bx=By=Bz=8 uses (10^3)*4 = 4000 bytes, which is comfortable. Larger tiles can improve reuse but may increase register pressure and reduce occupancy.

The loading pattern should be simple and regular. One approach is to have each thread load one or more shared elements using a linear index over the shared array. The key is to cover the entire shared tile, including halo, before computing.

Example: Tiled 3D Stencil Kernel Skeleton

__global__ void stencil3d_tiled(
  const float* __restrict__ in,
  float* __restrict__ out,
  int nx, int ny, int nz,
  float c0, float c1)
{
  const int Bx = blockDim.x, By = blockDim.y, Bz = blockDim.z;
  __shared__ float sh[(8+2)*(8+2)*(8+2)];

  int tx = threadIdx.x, ty = threadIdx.y, tz = threadIdx.z;
  int x0 = blockIdx.x * Bx;
  int y0 = blockIdx.y * By;
  int z0 = blockIdx.z * Bz;

  int sx = tx + 1, sy = ty + 1, sz = tz + 1;
  int gx = x0 + tx, gy = y0 + ty, gz = z0 + tz;

  // Load center into shared
  if (gx < nx && gy < ny && gz < nz)
    sh[(sx) + (Bx+2)*(sy + (By+2)*sz)] = in[gx + nx*(gy + ny*gz)];

  __syncthreads();

  // Interior-only compute example
  if (gx > 0 && gx < nx-1 && gy > 0 && gy < ny-1 && gz > 0 && gz < nz-1) {
    float center = sh[sx + (Bx+2)*(sy + (By+2)*sz)];
    float xm = sh[(sx-1) + (Bx+2)*(sy + (By+2)*sz)];
    float xp = sh[(sx+1) + (Bx+2)*(sy + (By+2)*sz)];
    float ym = sh[sx + (Bx+2)*((sy-1) + (By+2)*sz)];
    float yp = sh[sx + (Bx+2)*((sy+1) + (By+2)*sz)];
    float zm = sh[sx + (Bx+2)*(sy + (By+2)*(sz-1))];
    float zp = sh[sx + (Bx+2)*(sy + (By+2)*(sz+1))];
    out[gx + nx*(gy + ny*gz)] = c0*center + c1*(xm+xp+ym+yp+zm+zp);
  }
}

This skeleton shows the indexing discipline and the shared-only compute phase. In a complete implementation, you must load halo values into sh as well, not just the centers. A common method is to loop over shared indices and map them back to global coordinates with bounds checks.

Boundary Handling Without Sloppy Branching

Branching inside the hot compute loop can cost more than it saves. A clean approach is:

• Run a fast interior kernel that assumes x in [1, nx-2], y in [1, ny-2], z in [1, nz-2].
• Run a separate boundary kernel (or handle boundaries on the host) for the outer layers.

If you must keep one kernel, restrict conditionals to a single outer if that encloses the compute, and keep shared loading consistent.

Tuning Checklist That Actually Moves the Needle

1. Confirm shared reuse by checking that global load counts drop sharply versus the baseline.
2. Balance tile size and occupancy so shared memory and registers do not throttle active warps.
3. Keep indexing cheap by precomputing strides and using __restrict__ pointers.
4. Use fixed stencil coefficients to encourage compiler optimization and reduce redundant loads.
5. Validate numerics with a CPU reference on small grids, then scale up.

When tiling is correct, the kernel becomes limited by shared memory bandwidth and arithmetic rather than repeated global neighbor fetches. That shift is the whole point: fewer trips to global memory, more useful work per loaded value.

11.2 Optimizing a Reduction Kernel With Warp Level Primitives

Reduction kernels are the GPU equivalent of “sum everything, but do it fast and without losing correctness.” The baseline approach—one thread loads one element, then block-level reductions combine partial sums—works, but it often leaves performance on the table. Warp-level primitives help because they reduce communication overhead inside a warp, where threads execute in lockstep.

Foundational Setup: What We Reduce and How Threads Map

Assume we want a sum of N floats. A common mapping is grid-stride looping: each thread processes multiple elements separated by the total number of threads. This keeps the kernel simple and avoids requiring N to match the grid size.

Each thread accumulates a local sum in a register. That local accumulation is the first optimization: it reduces the number of values that must participate in the reduction tree.

Baseline Block Reduction with Shared Memory

A typical baseline reduces within a block using shared memory and a power-of-two stride. It usually looks like: write each thread’s partial to shared memory, synchronize, then halve the active range repeatedly.

This works, but it pays two costs repeatedly: (1) shared memory traffic and (2) __syncthreads() barriers. Barriers are especially expensive when the reduction stage is already confined to a warp.

Warp-Level Primitives: Reduce Without Shared Memory

Warp-level reduction uses shuffle operations to exchange values between lanes without shared memory. The key idea is that a warp’s lanes can read each other’s registers directly.

A systematic approach:

1. Perform a warp-level reduction on each lane’s register value.
2. Only one lane per warp (often lane 0) writes the warp’s result to shared memory.
3. Synchronize once per block.
4. Let the first warp reduce those warp results.

This reduces synchronization frequency and shrinks shared memory usage to “one value per warp,” not “one value per thread.”

Example: Warp Shuffle Reduction Kernel Skeleton

__inline__ __device__ float warp_reduce_sum(float v) {
    for (int offset = 16; offset > 0; offset >>= 1) {
        v += __shfl_down_sync(0xffffffff, v, offset);
    }
    return v;
}

__global__ void reduce_sum(const float* x, float* blockSums, int N) {
    int tid = threadIdx.x;
    int gid = blockIdx.x * blockDim.x + tid;

    float sum = 0.0f;
    for (int i = gid; i < N; i += gridDim.x * blockDim.x) {
        sum += x[i];
    }

    sum = warp_reduce_sum(sum);

    __shared__ float warpSums[32];
    int lane = tid & 31;
    int warpId = tid >> 5;

    if (lane == 0) warpSums[warpId] = sum;
    __syncthreads();

    if (warpId == 0) {
        float v = (lane < (blockDim.x + 31) / 32) ? warpSums[lane] : 0.0f;
        v = warp_reduce_sum(v);
        if (lane == 0) blockSums[blockIdx.x] = v;
    }
}

This kernel produces one partial sum per block. A second kernel (or a final reduction on the CPU) combines blockSums until a single value remains.

Correctness Details That Matter

1. Non-multiple sizes: The grid-stride loop naturally skips out-of-range indices.
2. Warp size assumptions: The shuffle loop assumes 32-lane warps. That matches current NVIDIA GPUs.
3. Floating-point order: Warp shuffles change the reduction order compared to a shared-memory tree. For strict reproducibility, you may need a deterministic reduction strategy, but for typical engineering workloads, the result is usually within expected floating-point error bounds.

Performance Reasoning: Why This Is Faster

Shared-memory reductions require multiple rounds of synchronization and reads/writes. The warp shuffle approach keeps most arithmetic in registers and removes barriers inside the warp. Shared memory is used only for warp leaders, so the shared memory footprint and traffic drop sharply.

A practical checklist:

• Ensure the input is accessed contiguously by threads in the first iteration to get coalesced loads.
• Keep block size reasonable so the number of warps fits the warpSums array.
• Watch register pressure; if local accumulation uses too many registers, occupancy can drop and performance can suffer.

Example: Verifying Against a CPU Reference

Compute a CPU sum in double precision and compare the GPU result with a tolerance. Use the same input data and ensure the GPU kernel is fully synchronized before reading results. If the difference is larger than expected, inspect reduction order sensitivity and confirm that the second-stage reduction is correct.

Summary of the Optimization Ladder

Start with grid-stride accumulation in registers, reduce within each warp using shuffle primitives, store only warp results to shared memory, then finish the block reduction using a single warp. This sequence systematically removes the expensive parts of the baseline reduction while keeping the kernel structure easy to reason about.

11.3 Optimizing Sparse Matrix Multiply With Format Selection

Sparse matrix multiply (SpMM) computes C = A × B, where A is sparse and B is dense. The core engineering problem is that “sparse” can mean many different structures, and the best storage format depends on how nonzeros are distributed and how you want to traverse them. Format selection is not a cosmetic choice; it changes memory access patterns, branch behavior, and the amount of useful work per memory transaction.

From Structure to Format Choice

Start by classifying A along two axes: row density and regularity.

• Row density: how many nonzeros per row. If most rows have similar counts, you can use formats that assume a predictable row structure.
• Regularity: whether nonzeros cluster in a few columns or are scattered.

A practical workflow is to compute quick statistics on the host: average nnz per row, variance of nnz per row, and the distribution of column indices per row. Then map those observations to a format.

Format Options and What They Trade

CSR (Compressed Sparse Row) stores row pointers plus column indices and values. It is a strong default when you need straightforward row-wise traversal.

• Strength: efficient row iteration and easy integration with row-based kernels.
• Cost: irregular column indices can cause scattered reads from B.

ELL (ELLPACK) stores a fixed number of entries per row, padding shorter rows. It works best when row lengths are similar.

• Strength: more regular memory access and fewer branches.
• Cost: padding overhead when row lengths vary.

COO stores explicit (row, col) triplets.

• Strength: simple construction.
• Cost: slower SpMM kernels because you typically need extra work to group by row.

A useful rule of thumb: if nnz per row has low variance, ELL often improves throughput by reducing control divergence. If nnz per row varies widely, CSR usually avoids padding waste.

Kernel Mapping That Makes Format Matter

For SpMM, each output element C[i, k] is a dot product over nonzeros in row i of A:

• For each row i, iterate over nonzeros (i, j) in A.
• Multiply A[i, j] × B[j, k] and accumulate.

A common GPU mapping is one thread block per row (or per group of rows), with threads covering the k dimension of B. This is where format selection pays off: the inner loop over nnz must be efficient, and the reads from B must have enough locality to avoid turning the kernel into a random-access benchmark.

Example: CSR SpMM with Row Blocks

Assume CSR arrays: row_ptr, col_idx, val. A CSR kernel typically does:

• block selects a row i from blockIdx.x
• each thread handles a subset of columns k in B
• loop over p from row_ptr[i] to row_ptr[i+1]
• load j = col_idx[p], a = val[p]
• accumulate a * B[j, k]

To reduce wasted work, guard the k range and keep the inner loop tight. Also, if B is stored so that B[j, k] for consecutive k is contiguous, threads in a warp will read contiguous addresses, which improves memory efficiency.

Example: ELL SpMM When Row Lengths Are Uniform

ELL stores ell_col[ row ][ t ] and ell_val[ row ][ t ] for t in [0, max_nnz_per_row). If you choose max_nnz_per_row based on a percentile (for example, the 90th percentile), you reduce padding while keeping rows mostly aligned.

In the ELL kernel:

• each thread block processes a row
• threads iterate t from 0 to max_nnz_per_row
• if an entry is padded, skip using a sentinel or a separate mask
• accumulate into C

The key benefit is that the loop trip count is uniform, which reduces divergence. The key risk is padding overhead: if many rows are much shorter than max_nnz_per_row, you spend time multiplying by zeros (or checking sentinels).

Choosing Between CSR and ELL with a Concrete Decision

Compute:

• mean_nnz = total_nnz / num_rows
• var_nnz or standard deviation
• max_nnz and a trimmed p90_nnz

Decision logic:

• If std_nnz / mean_nnz is small (rows are similar), set ELL’s max_nnz_per_row near p90_nnz and use ELL.
• If row lengths vary a lot, use CSR to avoid padding.

This decision is simple enough to automate, and it directly targets the two big costs: divergence (CSR) versus padding (ELL).

Practical Notes That Prevent “It Runs, but It’s Slow”

• Ensure B’s memory layout matches the thread mapping. If threads iterate over k, store B so k is contiguous.
• Keep index and value loads coalesced. CSR’s col_idx and val are contiguous within a row; mapping rows to blocks helps.
• Avoid extra indirection. COO often needs grouping into CSR/ELL before a fast SpMM kernel.

Case Study: Format Selection for a Mixed Sparsity Matrix

Suppose A has 100k rows. You measure:

• mean nnz per row = 12
• standard deviation = 2.5
• p90 nnz per row = 15

This indicates fairly uniform row lengths. ELL with max_nnz_per_row = 15 will have limited padding, and uniform loop counts reduce divergence. If instead standard deviation were 10 and p90 were 40, padding would dominate, and CSR would likely be faster because it only iterates over actual nonzeros.

In both cases, the winning format is the one that minimizes wasted work in the inner loop: either wasted iterations (ELL padding) or wasted control flow (CSR divergence).

11.4 Optimizing Data Layout for Multidimensional Transformations

Multidimensional transformations—like 2D/3D FFT-like stages, separable filters, and tensor contractions—often fail for the same reason: the computation is fast, but the data walks around like it has somewhere better to be. Optimizing data layout means arranging memory so that threads read and write contiguous (or at least predictable) regions, while keeping index math simple enough that it doesn’t become the hidden tax.

Core Idea: Map Indices to Memory Once

Start with a clear rule: choose a layout first, then write kernels that follow it. For an array A[x][y][z] stored in row-major order, the last index varies fastest in memory. That means threads that vary the last index should get coalesced accesses.

A practical workflow:

1. Pick the dimension that your kernel’s threads will iterate over most tightly.
2. Ensure that dimension is the fastest-changing one in memory.
3. For transforms that require different access patterns per stage, plan a layout change between stages (or use a transpose kernel).

Layout Choices for Multidimensional Arrays

Row-Major Versus Column-Major

In C/C++ with row-major storage, A[x][y][z] is contiguous in z. If your kernel uses threadIdx.x to cover z, each warp tends to hit consecutive addresses.

If instead threadIdx.x covers x, then each thread jumps by a large stride, and global memory transactions multiply.

Structure of Arrays Versus Array of Structures

When each element carries multiple fields (e.g., real and imaginary parts, or multiple channels), prefer Structure of Arrays (SoA):

• real[idx], imag[idx] rather than complex[idx] with interleaved fields.

SoA makes it easier to load only what a stage needs, and it improves coalescing because each array is a dense block.

Designing Kernels Around the Layout

Use Contiguous Thread Dimensions

For a 2D transform A[i][j], if you launch threads so that j changes fastest, then reads and writes align with contiguous memory. A common pattern is:

• threadIdx.x maps to j
• blockIdx.x maps to i

Then each block handles one row (or a tile of rows), and each warp reads a contiguous segment.

Tile for Shared Memory Staging

For operations that need neighbor values (stencils, local windows, separable filters), tile the input into shared memory. The tile shape should match the access pattern:

• If each thread reads A[i][j + k], then shared memory should store a contiguous span in j.

This reduces repeated global reads and turns many scattered accesses into shared-memory hits.

Transpose as a Layout Tool

Many separable transforms alternate between “row-friendly” and “column-friendly” stages. Instead of forcing every stage to fight the layout, insert a transpose between stages.

A transpose kernel is usually worth it when:

• The next stage has a fundamentally different access direction.
• The transpose cost is amortized over many operations per element.

Avoiding Bank Conflicts in Shared Memory

When transposing tiles, shared memory indexing can cause bank conflicts. A simple fix is to pad the shared tile width by 1 element so that columns don’t map to the same bank repeatedly.

Example: 2D Separable Transform with Two Layouts

Assume A is NxM in row-major order. Stage 1 processes rows, stage 2 processes columns.

1. Stage 1 kernel reads contiguous A[i][j].
2. Transpose B[j][i] = A[i][j].
3. Stage 2 kernel now treats B[j][i] rows as the original columns.

This keeps both stages coalesced without complicated strided loads.

Minimal Indexing Pattern

Let idx = i*M + j for A[i][j].

• Stage 1: threads vary j, so idx increments by 1 across a warp.
• Stage 2: after transpose: threads vary the new fast index, again producing contiguous accesses.

Example: 3D Transform and Dimension Reordering

For A[x][y][z] with row-major storage, z is contiguous. If your transform applies the same 1D operation along z first, you’re in luck: threads that vary z get coalesced reads.

If the next stage operates along x, you have two options:

• Reorder the tensor so that x becomes the fastest-changing dimension before that stage.
• Or transpose between stages to bring the working dimension into the contiguous position.

In both cases, the goal is identical: make the dimension your threads iterate over match the memory’s fastest-changing index.

Practical Checklist for Engineers

• Confirm which index is contiguous in memory for your chosen layout.
• Ensure the fastest-changing index is the one mapped to threadIdx.x (or the closest equivalent).
• Use SoA when kernels consume only subsets of fields.
• Insert transpose kernels when stage access patterns differ.
• Tile and pad shared memory to keep neighbor-heavy kernels efficient.

When these steps are followed, multidimensional transformations stop being a memory-access puzzle and start behaving like the arithmetic they were meant to perform.

11.5 End-to-End Optimization of an Iterative Solver Pipeline

An iterative solver pipeline is a loop with three recurring costs: (1) applying an operator, (2) reducing global quantities, and (3) updating vectors. Optimization works best when you treat the pipeline as a system rather than a collection of kernels.

Pipeline Map and Data Flow

The core idea is to keep data resident on the GPU, minimize synchronization points, and ensure each kernel has enough work to hide latency. A typical flow for a Krylov method looks like this:

• Start with vectors on device: x, b, r, p.
• Each iteration computes Ap = A*p.
• Compute scalars via reductions, then update x, r, and p.
• Stop when a norm criterion is met.

Step 1: Establish a Baseline That Matches Reality

Use a representative matrix and right-hand side size, not a toy case. Record three numbers per iteration: time spent in operator application, time spent in reductions, and time spent in vector updates. If reductions dominate, you will chase the wrong bottleneck.

A practical baseline checklist:

• Confirm all vectors are allocated on device once.
• Ensure the loop does not force device-to-host transfers for convergence.
• Use the same stopping criterion for both CPU and GPU runs.

Step 2: Optimize the Operator Application Kernel

The operator kernel usually determines the memory traffic. For a stencil, tiling with shared memory reduces redundant global loads. For sparse CSR, the main win comes from choosing a layout that improves locality and from using a kernel that matches the sparsity structure.

Concrete example: a CSR SpMV kernel often benefits from:

• Assigning one thread block per row chunk to reduce divergence.
• Using a consistent indexing strategy so reads of colIdx and val are predictable.
• Keeping x reads as efficient as possible, even though reuse is limited.

After changes, verify that the kernel still produces the same Ap within tolerance. A fast kernel that is wrong is just a faster way to fail.

Step 3: Reduce Scalars Without Turning the GPU into a Waiting Room

Dot products and norms are reduction-heavy and can introduce frequent synchronization. The goal is to reduce the number of global synchronization events and to keep intermediate results on the device.

A common pattern:

• Each block computes a partial sum in shared memory.
• A second stage reduces partial sums.
• The final scalar stays on device for subsequent updates.

To avoid repeated kernel launches, consider fusing reduction outputs into the update step when the data dependencies allow it. If fusion is not possible, at least ensure the reduction pipeline uses the same stream and does not force a host sync.

Step 4: Vector Updates with Minimal Temporaries

Vector updates are often bandwidth bound. The best improvements usually come from:

• Using a single kernel that performs multiple AXPY-like operations when possible.
• Avoiding extra intermediate arrays like tmp1, tmp2 when you can compute in registers.
• Ensuring contiguous memory access by using a structure-of-arrays layout for multiple right-hand sides.

Concrete example: if your iteration computes

• x = x + alpha * p
• r = r - alpha * Ap

you can implement both in one kernel by reading p and Ap once per element and writing x and r once per element.

Step 5: Convergence Check Without Host Bottlenecks

If you copy a scalar to the host every iteration, you pay a synchronization tax. Instead, keep the convergence scalar on the device and structure the loop so the host only learns the result when you must stop.

A simple engineering approach:

• Compute the norm on device.
• Store it in a device scalar.
• Use a device-side decision mechanism or batch the checks so the host sync frequency is reduced.

Step 6: Measure, Change One Thing, Re-Validate

Optimization is iterative. After each change, re-run the same iteration count and compare:

• Residual history shape (not just the final residual).
• Kernel times for operator, reductions, and updates.
• GPU memory throughput and achieved occupancy.

Example: One Iteration Timing Breakdown

Iteration 1
- Operator application: 42% of time
- Reductions: 35% of time
- Vector updates: 23% of time

Action
- If reductions stay high, focus on reduction staging and fusion.
- If operator application grows, re-check memory access and tiling.

Step 7: Keep Correctness Tight While Performance Moves

Use a CPU reference for a small problem and a GPU run for the same inputs. Check:

• Dot product agreement within tolerance.
• Residual norm monotonicity where expected.
• No indexing mistakes at boundaries.

When you change kernel fusion or reduction order, numerical differences can appear. That’s normal; what matters is that the solver still converges and the residual behavior remains consistent.

12.1 Project Structure With CMake and CUDA Compilation Targets

A maintainable CUDA project starts with a clean separation of concerns: host code, device code, and shared headers. CMake then becomes the “wiring diagram” that tells the compiler which files belong to which target, which flags apply, and how artifacts are organized. The goal is simple: when you add a new kernel file or a new executable, you should not edit five unrelated places.

Core Layout That Scales Past Hello World

A practical directory structure keeps responsibilities obvious:

• src/ for host and device implementation files
• include/ for headers shared across translation units
• kernels/ for CUDA .cu files when you want a clear boundary
• tests/ for correctness checks
• cmake/ for helper CMake modules

A common pattern is to compile CUDA code into a library target, then link that library into one or more executables. This avoids duplicating compilation flags and makes it easier to test kernels through a single interface.

Target Strategy That Keeps Flags Local

Use three kinds of targets:

1. A CUDA library for reusable kernels and device utilities.
2. One or more executables for applications and benchmarks.
3. Optional test executables that link the same library.

Each target should own its compile options. If you set -lineinfo or --use_fast_math, you want to know exactly which target got it.

CMake Essentials for CUDA Compilation

At minimum, you need:

• project(... LANGUAGES CXX CUDA) to enable CUDA language support
• add_library(... LANGUAGES CUDA) or add_library with .cu sources
• target_include_directories for include/
• target_compile_options for target-specific flags
• set_target_properties for CUDA-specific properties when needed

Example Mind Map

Mind maps are useful here because they show the flow from files to targets to build outputs.

A Minimal, Correct CMake Skeleton

This example shows a library that compiles CUDA sources and an executable that links it. Keep it small so the structure is the lesson.

cmake_minimum_required(VERSION 3.25)
project(cuda_engineering LANGUAGES CXX CUDA)

set(CMAKE_CXX_STANDARD 17)
set(CMAKE_CUDA_STANDARD 17)

add_library(cuda_kernels STATIC
  kernels/vector_add.cu
  src/device_utils.cu
)

target_include_directories(cuda_kernels PUBLIC include)

add_executable(vector_add_app src/main.cpp)

target_link_libraries(vector_add_app PRIVATE cuda_kernels)

Example: Per-Target Flags Without Global Chaos

If you set flags globally, you eventually regret it. Prefer target_compile_options.

target_compile_options(cuda_kernels PRIVATE
  $<$<CONFIG:Debug>:-G>
  $<$<CONFIG:RelWithDebInfo>:--generate-line-info>
)

target_compile_options(vector_add_app PRIVATE
  $<$<CONFIG:Debug>:-O0>
  $<$<CONFIG:RelWithDebInfo>:-O2>
)

Organizing Device Code Interfaces

A clean interface keeps host code from knowing too much about device internals. A typical approach:

• include/ declares host-callable functions like void launch_vector_add(...);
• kernels/ implements the kernels and the launch wrapper
• src/ contains host-side orchestration and validation

This separation matters because it reduces recompilation. When you change a kernel implementation, only the CUDA library rebuilds; the executable relinks.

Common Pitfalls and How the Structure Prevents Them

1. Mixing headers and compilation units unintentionally: headers should not contain heavy device code unless you intentionally use header-only device utilities.
2. Forgetting include visibility: if the executable cannot find headers, you’ll see it immediately. Use PUBLIC for library headers.
3. Overusing one giant target: a single target works for tiny projects, but it makes flag changes and incremental builds harder.
4. Debug flags applied everywhere: -G slows compilation and execution. Keep it scoped to the CUDA library target.

Mind Map: Build Flow

Practical Checklist for Adding a New Kernel File

When you add kernels/new_kernel.cu:

• Add the file to cuda_kernels sources.
• Keep the public launch wrapper declaration in include/.
• Ensure the executable calls only the wrapper, not the kernel directly.
• Run the existing test target or the smallest executable that exercises the new path.

This workflow keeps the build predictable and the codebase easy to navigate, even when the kernel count grows.

// Pseudocode for test logic
for (int N : {0, 1, 7, 256, 257, 1023, 1024, 1025}) {
  for (int block : {32, 64, 128, 256}) {
    int grid = (N + block - 1) / block;
    launch_kernel<<<grid, block>>>(d_in, d_out, N);
    cudaDeviceSynchronize();
    copy_back(out);
    ref = cpu_reference(in, N);
    assert_close(out, ref, tolerance);
  }
}

for (int N : {0, 1, 2, 3, 31, 32, 33, 1000}) {
  fill(in, pattern);
  launch_reduction<<<grid, block, shared>>>(d_in, d_partial, N);
  launch_final<<<1, final_block>>>(d_partial, d_out, N);
  copy_back(out);
  ref = cpu_reference_sum(in, N);
  assert_close(out, ref, 1e-6f);
}

runs/
  run_2026-04-15_1530/
    metadata.json
    timings.csv
    nsys_report.qdrep
    ncu_report.ncu-rep

// Pseudocode for repeatable timing
init_device(gpu_id);
set_deterministic_seed(seed);
prepare_inputs_once();
create_streams_consistently();

warmup(stream, warmup_iters);
for (i = 0; i < measure_iters; i++) {
  cudaEventRecord(start, stream);
  launch_kernels(stream, fixed_config);
  cudaEventRecord(stop, stream);
  cudaEventSynchronize(stop);
  times[i] = elapsed_ms(start, stop);
}

stats = compute_median_iqr(times);
write_metadata_and_csv(stats, run_id);
if (stats.regressed) {
  run_nsys_capture(run_id);
  run_ncu_capture(run_id, kernel_name);
}

// Device kernel
template<int TILE_X, int TILE_Y>
__global__ void stencil2d_kernel(
    const float* __restrict__ in,
    float* __restrict__ out,
    int nx, int ny,
    int inPitch, int outPitch)
{
    int x = blockIdx.x * TILE_X + threadIdx.x;
    int y = blockIdx.y * TILE_Y + threadIdx.y;
    if (x <= 0 || x >= nx-1 || y <= 0 || y >= ny-1) return;

    int idxIn  = y * inPitch  + x;
    int idxOut = y * outPitch + x;

    float c = in[idxIn];
    float l = in[idxIn - 1];
    float r = in[idxIn + 1];
    float u = in[idxIn - inPitch];
    float d = in[idxIn + inPitch];

    out[idxOut] = 0.25f * (l + r + u + d) + 0.0f * c;
}

struct Stencil2DParams {
    const float* in;
    float* out;
    int nx, ny;
    int inPitch, outPitch;
};

template<int TILE_X, int TILE_Y>
void launch_stencil2d(cudaStream_t s, const Stencil2DParams& p) {
    dim3 block(TILE_X, TILE_Y);
    dim3 grid((p.nx + TILE_X - 1) / TILE_X,
              (p.ny + TILE_Y - 1) / TILE_Y);
    stencil2d_kernel<TILE_X, TILE_Y><<<grid, block, 0, s>>>(
        p.in, p.out, p.nx, p.ny, p.inPitch, p.outPitch);
}

{
  "run_id": "2026-03-11T10:00:00Z",
  "gpu": "NVIDIA_XX",
  "cuda_version": "12.x",
  "commit": "abcdef123",
  "workload": {"nx": 256, "ny": 256, "nz": 128},
  "kernels": [
    {"name": "stencil3d",
     "block": [16,8,1],
     "grid": [16,32,128],
     "time_ms": {"mean": 2.12, "p95": 2.35},
     "dram_gbs": 910.4}
  ]
}