GPU Kernel Performance Engineering Lab

You change four lines of CUDA, recompile, and watch the kernel runtime drop by 40%. It feels like an instant win until your verification pass flags thousands of stale NaN values across the output matrix. The optimization didn't win anything; it just computed garbage at breakneck speed.

Real GPU performance engineering isn't about guessing clever tricks or chasing microbenchmarks in isolation. It's a disciplined, evidence-driven loop: form a concrete hypothesis about hardware bottlenecks, implement a single causal code modification, prove bit-exact correctness against edge-case matrices, and verify that Nsight Compute counters move for the exact mechanical reason you predicted.

In this lab, you'll optimize an out-of-place single-precision ($N \times N$) matrix transpose kernel. Transpose is the canonical proving ground for mastering GPU memory hierarchy because it performs zero arithmetic beyond indexing: exactly one read and one write per element. Every nanosecond spent in the kernel reflects global memory bus transactions, L1/L2 cache line utilization, shared memory bank conflicts, and instruction-level latency hiding.

Three prerequisites supply the vocabulary. CUDA for ML Training established grids, blocks, warps, asynchronous launches, and memory hierarchy. Accelerator Architecture Field Guide connected arithmetic intensity and resource ceilings to workload bounds. GPU Profiling, Correctness, and Benchmarking separated host timelines, kernel counters, and correctness evidence. Here those ideas control one kernel from first measurement through recovery.

Keep one experiment stable

The workload transposes an $N \times N$ row-major float32 matrix into a separate output allocation. This is an out-of-place transpose; passing overlapping input and output would violate the contract. If $A[r,c]$ lives at input[r * N + c], the output satisfies:

$$B[c,r] = A[r,c]$$

No floating-point arithmetic changes a value. The checker compares each float's 32-bit representation via bitwise equality, not a tolerance or checksum. It also reports maximum absolute error as a diagnostic; the mismatch count decides correctness. Positive and negative zero differ in bits even though their absolute difference is zero.

The harness limits N to 4096. In IEEE 754 single-precision float32, the significand has 24 bits of precision (23 stored bits plus 1 implicit leading bit). Integers from 0 through $2^{24} = 16{,}777{,}216$ have exact, unique representations. For $N=4096$, the matrix contains exactly $4096^2 = 16{,}777{,}216$ elements, meaning every single element receives a unique float value that avoids rounding ambiguity. Before each standalone validation launch, the harness fills the output with NaN bit patterns (0xff). A skipped store fails immediately instead of inheriting a correct value from a previous variant.

The supplied harness contains five variants:

copy isn't a transpose baseline. It performs a different operation with the same logical read and write volume. Treat it as context for the memory path, not as a denominator for a claimed transpose speedup.

Write the hypothesis before editing, reject incorrect candidates, then compare repeated timings under identical conditions. A surprising counter should change the next experiment, not become a story invented after the result.

Build the lab on an NVIDIA GPU

You need an NVIDIA GPU, a CUDA Toolkit with nvcc and Compute Sanitizer, and Nsight Compute's ncu command. Tool versions and GPU identity belong beside every result because compiler decisions, architecture, clocks, and profiler metric availability can change evidence.

Save the CUDA source and its host-checking header in the same directory. Save the recovery template for the deliberate failure later. The shell commands below target Bash on a CUDA-equipped machine. A CPU or Apple GPU can't execute these CUDA kernels; the later Python exercise checks indexing only.

Create an artifact directory and capture the environment. Keep pipefail enabled so a compiler or CUDA failure isn't hidden by a successful tee:

Compile optimized C++17 code with line information for attribution. -Xptxas=-v reports kernel resource use, including registers per thread, static shared memory per block, and local memory spills. -arch=native targets the GPU visible during compilation; this binary isn't a portable multi-architecture artifact. For a separate build host, select an explicit architecture supported by that Toolkit and the target GPU.^{[1]Reference 1NVIDIA CUDA Compiler Driver NVCChttps://docs.nvidia.com/cuda/cuda-compiler-driver-nvcc/}

Inspect the build log rather than guessing resource consumption:

• the compiler names each instantiated kernel;
• transpose_tiled<0> reserves $32 \times 32 \times 4 = 4096$ bytes of static shared memory;
• transpose_tiled<1> reserves $32 \times 33 \times 4 = 4224$ bytes;
• register count is recorded directly from ptxas;
• verify that spill stores and spill loads are both zero;
• no compilation error is present.

First exercise boundary shapes, including tiles that aren't full. Each command must exit zero:

Then run the fixed $4096 \times 4096$ performance case for 100 measured launches per normal variant:

