A CUDA Kernel from First Principles

A CUDA kernel is a loop with a coordinate system.

That sentence is not a metaphor. The hardware runs many copies of the same function, and each copy needs to know which part of the data it owns. Everything else — blocks, warps, shared memory, occupancy — follows from how you answer that question.

This note builds two kernels from scratch: a vector addition, which is the smallest useful kernel, and a block reduction, which introduces shared memory and synchronisation. The code is deliberately plain. The point is to see the machine, not to hide it behind a library.

The smallest useful kernel#

The task: given arrays a, b, and c of n floats, compute c[i] = a[i] + b[i].

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

Four things are happening in that index expression.

• threadIdx.x is the thread's position inside its block.
• blockDim.x is the number of threads in the block.
• blockIdx.x is the block's position in the grid.
• Their combination gives a unique global index for every thread.

The if (i < n) is not defensive programming. It is required. If n is not a multiple of the block size, the last block has more threads than elements. Without the check, those threads write past the end of the array.

The host side is ordinary CUDA:

const int n = 1 << 24;                  // 16,777,216 floats
const size_t bytes = n * sizeof(float);
 
float *h_a = new float[n];
float *h_b = new float[n];
float *h_c = new float[n];
// ... fill h_a and h_b ...
 
float *d_a, *d_b, *d_c;
cudaMalloc(reinterpret_cast<void**>(&d_a), bytes);
cudaMalloc(reinterpret_cast<void**>(&d_b), bytes);
cudaMalloc(reinterpret_cast<void**>(&d_c), bytes);
 
cudaMemcpy(d_a, h_a, bytes, cudaMemcpyHostToDevice);
cudaMemcpy(d_b, h_b, bytes, cudaMemcpyHostToDevice);
 
const int threads = 256;
const int blocks = (n + threads - 1) / threads;
vector_add<<<blocks, threads>>>(d_a, d_b, d_c, n);
 
cudaMemcpy(h_c, d_c, bytes, cudaMemcpyDeviceToHost);

The launch configuration <<<blocks, threads>>> is a grid of blocks blocks, each with threads threads. On current NVIDIA hardware, a block can have up to 1024 threads. A warp is 32 threads, and the scheduler issues instructions per warp. A block size of 256 gives eight warps per block, which is a reasonable default for a memory-bound kernel.

The grid-stride loop#

The version above assumes one thread per element. That is fine for a vector addition, but it ties the grid size to n. For very large arrays, or for kernels that do more work per thread, a grid-stride loop is more flexible:

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

Now the grid can be smaller than the data. Each thread handles several elements, striding by the total number of threads in the grid. This is the standard pattern for kernels that need to saturate a GPU without launching millions of blocks.

A common launch heuristic is to cap the grid at a multiple of the number of streaming multiprocessors:

int device;
cudaGetDevice(&device);
int sms;
cudaDeviceGetAttribute(&sms, cudaDevAttrMultiProcessorCount, device);
 
const int threads = 256;
const int blocks = std::min((n + threads - 1) / threads, sms * 32);
vector_add_grid_stride<<<blocks, threads>>>(d_a, d_b, d_c, n);

The exact multiple is workload-dependent. The important part is that the grid is a scheduling decision, not a direct reflection of the array length.

Memory is the real boss#

Vector addition performs one floating-point addition for every twelve bytes moved: four bytes for a[i], four for b[i], four for c[i]. Its arithmetic intensity is about

A GPU with 1 TB/s of memory bandwidth and 30 TFLOP/s of peak compute has a machine balance of about 30 flop/byte. Vector addition is memory-bound by a wide margin. Making the arithmetic faster will not help; moving fewer bytes, or moving them more efficiently, will.

The first efficiency rule is coalescing. Consecutive threads should access consecutive addresses. In the kernel above, thread 0 reads a[0], thread 1 reads a[1], and so on. A warp of 32 threads touches a contiguous 128-byte region, which the memory system can service in a small number of transactions.

The opposite pattern is a stride:

// Don't do this for a simple elementwise kernel.
int i = (threadIdx.x * 32 + blockIdx.x);
c[i] = a[i] + b[i];

Here neighbouring threads are 32 elements apart. Each thread touches a different cache line, and the warp issues many more memory transactions for the same amount of useful data. The kernel becomes latency-bound. The fix is usually to change the indexing, not to add more threads.

