GPU Kernel Engineering: Why C = A + B Is Not Enough (Series Part 2)

A single line of arithmetic, C[i] = A[i] + B[i], can keep an NVIDIA H100 at well under one percent of its advertised floating-point speed, and no amount of clever mathematics will fix it. The chip is not slow. It is hungry, and the kitchen that feeds it is far too small for the dining room. Closing that gap is what GPU kernel engineering is about: writing the small program that runs on thousands of threads so that data arrives at the arithmetic units before they go idle.

This matters in 2026 because every large language model, every computational photography pipeline and every simulation workload is, at the bottom, a stack of kernels. The difference between a kernel that moves bytes well and one that does not is routinely a factor of ten, and it shows up directly as GPU rental cost.

This is Part 2 of the series “GPU Compute: From iPhone to H100”. In Part 1, the anatomy of the SM, block, warp and thread hierarchy, we built the hardware picture. Here we write the software that runs on it, and find out why the obvious program starves.

What this covers: what a kernel is (5W1H), thread indexing and launch configuration with compilable CUDA C++, arithmetic intensity and the roofline, data starvation, latency hiding, coalescing, divergence, and a one-page cheat sheet.

Context and Background

Part 1 ended with a physical hierarchy. An H100 SXM5 chip contains 132 Streaming Multiprocessors (SMs, the independent compute tiles of an NVIDIA GPU). Each SM can hold up to 2,048 resident threads, grouped into warps of 32 threads that execute together, and a block may contain at most 1,024 threads. Those limits come from NVIDIA’s compute capability 9.0 table in the CUDA Programming Guide, which also lists 64 resident warps and 32 resident blocks per SM, 64K 32-bit registers, and up to 228 KB of shared memory per SM.

Hardware alone does nothing. Someone must decide what each of those roughly 270,000 resident threads does, which memory address it touches, and how the work is carved into blocks. That decision is the kernel, and the person making it is practising kernel engineering. Historically this lived in the high-performance computing (HPC) community; today it is the job of the engineers behind libraries such as cuBLAS, FlashAttention and the inference engines compared in our vLLM, SGLang and TensorRT-LLM benchmark on H100.

The economic stakes sit in the software layer as much as in silicon. CUDA, NVIDIA’s programming model, is widely described as the moat that keeps rivals at bay, because years of tuned kernels are hard to replicate. We examine that argument in the Qualcomm and Modular CUDA moat analysis. Understanding what a kernel actually does is the quickest way to judge whether that moat is real.

A word on the phone side of our series title. Apple GPUs run the same idea under different names: Apple documents threadgroups and SIMD-groups, and Metal kernels use thread positions in a grid. The analogy to CUDA is close but not one-to-one, so in this post we stay with CUDA C++, where the vocabulary and the published limits are most precisely documented.

The School Analogy Continued: What a Kernel Is

In Part 1 we used a school. The SM is a department. A block is a classroom of 256 students. A warp is a row of 32 desks that the teacher addresses at once. The warp scheduler is that teacher, who looks around the room and picks which row to teach next. Now add the missing piece: a kernel is the homework sheet. Every student receives the same sheet, but each writes down an answer for a different question number, because each student knows their own seat number.

That sentence contains the whole programming model. A kernel is one function. You launch it once from the host (the CPU), and the hardware instantiates it across thousands of threads. Each thread discovers its own identity from built-in variables and uses that identity to choose its data. There is no loop over elements in the kernel; the loop is the launch.

The analogy also predicts the failure we are about to study. Students can only answer when they have the textbook page in front of them. The textbook lives in the library across campus (device memory, called HBM for High Bandwidth Memory, a stack of DRAM chips beside the die). A student’s desk drawer is a register, and the classroom shelf is shared memory. If every student has to walk to the library for every single question, the teacher stands idle while the room waits. That is data starvation, and for simple kernels it is the default, not the exception.

What, Why, Who, When, Where, How of GPU Kernel Engineering

The 5W1H frame is a quick way to pin the idea down before the details.

What. GPU kernel engineering is the practice of designing, launching and tuning the functions that execute across GPU threads, with the explicit goal of matching three things: the shape of the parallel work, the shape of the hardware, and the shape of the memory traffic.