The harness prints CSV columns variant, correct, mismatches, max_abs_error, average_ms, and effective_gb_s. Device and launch shape go to standard error. A fresh-buffer check runs before warmup and another after timing. Failure suppresses performance fields and returns exit code 2. The known-racy broken variant always returns 2 and is never timed, even if its sampled output matches. Runtime or argument errors return 1.

Pass conditions for the four normal variants are:

Don't rank variants yet. A single process run can be affected by clocks, device temperature, competing work, or startup state. Repeat separate processes under the same conditions, alternate order by round, and preserve every row:

These are repeated-buffer, warm-cache measurements, not streaming cold-DRAM measurements. The event interval excludes allocation, host transfers, and validation, but includes any device idle gaps between host-submitted launches. For tiny matrices, launch overhead can dominate. Profiler replay and cache-control settings can produce different cache conditions, so don't equate its counters with the unprofiled timing run automatically.

Compare repeated samples, not only the best row. Overlapping samples don't by themselves prove equality; report the variation and collect enough evidence to resolve a distinct difference. Don't claim a win from one favorable sample.

Read the naive memory path

CUDA coalesces a warp's global-memory requests into cache-line transactions that cover requested addresses. Modern NVIDIA architectures route global loads and stores through the L1/Texture cache and L2 cache down to DRAM channels. Hardware organizes memory traffic into 128-byte cache lines, which are subdivided into four 32-byte sectors.^{[2]Reference 2CUDA Programming Guide.https://docs.nvidia.com/cuda/cuda-programming-guide/} When threads in a warp access consecutive words within an aligned 128-byte segment, the hardware services all 32 threads with a single 128-byte cache request.

The naive kernel reads input with ideal spatial order, then scatters that order on output. Here is the indexing logic from transpose_naive:

Look closely at what happens across a single warp (32 threads where threadIdx.y is constant and threadIdx.x runs from 0 through 31):

1. Input read: input_index = (y + offset) * n + x. Because x = blockIdx.x * 32 + threadIdx.x, consecutive threads access consecutive float32 values in row y + offset. The 32 threads request $32 \times 4 = 128$ contiguous bytes. The memory subsystem serves all 32 lanes with a single 128-byte transaction (four adjacent 32-byte sectors). Coalescing efficiency is 100%.
2. Output store: output_index = x * n + (y + offset). Here x multiplies the stride $N$. Lane 0 writes to row 0, Lane 1 writes to row 1, Lane 2 writes to row 2, and Lane 31 writes to row 31. For $N=4096$, consecutive lanes write to addresses spaced by $4096 \times 4 = 16{,}384$ bytes ($16\text{ KB}$).
3. The transaction explosion: Because each thread's 4-byte float lands in a completely different 32-byte sector and 128-byte cache line, the memory controller can't coalesce them. It issues 32 separate 32-byte sector transactions to write just 128 bytes of payload.

That means transferring 1024 bytes across the memory bus to commit 128 bytes of data, wasting 87.5% of requested transfer volume. The write pipeline chokes on memory controller queues, and kernel throughput collapses.

Transpose does almost no useful floating-point arithmetic. Its logical data volume per launch is one read and one write for every element:

$$\text{logical bytes} = 2N^2 \cdot \operatorname{sizeof}(\text{float})$$

For $N=4096$, that's $134{,}217{,}728$ bytes, or 128 MiB. The harness reports:

$$\text{effective GB/s} = \frac{2N^2 \cdot 4}{\text{average time in seconds} \cdot 10^9}$$

Effective bandwidth isn't physical DRAM traffic. Cache hits, excess transactions, error-correction traffic, and architecture-specific paths make hardware traffic differ. Use it as a workload-normalized rate for comparing the same operation on the same setup.

A $4096 \times 4096$ transpose averages $0.250$ ms. What effective rate should the ledger record, and what hardware claim remains unsupported?

Answer

The logical volume is $2 \times 4096^2 \times 4 = 134{,}217{,}728$ bytes. Dividing by $0.000250$ seconds and then by $10^9$ gives about $536.9$ GB/s. That is a workload-normalized effective rate. It doesn't establish physical DRAM traffic or percentage of peak bandwidth without profiler counters and the matching hardware conditions.

Profile before editing

