CUDA Kernel

A CUDA kernel is the unit of work you write for the GPU: one function, launched once from the host, executed by many threads in parallel. If the CPU is a careful craftsman, the GPU is a factory floor, and the kernel is the procedure every worker runs at once.

This post is a production-oriented walkthrough of what a kernel is, how it maps onto NVIDIA’s programming model, how to launch and index work correctly, and how to think about memory and performance when you write one. Primary sources are NVIDIA’s current CUDA Programming Guide and the CUDA C++ Programming Guide.

---

1. What a kernel is (and is not)

In CUDA, device code is the code that runs on the GPU. A function that the host invokes for execution on the device is called a kernel, historically named that way in the programming model NVIDIA documents today (Programming Model).

Key properties:

Property	Kernel behavior
Declaration	Marked __global__ in CUDA C++
Return type	Always void (results go through memory)
Call site	Host (or another device API path); not a normal C++ call
Execution	One launch → many concurrent thread executions
Lifetime	Starts at launch; ends when all threads in the grid finish

A kernel is not:

• a single GPU “job” that runs once sequentially;
• a CPU thread with a different ABI;
• something you return values from like a normal function.

You launch it once; thousands (often millions) of threads each run the same code with a different identity, usually so each thread owns a different slice of data.

---

2. Host vs device: who does what

CUDA programs split work across two address spaces and instruction streams:

• Host, CPU + system memory. Allocates device buffers, copies data, launches kernels, synchronizes, reads results.
• Device, GPU + device memory. Runs kernels; threads share Device Hardware resources (SMs, caches, DRAM).

The usual lifecycle for a simple kernel:

1. Allocate device memory (cudaMalloc / modern wrappers).
2. Copy inputs host → device (cudaMemcpy or async variants).
3. Launch the kernel with an execution configuration.
4. Synchronize if needed (cudaDeviceSynchronize or stream sync).
5. Copy outputs device → host.
6. Free device memory.

Async streams and graphs complicate the schedule, but the mental model stays the same: the kernel is the device-side compute step.

---

3. Defining a kernel in CUDA C++

Kernels use the __global__ qualifier. That tells nvcc to compile the function for the device and to emit a host-callable launch stub (Intro to CUDA C++):

__global__ void saxpy(int n, float a, const float* x, float* y) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < n) {
        y[i] = a * x[i] + y[i];
    }
}

Notes that matter in practice:

• Prefer const on read-only pointers when the API allows, it documents intent and can help the compiler.
• Always guard out-of-range indices. Grid sizes are usually rounded up; leftover threads must no-op.
• Avoid large stack frames and recursion; keep per-thread state small (registers are precious).

Related qualifiers you will see next to kernels:

• __device__, callable from device code only (helpers used by kernels).
• __host__, CPU only (default for unmarked functions).
• __host__ __device__, compiled for both (common for small math helpers).

---

4. Launch configuration: <<<grid, block>>>

Launching a kernel uses CUDA’s execution configuration (the triple-chevron syntax):

saxpy<<<gridDim, blockDim>>>(n, a, d_x, d_y);

The first two arguments mean:

1. Grid dimensions, how many thread blocks.
2. Block dimensions, how many threads per block.

You may pass integers for 1D launches, or dim3 for 2D/3D (Writing CUDA Kernels):

dim3 block(16, 16);          // 256 threads per block
dim3 grid((N + 15) / 16, (N + 15) / 16);
matMul<<<grid, block>>>(d_A, d_B, d_C, N);

Optional further parameters (stream, shared memory bytes, cooperative/cluster launches) appear in extended launch APIs such as cudaLaunchKernelEx, use those when you outgrow the basic chevron form (CUDA C++ Programming Guide, Execution Configuration).

Rule of thumb for block size: multiples of the warp size (32 on current NVIDIA GPUs). Common choices: 128, 256, 512. Hard upper bound on threads per block is 1024 on current devices (check your compute capability’s limit in the Nsight Compute Occupancy Calculator / device props).

Ceiling division for covering n elements with 1D blocks of size B:

int B = 256;
int G = (n + B - 1) / B;
kernel<<<G, B>>>(...);

---

5. Thread hierarchy: grid → block → thread

When you launch a kernel, CUDA organizes threads into a hierarchy (Programming Model):

Grid (all threads for this kernel launch)
 └── Thread block 0 .. gridDim-1
      └── Thread 0 .. blockDim-1
           └── (hardware) Warp of 32 threads, scheduling unit on the SM