Why. Because peak numbers are not delivered numbers. The H100 SXM page on NVIDIA’s site lists 67 teraFLOPS of FP32 (32-bit floating point) and 3.35 TB/s of memory bandwidth. Whether a given program sees either number depends entirely on how it is written.

Who. Library authors at NVIDIA and in the open-source community, inference-engine teams, scientific programmers, and increasingly application engineers who write one custom fused operation to remove a bottleneck.

When. When a profiler shows an existing library call is not good enough: an unusual shape, an operation sequence that library calls cannot fuse, or a memory-bound step that dominates the timeline. Writing a kernel before profiling is usually premature.

Where. On the device. Kernels run on the GPU, are launched from host code, and operate on device memory. The data transfers over PCIe or NVLink between host and device are a separate bottleneck that kernel tuning cannot remove.

How. By choosing an indexing scheme, a launch configuration, a memory access pattern, and an amount of work per thread, then iterating against measurements. The rest of this post is that “how”.

Anatomy of a CUDA Kernel: From Launch to Index

Figure 1: How a CUDA kernel launch becomes work on the hardware. The grid is a software idea; blocks are assigned to SMs; each block is split into warps that the scheduler issues.

Figure 1 summarises the path. The host code names a grid size and a block size; the hardware distributes blocks to whichever SMs have room; each block is sliced into warps of 32; and the warp scheduler issues instructions to whichever warp is ready.

Here is the same hierarchy in plain text, with the indexing variables that each level exposes.

HOST: vecAdd<<<gridDim, blockDim>>>(A, B, C, N)
        |
        v
GRID  (gridDim.x blocks, one kernel launch)
 |-- Block 0   blockIdx.x = 0   (blockDim.x threads, e.g. 256)
 |     |-- Warp 0   threads  0..31   (lanes 0..31)
 |     |-- Warp 1   threads 32..63
 |     |-- ...      Warp 7  threads 224..255
 |-- Block 1   blockIdx.x = 1   -> runs on SOME SM, any order
 |-- ...
 +-- Block G-1

Global index of a thread:
   i = blockIdx.x * blockDim.x + threadIdx.x

Hardware view (H100 SXM5, 132 SMs):
 Chip --> SM (up to 32 blocks, 64 warps, 2048 threads resident)
          --> 4 warp schedulers, each picks a READY warp per cycle
               --> 32 lanes execute one instruction together

The four warp schedulers per SM come from NVIDIA’s Hopper architecture documentation; the launch syntax and indexing variables come from the CUDA Programming Guide. The global index formula is the most-typed line in GPU programming, and the next section explains why it is safe.

Three rules the model gives you

First, blocks are independent. The runtime may execute them in any order and on any SM, so a correct kernel never assumes block 1 runs after block 0. This freedom is what lets the same binary scale from a small laptop GPU to an H100: a bigger chip simply runs more blocks at once.

Second, threads within a block can cooperate. They share fast on-chip memory and can synchronise with __syncthreads(). Threads in different blocks cannot synchronise inside a kernel in the ordinary model, apart from cooperative-launch features and Hopper’s thread block clusters.

Third, a warp executes in lockstep for ordinary instructions. All 32 lanes issue the same instruction. If lanes disagree about a branch, the warp serialises the paths. This is the root of thread divergence, which we return to below.

Launch Configuration: Choosing Grid and Block Size

The expression between the triple angle brackets is the kernel launch configuration. <<<grid, block>>> sets how many blocks and how many threads per block. The total thread count is the product, and the number of elements usually dictates the grid size.

Why 256 is the folk default

Block size has hard bounds and soft preferences. The hard bound on current NVIDIA GPUs is 1,024 threads per block. The soft preference is a multiple of the warp size of 32, because a block of 100 threads still occupies four warp slots and wastes 28 lanes in the last warp.

Beyond that, the choice is a trade. Larger blocks amortise per-block overhead and allow more cooperation through shared memory, but they are coarser to schedule and consume registers and shared memory in bigger chunks. With a limit of 2,048 resident threads per SM, 256-thread blocks allow up to eight blocks per SM, enough that when one block waits on a barrier or a memory load, others fill in. A block of 1,024 threads allows only two. So 128 to 256 is a good first guess, and measurement decides the rest.