Nsight Compute groups kernel evidence into named sections. LaunchStats reports grid, block, registers, and shared memory. SpeedOfLight compares compute and memory resource throughput with device ceilings. MemoryWorkloadAnalysis traces traffic through the memory hierarchy. Occupancy reports active-warps capacity and its resource limiters.^{[3]Reference 3Nsight Compute Profiling Guidehttps://docs.nvidia.com/nsight-compute/ProfilingGuide/}

Collect the first measured naive launch, skipping its one validation launch and five warmups. The full section set can require replay, so use the unprofiled harness for timing. CLI filters and launch selection are documented separately from the profiling metrics.^{[4]Reference 4Nsight Compute CLIhttps://docs.nvidia.com/nsight-compute/NsightComputeCli/}

Read the report in this order:

1. Confirm kernel name, matrix size, block 32 x 8, and one profiled launch.
2. Check whether memory throughput is closer to its ceiling than compute throughput (dram__throughput.avg.pct_of_peak_sustained_elapsed).
3. Inspect global load and store sector efficiency (smsp__sass_average_data_bytes_per_sector_mem_global_op_st.pct). In transpose_naive, store efficiency plunges near 12.5% because only 4 bytes are utilized per 32-byte sector.
4. Check static shared memory is zero for the naive variant.
5. Read the occupancy limiter, but don't optimize it unless evidence connects the limiter to latency or throughput.

Low arithmetic intensity makes the memory path the primary suspect, not a confirmed bottleneck. Small launches can instead be dominated by overhead, and even large kernels can hit memory-instruction throughput before DRAM bandwidth. Roofline is a ceiling model, not a promise that every low-intensity kernel reaches peak bandwidth.^{[5]Reference 5Roofline: An Insightful Visual Performance Model for Multicore Architectureshttps://doi.org/10.1145/1498765.1498785} Use the profiler to identify which resource limits this launch.

Write one bound statement before editing:

If the report doesn't support a memory-path limit, revise the hypothesis before crediting a shared-memory rewrite with a speedup.

Coalesce both global directions

Shared memory is on-chip SRAM allocated per Streaming Multiprocessor (SM). It provides low latency and terabytes per second of aggregate bandwidth. The tiled transpose loads an input tile into shared memory with adjacent global reads. After a block barrier, threads read the tile with swapped indexes and write adjacent output words. Official CUDA guidance uses shared memory for this reordering pattern.^{[6]Reference 6CUDA C++ Best Practices Guidehttps://docs.nvidia.com/cuda/cuda-c-best-practices-guide/}

The core kernel below accepts Padding as a compile-time tile-width change. Variant tiled instantiates Padding = 0:

Notice the symmetry in coordinates:

• During input loading, thread threadIdx.x reads global index input_y * n + input_x and writes to tile[threadIdx.y + offset][threadIdx.x]. Because input_x advances with threadIdx.x, global reads are 100% coalesced.
• __syncthreads() guarantees all 256 threads in the thread block finish writing to tile before any thread reads from it.
• During output storing, thread threadIdx.x writes to global index output_y * n + output_x. Because output_x = blockIdx.y * kTile + threadIdx.x, output_x advances with threadIdx.x! The global store is now 100% coalesced into a single 128-byte cache transaction.

Keep the barrier outside bounds checks so every thread in the block reaches it, including threads processing edge tiles.

Check the edge-tile mapping without a GPU

For N=33, most threads in the last tile are out of bounds. Does an admitted output read still have a producer? Its shared coordinate is (threadIdx.x, threadIdx.y + offset). The corresponding input coordinate is (blockIdx.y * 32 + threadIdx.x, blockIdx.x * 32 + threadIdx.y + offset), exactly the transpose of the admitted output coordinate. The two bounds conditions agree.

This CPU model enumerates those coordinates, requires every shared read to have a producer, and checks that every output location is written exactly once. It also checks the bank arithmetic. It doesn't execute CUDA, model warp scheduling, or measure a memory transaction.

Run correctness and event timing first, then collect the profiler report:

Expected evidence shows two distinct shifts:

• global reads and writes now use adjacent lane addresses, curing global store serialization;
• the shared-memory column read triggers a severe bank conflict.

Shared memory is split into 32 hardware banks, each 4 bytes wide (32 bits). Successive 32-bit words map cyclically to successive banks: $\operatorname{bank}(k) = k \bmod 32$. When threads in a warp access distinct banks, requests complete in parallel in a single clock cycle. If multiple threads request different words from the same bank, the hardware serializes access across multiple cycles.

In the unpadded tile tile[32][32], row $r$, column $c$ sits at 1D float offset $32r + c$. When reading down column $c$, thread lane $i$ reads row $r = i$:

$$\operatorname{bank}(i) = (32i + c) \bmod 32 = c$$

Every single lane $0, 1, 2, \dots, 31$ requests a distinct word from Bank $c$. The hardware serializes the request into a 32-way bank conflict, taking 32 cycles instead of 1 cycle. The strided access penalty moved from global DRAM into on-chip shared memory.

Query the installed Nsight Compute metrics to inspect shared memory conflicts:

In the detailed Nsight report for transpose_tiled<0>, look for l1tex__data_bank_conflicts_pipe_lsu_mem_shared_op_ld.sum. The counter will report substantial shared load bank conflicts attributed directly to the line reading tile[threadIdx.x][threadIdx.y + offset].

Change the shared-memory stride

The next candidate changes only the tile declaration width:

By adding a single padding float to each row, the row stride increases from 32 to 33. The column element for lane $i$ at column $c$ now sits at word offset $33i + c$:

$$\operatorname{bank}(i) = (33i + c) \bmod 32 = (32i + i + c) \bmod 32 = (i + c) \bmod 32$$

Because $i$ ranges from 0 through 31, $(i + c) \bmod 32$ generates a complete cyclic permutation of the integers 0 through 31. Lane 0 accesses Bank $c$, Lane 1 accesses Bank $(c+1) \bmod 32$, and Lane 31 accesses Bank $(c+31) \bmod 32$.

Every lane hits a unique bank. The 32-way bank conflict vanishes completely, allowing all 32 lanes to complete their shared memory read in a single clock cycle. Padding doesn't alter matrix shape, output buffers, or logical global traffic. It changes address stride inside on-chip SRAM.^{[6]Reference 6CUDA C++ Best Practices Guidehttps://docs.nvidia.com/cuda/cuda-c-best-practices-guide/}

Would float tile[32][34] remove the column conflict as completely as float tile[32][33]? Derive the bank pattern before answering.

Answer

No. A 34-word row gives $\operatorname{bank}(i) = (34i + c) \bmod 32 = (2i + c) \bmod 32$. Because $\gcd(34, 32) = 2$, lanes 0 through 15 reach 16 distinct banks, then lanes 16 through 31 repeat those exact same banks. Each bank receives two different words, leaving a two-way conflict. A 33-word row works because $\gcd(33, 32) = 1$, so its stride is one modulo 32.

Check the padded result before timing or profiling it:

Compare identical evidence fields across tiled and padded:

Don't write that padding is faster without citing your hardware, software, shape, precision, timing method, baseline, and correctness record. Architecture and toolchain decide the size of the benefit, so a universal speedup claim would be unfounded.

A candidate is 6% faster, but its edit both pads the tile and changes the block from 32 x 8 to 32 x 16. What can the evidence ledger conclude, and how should the experiment recover?

Answer

The run shows that the combined candidate differs, but it can't attribute the change to padding or launch geometry. Restore 32 x 8 and measure padding alone against the fixed baseline. Then hold tile width fixed and test 32 x 16 separately. Each branch still needs the same exact-output, profiler, sanitizer, and repeated-timing gates.

Loop unrolling and instruction-level parallelism

Hardware needs work in flight to hide memory latency. On modern GPUs, fetching a line from off-chip DRAM takes between 200 and 800 clock cycles. There are two primary mechanisms to tolerate that latency: warp-level parallelism (having enough active warps on the SM to switch between) and instruction-level parallelism (ILP, having multiple independent instructions in flight within the same thread).

In transpose_tiled, each block contains $32 \times 8 = 256$ threads, but the tile holds $32 \times 32 = 1024$ floats. Each thread processes four rows separated by kBlockRows = 8:

The #pragma unroll directive instructs the compiler to unroll all four iterations into straight-line SASS instructions. This transformation delivers two distinct benefits:

1. Eliminating loop overhead: It removes induction variable increments, compare instructions, and branch jumps.
2. Exposing independent loads to the scheduler: The compiler issues four independent global load instructions (LDG.E) back-to-back before any dependent write to shared memory. The memory pipeline dispatches these four memory requests concurrently, overlapping their latency windows.

Unrolling increases register pressure. Each in-flight memory operation requires dedicated registers to hold target addresses and loaded values. Inspect the compiler resource output from nvcc -Xptxas=-v:

If register usage climbs too high (for example, above 32 or 48 registers per thread depending on architecture), the GPU can't fit as many concurrent warps on the SM. If registers exceed the allocation limit, the compiler spills variables to local memory (spill stores and spill loads). Local memory resides in off-chip DRAM, introducing cache thrashing and destroying throughput. Always verify that spill stores and loads remain zero.

Branch divergence is minimal in this kernel. For interior tiles where $x < n$ and $y + 32 \le n$, every thread in the warp evaluates the bounds check to true. On boundary tiles ($N=33$ or $N=65$), threads beyond matrix limits take the false branch. Modern architectures with Independent Thread Scheduling execute divergent branches sequentially using thread masks, so keeping work outside bounds guards ensures maximum SIMT lane utilization.

Occupancy is a constraint, not a score

Occupancy is the ratio of active warps on a Streaming Multiprocessor (SM) to the hardware maximum. More resident warps can hide latency, but highest occupancy doesn't automatically deliver highest performance. Registers, shared memory, threads per block, and barriers can limit resident blocks, while forcing resource use down can introduce spills or extra instructions.^{[6]Reference 6CUDA C++ Best Practices Guidehttps://docs.nvidia.com/cuda/cuda-c-best-practices-guide/[3]Reference 3Nsight Compute Profiling Guidehttps://docs.nvidia.com/nsight-compute/ProfilingGuide/}

The lab block uses 32 x 8 = 256 threads, or eight 32-thread warps. The concrete resource comparison across variants shows:

The 128-byte increase from padding may leave occupancy unchanged or cross an allocation boundary for a particular resource configuration. Read the report. Don't infer resident blocks from source alone.

Use these diagnoses when evaluating profiler evidence:

Keep the padded variant when it removes the measured conflict and improves repeat timing without correctness or resource regression. Revert when the counter improves but end-to-end kernel timing doesn't. Counter movement isn't a product outcome by itself.

Deliberately remove the barrier

Performance edits often break synchronization before they break memory bounds. Variant broken uses the padded tile but omits __syncthreads() between shared writes and shared reads. Some warps read shared memory locations before producer warps finish writing them.

Run a smaller matrix under the tools. Compute Sanitizer documentation recommends running memcheck before racecheck because racecheck doesn't check invalid memory addresses.^{[7]Reference 7Compute Sanitizerhttps://docs.nvidia.com/compute-sanitizer/ComputeSanitizer/} The harness intentionally exits 2 for broken; capture that expected application failure without aborting the diagnostic sequence. Keep the full logs: --error-exitcode 99 applies to detected sanitizer errors when the application itself would otherwise succeed. A nonzero exit alone isn't evidence of a race: tool startup failures and application errors also fail the command. Identify the reported shared-memory hazard in the racecheck log.

Expected diagnostic split:

• memcheck reports no illegal addresses because every index stays within allocation bounds;
• racecheck reports shared-memory read-after-write (RAW) hazards around the missing barrier;
• synccheck reports no misuse because the kernel omitted a barrier instead of executing an invalid barrier;
• output checker may fail, but an accidental correct output on one run doesn't clear the race.

That last point is critical. Scheduler timing can hide data races during an ordinary run. Dynamic race evidence invalidates a kernel even when sampled output happens to match.

Recover by restoring the block barrier. Validate every candidate you intend to keep, including a partial tile and a full-tile case. These checks run outside the performance measurement:

The recovery record connects the symptom to the ordering requirement:

Fill recovery-template.md and keep it with the raw reports. A reviewer should be able to trace the failure, repair, tested shapes, and post-recovery cost. Clean dynamic checks cover those executions, not every possible input or schedule.

Decide between transpose, reduction, and scan

Shared-memory tiling works well for spatial reordering. It isn't a universal answer for every memory-bound kernel. Ask what output dependency requires before choosing a primitive:

Consider moving the correctness check from CPU back to GPU. If you need only total mismatch count and maximum absolute error, use reductions. Both outputs collapse $N^2$ comparisons into one or two scalars. A scan would do extra work by producing a prefix result for every comparison.

Now consider compacting indexes of mismatched elements. Each failing element needs a destination position based on the number of failures before it. That is a prefix dependency, so a scan computes offsets. A reduction can count total failures but can't assign each failure a unique compacted slot.

The decision rule remains mechanical:

1. Same number of values, different layout: transpose.
2. Fewer summary values: reduction.
3. Same number of prefix states or compaction offsets: scan.

Tiling can appear inside all three implementations, but the dependency graph picks the algorithm. Don't start by copying a transpose tile into a reduction or scan kernel.

Report the result you measured

Complete the comparison only after the variants share the same workload and validity contract. The entries below are hypotheses and required observations, not supplied benchmark results:

One row isn't complete if its correctness cell says only that it looks right. Keep raw CSV, .ncu-rep files, sanitizer logs, compiler resource output, environment files, source hash, and the filled recovery record.

You've finished when you can defend all six claims:

• workload and timing method stayed fixed;
• baseline bound came from profiler evidence;
• each candidate changed one causal mechanism;
• output check matched transpose exactly;
• sanitizer results are clean for the accepted kernels and tested shapes;
• reported performance includes conditions and distribution, not unsupported headline speedup.