Important scheduling facts from the guide:

• The runtime assigns thread blocks to Streaming Multiprocessors (SMs).
• You cannot control or query which block lands on which SM.
• Blocks must be independent: any interleaving of blocks must still be correct. Do not assume block order.

That independence is why kernels scale across GPUs of different sizes: more SMs → more blocks in flight, same source code.

Warps (groups of 32 threads) are the hardware scheduling unit. Divergence inside a warp (different control flow) serializes paths, a major performance topic once kernels grow past “embarrassingly parallel.”

---

6. Built-in indices: giving every thread an identity

Inside a kernel, CUDA provides built-in variables (Writing CUDA Kernels):

Built-in	Meaning
threadIdx.{x,y,z}	Thread index within its block (0-based)
blockIdx.{x,y,z}	Block index within the grid
blockDim.{x,y,z}	Threads per block (launch config)
gridDim.{x,y,z}	Blocks in the grid (launch config)

Unspecified dimensions default to 1.

1D global index (most common)

int i = blockIdx.x * blockDim.x + threadIdx.x;

2D global indices (images, matrices)

int col = blockIdx.x * blockDim.x + threadIdx.x;
int row = blockIdx.y * blockDim.y + threadIdx.y;

Flattened 2D → 1D (row-major)

int idx = row * width + col;

Getting indexing wrong is the #1 source of silent corruption and illegal memory accesses (cudaErrorIllegalAddress). Print a few (blockIdx, threadIdx) values in a tiny debug kernel when you are unsure.

---

7. Memory a kernel can see

Kernels typically read/write global memory (device DRAM) via pointers passed from the host. Other spaces matter as you optimize:

Space	Scope	Latency / notes
Registers	Per thread	Fastest; spilled to local memory if exhausted
Local memory	Per thread	Backed by device memory; used for spills/large arrays
Shared memory	Per block	On-chip, explicitly managed; __shared__
Global memory	Grid-visible	High bandwidth, high latency; coalesce accesses
Constant / texture	Special caches	Good for read-only broadcast patterns

Synchronization inside a block:

__syncthreads();  // all threads in the block must reach this

There is no safe global barrier across the entire grid in ordinary kernels (cooperative groups / grid sync exist as specialized features). Design algorithms so blocks do not need to wait on each other mid-kernel, or split work across multiple kernel launches.

---

8. Worked example: vector addition

Canonical first kernel, each thread computes one output element:

__global__ void vectorAdd(const float* a, const float* b, float* c, int n) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < n) {
        c[i] = a[i] + b[i];
    }
}

// Host launch
int n = 1 << 20;
int threads = 256;
int blocks = (n + threads - 1) / threads;
vectorAdd<<<blocks, threads>>>(d_a, d_b, d_c, n);

Why the if (i < n)? Because blocks * threads may exceed n. Extra threads must not touch memory.

This kernel is memory-bound: one FLOP-ish of work per several bytes loaded. That is fine for learning indexing; real kernels chase higher arithmetic intensity (FLOPs per byte) so the SM’s math pipes stay busy relative to DRAM.

---

9. Worked example: matrix multiply (naive → tiled)

Matrix multiply is the “hello world” of GPU compute because it exposes the tension between global memory traffic and on-chip reuse.

9.1 Naive: one thread per output element

Each thread computes C[row][col] by streaming a row of A and a column of B from global memory:

__global__ void matmulNaive(const float* A, const float* B, float* C, int N) {
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;

    if (row < N && col < N) {
        float sum = 0.0f;
        for (int k = 0; k < N; ++k) {
            sum += A[row * N + k] * B[k * N + col];
        }
        C[row * N + col] = sum;
    }
}

Problem: for every multiply-add, the kernel reloads from global memory. On modern GPUs, peak FLOP/s ≫ peak DRAM bandwidth, so this pattern leaves math units idle. See NVIDIA’s discussion of the programming model and memory hierarchy in the CUDA C++ Programming Guide, and classic treatments in Programming Massively Parallel Processors (Kirk & Hwu).

9.2 Tiled: reuse via shared memory

Load tiles of A and B into __shared__ memory once per block phase, then compute many MACs from that on-chip data before moving to the next tile:

constexpr int TILE = 16;