Computing the grid size correctly

For N elements and a block of B threads, the grid needs ceil(N / B) blocks. In integer arithmetic that is (N + B - 1) / B. The last block usually contains threads whose index lands beyond N, which is why every well-formed kernel contains a bounds check. Skipping it writes to memory you do not own, and the symptoms are the cruel kind: wrong answers on some sizes only.

Worked Example: The Naive Vector Add in Real CUDA C++

Here is a complete, compilable program. Save it as vecadd.cu and build with nvcc -O3 -o vecadd vecadd.cu. It uses an error-checking macro, because silent CUDA failures are a common source of bad benchmarks.

#include <cstdio>
#include <cstdlib>
#include <cuda_runtime.h>

#define CUDA_CHECK(call)                                                   \
    do {                                                                   \
        cudaError_t err__ = (call);                                        \
        if (err__ != cudaSuccess) {                                        \
            std::fprintf(stderr, "CUDA error %s at %s:%d\n",               \
                         cudaGetErrorString(err__), __FILE__, __LINE__);   \
            std::exit(EXIT_FAILURE);                                       \
        }                                                                  \
    } while (0)

// One thread, one element. Bounds check protects the ragged last block.
__global__ void vecAdd(const float* A, const float* B, float* C, size_t N) {
    size_t i = static_cast<size_t>(blockIdx.x) * blockDim.x + threadIdx.x;
    if (i < N) {
        C[i] = A[i] + B[i];
    }
}

int main() {
    const size_t N = 1ull << 28;               // 268,435,456 elements
    const size_t bytes = N * sizeof(float);    // 1 GiB per array

    float *hA = (float*)std::malloc(bytes);
    float *hB = (float*)std::malloc(bytes);
    float *hC = (float*)std::malloc(bytes);
    for (size_t i = 0; i < N; ++i) { hA[i] = 1.0f; hB[i] = 2.0f; }

    float *dA, *dB, *dC;
    CUDA_CHECK(cudaMalloc(&dA, bytes));
    CUDA_CHECK(cudaMalloc(&dB, bytes));
    CUDA_CHECK(cudaMalloc(&dC, bytes));
    CUDA_CHECK(cudaMemcpy(dA, hA, bytes, cudaMemcpyHostToDevice));
    CUDA_CHECK(cudaMemcpy(dB, hB, bytes, cudaMemcpyHostToDevice));

    const int block = 256;
    const unsigned grid = static_cast<unsigned>((N + block - 1) / block);

    cudaEvent_t t0, t1;
    CUDA_CHECK(cudaEventCreate(&t0));
    CUDA_CHECK(cudaEventCreate(&t1));

    vecAdd<<<grid, block>>>(dA, dB, dC, N);    // warm-up launch
    CUDA_CHECK(cudaGetLastError());

    CUDA_CHECK(cudaEventRecord(t0));
    vecAdd<<<grid, block>>>(dA, dB, dC, N);
    CUDA_CHECK(cudaEventRecord(t1));
    CUDA_CHECK(cudaEventSynchronize(t1));

    float ms = 0.0f;
    CUDA_CHECK(cudaEventElapsedTime(&ms, t0, t1));
    double gbMoved = 3.0 * bytes / 1e9;        // read A, read B, write C
    std::printf("time %.3f ms, effective bandwidth %.1f GB/s\n",
                ms, gbMoved / (ms / 1e3));

    CUDA_CHECK(cudaMemcpy(hC, dC, bytes, cudaMemcpyDeviceToHost));
    std::printf("C[0] = %.1f, C[N-1] = %.1f\n", hC[0], hC[N - 1]);

    cudaFree(dA); cudaFree(dB); cudaFree(dC);
    std::free(hA); std::free(hB); std::free(hC);
    return 0;
}

Three details deserve attention. The index is computed in size_t, because with 2^28 elements and larger arrays a 32-bit int product can overflow. The timing uses CUDA events and a warm-up launch, because kernel launches are asynchronous and the first launch pays one-time costs. And the program reports effective bandwidth, not FLOPS, which is the right metric for this kernel, as the next section shows.

Why C = A + B Is Not Enough: Arithmetic Intensity and the Roofline