Two small qualifiers also help the compiler: const promises that the input pointers are not written, and __restrict__ promises that they do not alias the output. Both can enable better load scheduling. They do not change the algorithm, but they give the optimizer permission to do its job.

A block reduction#

Vector addition has no communication between threads. A reduction does. The goal is to sum an array of n floats.

The simplest correct approach is a two-stage reduction. Each block reduces its slice into one partial sum, then a second kernel (or a single block) reduces the partials.

__global__ void reduce_sum(const float* __restrict__ x,
                           float* __restrict__ partial,
                           int n) {
    extern __shared__ float s[];
 
    int tid = threadIdx.x;
    int i = blockIdx.x * blockDim.x + tid;
 
    s[tid] = (i < n) ? x[i] : 0.0f;
    __syncthreads();
 
    for (int stride = blockDim.x / 2; stride > 0; stride >>= 1) {
        if (tid < stride) {
            s[tid] += s[tid + stride];
        }
        __syncthreads();
    }
 
    if (tid == 0) {
        partial[blockIdx.x] = s[0];
    }
}

There are three details worth slowing down for.

Dynamic shared memory. extern __shared__ float s[] declares a shared array whose size is supplied at launch time. The launch looks like this:

const int threads = 256;
const int blocks = (n + threads - 1) / threads;
 
float* d_partial;
cudaMalloc(reinterpret_cast<void**>(&d_partial), blocks * sizeof(float));
 
reduce_sum<<<blocks, threads, threads * sizeof(float)>>>(d_x, d_partial, n);

The third launch parameter is the number of bytes of dynamic shared memory per block. Here it is 256 * 4 = 1024 bytes, which is tiny compared with the 48–228 KB available per SM on modern GPUs.

Synchronisation. __syncthreads() is a barrier. Every thread in the block must reach it before any thread continues. The first barrier ensures that every element has been written to shared memory before the reduction reads it. The barrier inside the loop ensures that one step's writes are visible to the next step's reads.

Removing either barrier produces a race. The kernel may still pass a small test, because the timing happens to work out. That is not correctness; it is luck.

The reduction tree. On the first iteration, thread 0 adds s[128], thread 1 adds s[129], and so on. On the next iteration, the active range halves. After eight iterations for a 256-thread block, thread 0 holds the sum of the block.

The if (tid < stride) means half the threads become idle at each step. For a 256-thread block that is fine. For larger blocks, or for the final warp, a shuffle-based reduction avoids shared memory entirely:

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

__shfl_down_sync moves a value from a higher lane to a lower lane within the warp. No shared memory, no block barrier. The standard pattern is to reduce within each warp, write one value per warp to shared memory, then let the first warp reduce those values.

What to measure#

A kernel that looks efficient can still be slow for reasons the source code does not show: register pressure, low occupancy, uncoalesced accesses, or a launch configuration that leaves SMs idle.

Measure effective bandwidth first:

For vector addition, that is 3 * n * sizeof(float) divided by the kernel time. Compare it with the device's theoretical peak. If the kernel reaches a large fraction of peak, it is memory-bound and further arithmetic changes will not matter. If it does not, look at the access pattern.

Nsight Compute gives the details:

ncu --set full ./vector_add

The interesting sections are memory throughput, achieved occupancy, warp stall reasons, and the roofline chart. Nsight Systems is for the timeline: kernel launches, memory copies, and gaps between them.

Two common mistakes are worth naming. The first is timing the kernel with cudaMemcpy included; the copy can dominate. The second is timing a single launch on a small array, where launch overhead is larger than the work. Use a large enough problem and repeat the launch.

A checklist#

• Check every CUDA call. cudaGetLastError() after a launch is the minimum; a small CUDA_CHECK macro is better.
• Use cudaDeviceSynchronize() while debugging. Remove it from the hot path once the kernel is correct.
• Prefer coalesced access. Consecutive threads should touch consecutive addresses.
• Use const and __restrict__ where they are true.
• Size the grid for the machine, not for the array. A grid-stride loop is usually the right default.
• For reductions, use shared memory for the block and a second pass for the partials. Use warp shuffles for the final warp.
• Profile before optimizing. The bottleneck is often not where the source code suggests.

A kernel is a small program with a large machine underneath it. The syntax is easy. The mental model — threads, warps, blocks, memory transactions — is the part worth building.