__global__ void matmulTiled(const float* A, const float* B, float* C, int N) {
    __shared__ float As[TILE][TILE];
    __shared__ float Bs[TILE][TILE];

    int row = blockIdx.y * TILE + threadIdx.y;
    int col = blockIdx.x * TILE + threadIdx.x;
    float sum = 0.0f;

    for (int t = 0; t < (N + TILE - 1) / TILE; ++t) {
        int aCol = t * TILE + threadIdx.x;
        int bRow = t * TILE + threadIdx.y;

        As[threadIdx.y][threadIdx.x] =
            (row < N && aCol < N) ? A[row * N + aCol] : 0.0f;
        Bs[threadIdx.y][threadIdx.x] =
            (bRow < N && col < N) ? B[bRow * N + col] : 0.0f;

        __syncthreads();

        for (int k = 0; k < TILE; ++k) {
            sum += As[threadIdx.y][k] * Bs[k][threadIdx.x];
        }

        __syncthreads();
    }

    if (row < N && col < N) {
        C[row * N + col] = sum;
    }
}

What changed conceptually:

1. Cooperative loading, threads in a block fill shared tiles together.
2. __syncthreads(), everyone finishes the load before anyone multiplies; everyone finishes the MAC loop before the next load overwrites shared memory.
3. Higher arithmetic intensity, each global load fuels many FLOPs.

Production GEMM (cuBLAS / CUTLASS) goes much further: warp-level tiles, tensor cores, pipelined async copies (cp.async), epilogue fusion. The tiled kernel above is the pedagogical bridge.

---

10. Correctness checklist

Before you chase occupancy numbers:

1. Bounds checks on every global index.
2. Launch geometry covers the domain (ceil division).
3. No cross-block assumptions about ordering or shared state.
4. cudaGetLastError() right after launch + sync errors after cudaDeviceSynchronize / stream sync.
5. Initialize outputs if not every element is written.
6. Alignment / types, prefer vectorized loads (float4) only when addresses and counts allow.
7. Determinism, floating-point reduction order across threads is not associative; bit-identical CPU vs GPU results are not guaranteed.

For debugging: compute-sanitizer, Nsight Compute, and printf-from-device (sparingly) remain the standard toolkit (CUDA Compute Sanitizer).

---

11. Performance mindset for kernels

When a kernel is “slow,” classify the bottleneck:

Symptom	Likely cause	First moves
Low SM throughput, high DRAM traffic	Memory-bound	Tiling, coalescing, vectorized loads, reduce redundant loads
High warp stall on memory	Latency not hidden	Raise occupancy carefully; overlap; more independent work per thread
Low efficiency with branching	Warp divergence	Restructure predicates; sort/warp-uniform control flow
Register spills in SASS	Too much per-thread state	Spill less; split kernels; __launch_bounds__
Shared memory bank conflicts	Bad shared access pattern	Pad arrays; change indexing

Occupancy (active warps per SM vs maximum) is a capacity metric, not a goal by itself. Sometimes lower occupancy with more registers per thread wins. Use Nsight Compute’s roofline and memory charts rather than optimizing occupancy in isolation (Nsight Compute).

For LLM inference and other production GPU work, kernels are often hidden behind libraries (FlashAttention, cuBLASLt, custom Triton). Understanding the kernel model is still what lets you read those profiles and decide whether the bottleneck is launch overhead, memory, or math.

---

12. Mental model to keep

Host launches kernel with <<<grid, block>>>
        │
        ▼
   Grid of blocks  ──scheduled independently──►  SMs
        │
        ▼
   Threads compute index from threadIdx / blockIdx
        │
        ▼
   Load → compute → store (global / shared / registers)
        │
        ▼
   Grid completes → host observes results in device memory

A CUDA kernel is simply: same program, many identities, one launch. Everything else, tiling, warps, tensor cores, graphs, is refinement on that idea.

---

References

1. NVIDIA. CUDA Programming Guide, Programming Model.
2. NVIDIA. CUDA Programming Guide, Intro to CUDA C++.
3. NVIDIA. CUDA Programming Guide, Writing CUDA Kernels.
4. NVIDIA. CUDA C++ Programming Guide (execution configuration, memory, advanced launch).
5. NVIDIA. CUDA Refresher: The CUDA Programming Model.
6. NVIDIA. Compute Sanitizer; Nsight Compute.
7. Kirk, D. & Hwu, W. Programming Massively Parallel Processors (matrix multiply tiling chapters).
8. NVIDIA. CUDA Occupancy Calculator / Best Practices.