The vector add looks efficient: every thread is busy, there are no branches, and the grid fills the chip. Yet it is a textbook example of a memory-bound kernel. To see why, count what each thread does. It loads two 4-byte floats, performs one addition, and stores one 4-byte float. That is 12 bytes of traffic for 1 floating-point operation (FLOP).

The ratio of work to traffic is called arithmetic intensity, measured in FLOPs per byte. For fp32 vector add it is 1/12, about 0.083 FLOP per byte. Every other number in this section follows from that single fraction.

Figure 2: Sequence flow of one warp in the naive vector add. The warp issues two loads, then stalls for the round trip to HBM; only after both operands return does the single addition happen.

The roofline in four lines

The roofline model, introduced by Williams, Waterman and Patterson in a 2009 Communications of the ACM article, says attainable performance is the minimum of two ceilings: the compute peak, and the bandwidth multiplied by arithmetic intensity. A kernel left of the “ridge point” is bandwidth-limited; a kernel to the right is compute-limited.

Use NVIDIA’s published H100 SXM figures: 67 TFLOPS of FP32 and 3.35 TB/s of HBM3 bandwidth. The ridge point is the intensity at which the two ceilings meet, so 67e12 divided by 3.35e12 gives roughly 20 FLOP per byte (derived). The vector add sits at 0.083, about 240 times to the left of the ridge.

Attainable performance for the vector add is therefore 3.35e12 bytes/s times 1/12 FLOP/byte, about 2.8e11, or 279 GFLOPS (derived). Against a 67 TFLOPS peak, that is 0.42 percent. The arithmetic units could do 240 times more additions per second; there is simply no data to feed them. This is the precise sense in which the line C[i] = A[i] + B[i] is “not enough”: correctness is satisfied, and the hardware is 99.6 percent idle.

Worked numbers for the example program

With N = 2^28 = 268,435,456 elements, each array is 1 GiB (1,073,741,824 bytes) and the kernel moves 3 x 1.07 GB, about 3.22 GB. At the published 3.35 TB/s the floor is 3.22e9 / 3.35e12, about 0.96 ms (derived, theoretical). Real kernels do not reach the datasheet peak; sustained fractions in the range of roughly 80 to 90 percent are commonly reported for well-written streaming kernels, so expect on the order of one millisecond or slightly above. That is a typical-range statement, not a measurement of any specific system: run the program above and read your own number.

The key insight is what you would do next. Rewriting the addition in assembly makes no difference. Using the faster tensor cores makes no difference. The only levers are moving fewer bytes, moving them more efficiently, or reusing them once they are on chip.

Data Starvation: When the Teacher Has Nothing to Teach

Return to the classroom. The teacher (warp scheduler) is ready to issue an instruction every clock cycle, but each instruction needs a textbook page from the far library. If a page takes hundreds of cycles to arrive, then a teacher with only one row to teach sits idle for hundreds of cycles per addition. Data starvation is exactly this: functional units stall because operands are not yet present.

Two numbers define the problem: bandwidth, the width of the road to the library, and latency, the length of the walk. GPUs attack the latency half with a trick that CPUs mostly avoid, and it is the heart of kernel design.

Latency hiding through occupancy

A CPU hides memory latency with large caches and speculation on a few fast threads. A GPU hides it with a crowd. While warp 3 waits for its page, the scheduler issues instructions from warp 7, which already has its data; then warp 12; and so on. As long as some warp is ready every cycle, the arithmetic units never see the wait. Memory is slow, but the machine is never idle.

Occupancy is the measure of this crowd: the ratio of resident warps on an SM to the maximum the SM supports, which is 64 warps on compute capability 9.0. Higher occupancy gives the scheduler more candidates. But occupancy is bounded by three resources per SM: threads (2,048), registers (64K 32-bit registers), and shared memory (up to 228 KB). A kernel that uses 128 registers per thread can fit only 512 threads on an SM (65,536 divided by 128), which is 25 percent occupancy, even if the block size looks fine.

Little’s law sets the target

How many bytes must be in flight to saturate memory? Little’s law from queueing theory gives the answer: bytes in flight equal bandwidth multiplied by latency. Suppose, purely as an illustration, that a global memory access takes about 500 nanoseconds under load; this number is an assumption and not a published H100 figure. Then 3.35e12 bytes/s times 500e-9 s is about 1.7 MB in flight across the chip, or roughly 12.7 KB per SM across 132 SMs (derived, illustrative).

At 12 bytes per thread-element, about a thousand elements must be in flight per SM. That is comfortably within a 2,048-thread SM, which is why the naive vector add can reach a high fraction of bandwidth despite its naivety. Heavier kernels, which use many registers, cannot rely on this and need explicit prefetching or asynchronous copies.

Second-order thinking: what breaks when you scale up

The tempting conclusion is “more occupancy is always better”. That assumption breaks. Pushing registers down to raise occupancy causes register spilling to local memory, which is slow and itself consumes bandwidth. Higher occupancy also shrinks the cache and shared memory available per thread. Expert kernels often run at modest occupancy with more independent work per thread (instruction-level parallelism) to hide latency instead.

A second broken assumption: “the problem is the same on a bigger GPU”. Moving from an A100 to an H100 raises both bandwidth and compute, but the compute peak has historically grown faster than bandwidth, so the ridge point drifts right and more kernels become memory-bound on each generation. A kernel that was balanced last generation may be starved on the next. Plan for the ratio, not the absolute number.

Coalescing: How Warps Read Memory

The previous sections assumed bytes move at the peak. That only happens when the 32 threads of a warp read neighbouring addresses. Global memory is served in aligned segments, in units that on recent NVIDIA GPUs are 32-byte sectors, and the memory system merges the 32 per-thread loads of a warp into as few sectors as possible. This merging is memory coalescing.

Figure 3: A coalesced warp touches four consecutive 32-byte sectors for 128 useful bytes. A stride-32 pattern touches 32 different sectors and uses 4 bytes from each.

In the vector add, thread i reads A[i], so a warp reads 32 consecutive floats: 128 contiguous bytes, which is four 32-byte sectors, with 100 percent of the fetched bytes used. Now compare a kernel where thread i reads A[i * 32]. The same warp touches 32 different sectors, fetching 1,024 bytes to use 128. That is 12.5 percent efficiency, an eightfold waste of precious bandwidth, from a one-token change in an index expression (derived from the 32-byte sector size).

// Coalesced: consecutive threads touch consecutive floats.
__global__ void copyCoalesced(const float* in, float* out, size_t n) {
    size_t i = static_cast<size_t>(blockIdx.x) * blockDim.x + threadIdx.x;
    if (i < n) out[i] = in[i];
}

// Strided: consecutive threads are `stride` floats apart.
// Same arithmetic, up to ~8x more memory sectors fetched per useful byte.
__global__ void copyStrided(const float* in, float* out, size_t n, int stride) {
    size_t i = static_cast<size_t>(blockIdx.x) * blockDim.x + threadIdx.x;
    size_t j = (i * stride) % n;
    if (i < n) out[i] = in[j];
}

The pattern recurs everywhere in real software. A matrix stored row-major is read well by threads that walk along a row, and badly by threads that walk down a column. An array-of-structures layout gives each thread a strided view of one field; a structure-of-arrays layout makes the same access contiguous. The fix is almost always a data layout change, or a staged transpose through shared memory.

The Grid-Stride Loop: Decoupling Problem Size from Launch Size

The one-thread-per-element kernel ties the grid to N. That is fine for a demo, but production kernels often use a grid-stride loop: launch a grid sized to the hardware, and let each thread walk the array in steps of the total thread count.

__global__ void vecAddStride(const float* A, const float* B, float* C, size_t N) {
    size_t stride = static_cast<size_t>(gridDim.x) * blockDim.x;
    for (size_t i = static_cast<size_t>(blockIdx.x) * blockDim.x + threadIdx.x;
         i < N; i += stride) {
        C[i] = A[i] + B[i];
    }
}

// Launch sized to the device, not to N:
//   int dev = 0, sms = 0;
//   cudaGetDevice(&dev);
//   cudaDeviceGetAttribute(&sms, cudaDevAttrMultiProcessorCount, dev);
//   vecAddStride<<<sms * 8, 256>>>(dA, dB, dC, N);

Why does this help? First, the loop stays coalesced: in each iteration, consecutive threads still touch consecutive elements. Second, launch overhead and block-scheduling cost are paid once per SM-sized wave instead of once per few thousand elements. Third, the same kernel is correct for any N, including N smaller than the grid, and it gives you a natural place to add per-thread work such as loading a vector of four floats at a time.

The trade-off is subtle. Fewer, longer-lived blocks reduce scheduling overhead but can cause tail effects: if work per thread is uneven, the last blocks finish at different times and some SMs idle. For a uniform streaming kernel this does not matter. For irregular work it can, and a larger grid with dynamic scheduling does better. As with occupancy, the right answer is measured, not assumed.

Vectorised loads: fewer instructions for the same bytes

A thread can load 16 bytes in one instruction by treating four floats as a float4. This reduces the number of load instructions by four and tends to improve the number of bytes in flight per thread, helping with latency hiding. It requires the pointer to be 16-byte aligned (memory from cudaMalloc is aligned to at least 256 bytes) and the length to be handled when N is not a multiple of four. It does not change the 12 bytes per FLOP, so it does not move the kernel right on the roofline; it only helps it slide closer to the bandwidth ceiling.

Warp Divergence: When the Row of Desks Disagrees

Coalescing concerns where lanes read. Thread divergence concerns what lanes do. All 32 threads of a warp share one instruction stream. If an if statement sends some lanes one way and others another, the hardware executes the taken path with the non-participating lanes masked off, then executes the other path with the roles reversed. The warp pays for both paths.

In the classroom, the teacher says “students with even seat numbers do problem A, odd seat numbers do problem B”. The teacher cannot teach both at once to one row. She teaches A while the odd seats sit, then B while the even seats sit. Throughput for that stretch halves.

Figure 4: A warp that diverges executes both branch paths one after the other with lanes masked, then reconverges. Cost is the sum of the paths, not the maximum.

// Divergent: even and odd lanes of every warp take different paths.
__global__ void divergent(float* x, size_t n) {
    size_t i = static_cast<size_t>(blockIdx.x) * blockDim.x + threadIdx.x;
    if (i < n) {
        if (threadIdx.x % 2 == 0) x[i] = x[i] * 2.0f;
        else                      x[i] = x[i] + 1.0f;
    }
}

// Not divergent at warp granularity: the condition is uniform within a warp.
__global__ void warpUniform(float* x, size_t n) {
    size_t i = static_cast<size_t>(blockIdx.x) * blockDim.x + threadIdx.x;
    if (i < n) {
        if ((threadIdx.x / 32) % 2 == 0) x[i] = x[i] * 2.0f;
        else                             x[i] = x[i] + 1.0f;
    }
}

The second kernel does the same amount of branching in aggregate, but each warp is wholly on one side, so no warp executes both paths. The lesson is to arrange data so that condition outcomes are correlated within each group of 32 consecutive threads: sort or bucket records by type, partition work so each warp handles one case, or replace short branches with predicated arithmetic that the compiler can turn into a select instruction.

Note the bounds check if (i < n) is itself a branch. It diverges only in the final warp of the final block, which is why it is essentially free. Divergence matters when it is frequent, and when both paths are long. For a memory-bound kernel such as the vector add, a bit of divergence is often hidden under memory latency; for a compute-bound kernel, it is a direct tax. NVIDIA’s Volta architecture added independent thread scheduling, so lanes of a diverged warp can interleave and make progress independently, which fixes certain deadlock patterns but does not make divergence free.

Putting It Together: One Kernel’s Life, Step by Step

Here is the sequence in order, so you can place each concept on a timeline.

The host allocates device memory with cudaMalloc and copies inputs across PCIe or NVLink. This transfer is outside the kernel and can dwarf it for small problems.
The host launches vecAdd<<<grid, block>>>. The launch is asynchronous; the CPU continues while the GPU queues the work.
The block scheduler hands blocks to SMs with free thread, register and shared memory budget. Blocks that do not fit wait.
Each block splits into warps. Warp schedulers pick any warp whose operands are ready.
A warp issues two global loads. The memory system coalesces them into sectors and fetches from L2 or HBM; the warp stalls, and the scheduler switches to another warp.
When both operands arrive, one add executes across 32 lanes. A store goes back to memory.
When all blocks retire, the kernel is complete; the host synchronises (an event or cudaDeviceSynchronize) and copies results back.

Notice how little of that timeline is arithmetic. Step 6 is a single instruction wrapped in six steps of logistics. That ratio, not the speed of the add, is what kernel engineering optimises.

Diagnosing Which Wall You Hit

A practical question is how to tell starvation from other limits. Three measurements discriminate most cases. Compare achieved bandwidth to the peak: if you are at 80 percent or more of 3.35 TB/s, the kernel is saturating memory and only moving fewer bytes helps. Compare achieved FLOPS to the peak: if you are far below both, the problem is latency or occupancy, not either throughput ceiling. Finally, check the stall reasons that a profiler such as NVIDIA Nsight Compute reports per warp, which separates waiting on memory from waiting on barriers, dependencies or instruction fetch.

The decision flow is mostly mechanical. Low bandwidth and low occupancy suggests insufficient parallelism or too many registers. Low bandwidth with high occupancy suggests uncoalesced access. High bandwidth with low FLOPS suggests the kernel is memory-bound and wants fusion or reuse. High FLOPS near peak suggests you are done, or that you should consider the tensor cores for matrix work.

Trade-offs, Gotchas, and What Goes Wrong

Benchmarking traps. Timing a kernel with a CPU clock without synchronising measures only the launch. Timing a first launch includes initialisation. Measuring a small array that fits in the L2 cache reports cache bandwidth, not HBM bandwidth; the H100 has a large L2 (NVIDIA lists 50 MB for the SXM part), so an array under that size can look dramatically faster than the 3.35 TB/s HBM figure. Use arrays well beyond L2 when measuring memory throughput.

Peak numbers are conditions. The 67 TFLOPS FP32 figure assumes every lane issues a fused multiply-add every cycle at boost clock; the bandwidth figure assumes ideal streaming. Neither is a promise for your kernel, and power or thermal limits can lower clocks under sustained load.

Over-tuning to one chip. A block size, tile size or unroll factor tuned to an H100 may be wrong on another generation or on a data-centre part with a different SM count. Hardcoding 132 anywhere in your code is a bug waiting to happen; query the device attributes instead.

Correctness hazards. Forgetting the bounds check, using a 32-bit index on large arrays, racing on shared memory without __syncthreads(), and calling __syncthreads() inside divergent code (which can hang the block) are the classic bugs. Always run compute-sanitizer on new kernels.

Premature kernel writing. A hand-written kernel is a maintenance liability. If cuBLAS, cuDNN or a fused library kernel covers the case, use it. Custom kernels pay off for unusual shapes or for fusing steps whose intermediate results would otherwise round-trip through HBM, which is the subject of Part 3.

Ignoring the host side. If the host enqueues many tiny kernels, launch overhead and PCIe transfers dominate. The remedies there are CUDA Graphs, batching, and keeping data resident on the device.

Practical Recommendations

Begin with arithmetic intensity, not with code. Before touching a kernel, write down bytes moved and FLOPs performed, compute the intensity, and compare with the ridge point of your target chip (about 20 FLOP per byte for H100 FP32, from published figures). That one calculation tells you whether the target is bandwidth or compute, and it takes two minutes.

Then follow a short, ordered routine. Get a correct baseline with a bounds check and size_t indexing. Measure effective bandwidth with CUDA events. Fix access patterns before anything else, because uncoalesced loads cost integer factors. Only then tune block size, try a grid-stride loop with vector loads, and examine occupancy. Fuse operations last, because it changes the program structure the most.

Pick a block size that is a multiple of 32; start at 128 or 256.
Size grids from the device’s SM count when using grid-stride loops.
Keep consecutive threads on consecutive addresses.
Keep conditions uniform within groups of 32 threads where possible.
Measure with arrays far larger than L2, with warm-up and events.
Profile with Nsight Compute before changing anything.

One-Page Cheat Sheet

Concept	One-line meaning	School analogy	Number to remember (H100 SXM)	Typical fix
Kernel	One function run by many threads	Homework sheet	Launched via `<<<grid, block>>>`	Launch once, not in a loop
Global index	`blockIdx.x * blockDim.x + threadIdx.x`	Seat number across the school	Use `size_t` for big N	Bounds check always
Block size	Threads per block	Classroom size	Max 1,024; start 128-256	Multiple of 32
Warp	32 lockstep threads	Row of 32 desks	64 resident warps per SM	Avoid partial warps
Resident limits	What one SM can hold	Department capacity	2,048 threads, 32 blocks, 64K registers	Watch register use
Arithmetic intensity	FLOPs per byte moved	Work per page fetched	fp32 add: 1/12 FLOP/byte	Reuse data on chip
Ridge point	Intensity where compute meets bandwidth	Break-even point	About 20 FLOP/byte (derived)	Know which side you are on
Data starvation	Units idle awaiting data	Teacher with no textbooks	3.35 TB/s peak HBM3	Move fewer bytes
Occupancy	Resident warps over maximum	Rows ready to teach	Max 64 warps per SM	Balance with registers
Coalescing	Warp reads contiguous sectors	One trip to the library	32-byte sectors	Contiguous layout, SoA
Divergence	Lanes take different branches	Teacher teaches A then B	Cost is sum of paths	Make conditions warp-uniform
Grid-stride loop	Fixed grid walks any N	Students take turns on a long sheet	Grid near SMs x 8 blocks	Add vector loads

Conclusion: From Starved Kernels to the Physics of Memory

Three lessons carry forward. A kernel is one function instantiated across thousands of threads, and its identity comes from blockIdx and threadIdx. A launch configuration is a contract with the hardware about blocks, warps and resources. And the naive C = A + B is correct but starves the chip, because 12 bytes per FLOP puts it 240 times left of the H100 ridge, near 0.4 percent of peak compute.

If you missed the hardware side, Part 1 covers the SM, block, warp and thread hierarchy and the 120 FPS camera-pipeline anchor. The natural next question is the one this post kept deferring: if bytes are the bottleneck, why is moving them so expensive, and what can we do structurally? Part 3, GPU memory physics, DRAM versus SRAM, kernel fusion and FlashAttention answers it, showing how tiling and fusion turn the same hardware from starved to fed.

Frequently Asked Questions

What is GPU kernel engineering?

GPU kernel engineering is the practice of writing and tuning the functions that run across thousands of GPU threads. It covers thread indexing, launch configuration, memory access patterns, occupancy and branching behaviour. The goal is to deliver as much of the hardware’s published compute and bandwidth as the algorithm allows, rather than leaving most of it idle while threads wait for data from memory.

Why is vector addition memory-bound on a GPU?

Each element needs two 4-byte loads and one 4-byte store for a single addition, so 12 bytes move per FLOP. On an H100 with 3.35 TB/s of bandwidth, that caps throughput near 279 GFLOPS (derived), about 0.4 percent of the 67 TFLOPS FP32 peak. The arithmetic units finish instantly and then wait on memory.

What is the best block size for a CUDA kernel?

There is no universal answer, but a multiple of 32 between 128 and 256 is a sound starting point. The hard limit is 1,024 threads per block. Smaller blocks let more blocks reside on an SM, improving latency hiding, while larger blocks help cooperation through shared memory. Measure with a profiler on your kernel and device.

What is memory coalescing in CUDA?

Memory coalescing is the hardware merging the 32 loads or stores of a warp into as few memory transactions as possible. When consecutive threads access consecutive addresses, a warp reads 128 contiguous bytes in four 32-byte sectors with no waste. Strided or scattered access fetches many sectors and discards most bytes, wasting bandwidth by up to roughly eight times.

What is thread divergence and how do I avoid it?

Thread divergence happens when threads of one warp take different branches. The hardware runs each path in turn with other lanes masked, so cost is the sum of the paths. Avoid it by making branch conditions uniform within groups of 32 threads, sorting or bucketing data by case, or replacing short branches with predicated arithmetic.

Do I need to write custom kernels for AI workloads?

Usually not first. Libraries such as cuBLAS, cuDNN and fused attention kernels cover most standard operations, and inference engines wrap them. Custom kernels pay off for unusual shapes or when fusing steps to avoid round trips through HBM. Profile first, and write a kernel only when a measured memory-bound step justifies it.