Your First CUDA Kernel: Vector Addition End-to-End

Vector addition: $c_i = a_i + b_i$ for $i \in [0, N)$. On a CPU it's a one-line loop; on a GPU it's the canonical first kernel because every step is independent — perfect for one-thread-per-element parallelism.

The full program does six things, and that pattern repeats in every CUDA app:

1. Allocate input/output buffers on the host (malloc).
2. Allocate matching buffers on the device (cudaMalloc).
3. Copy inputs from host to device (cudaMemcpy).
4. Launch the kernel.
5. Copy the result back to the host.
6. Free both sides.

Master this six-step ritual and you've got the skeleton for everything else.

A single thread does one addition. The grid is sized so one thread exists per element.

__global__ void vecAdd(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];
}

Three things to notice:

• __global__ — kernel, called from host, runs on device.
• The 1D index formula is the same in every CUDA program.
• The if (i < n) guard matters: we'll launch slightly more threads than n so the count divides into whole blocks. Without the guard, the extra threads would write past the end of c.

The host and device have separate address spaces. A host pointer is not dereferenceable on the device, and vice versa. You need two allocations per buffer:

const int N = 1 << 20;        // ~1 million elements
size_t bytes = N * sizeof(float);

float *h_a = (float*)malloc(bytes);
float *h_b = (float*)malloc(bytes);
float *h_c = (float*)malloc(bytes);

float *d_a, *d_b, *d_c;
cudaMalloc(&d_a, bytes);
cudaMalloc(&d_b, bytes);
cudaMalloc(&d_c, bytes);

Convention: prefix host pointers with h_ and device pointers with d_. It's a tiny habit that saves you from cudaMemcpy(h_a, d_a, ...) mix-ups, which the compiler is happy to let through and the runtime is happy to segfault on.

cudaMemcpy is the bridge. Its fourth argument is the direction of the copy.

for (int i = 0; i < N; ++i) {
  h_a[i] = (float)i;
  h_b[i] = 2.0f * (float)i;
}

cudaMemcpy(d_a, h_a, bytes, cudaMemcpyHostToDevice);
cudaMemcpy(d_b, h_b, bytes, cudaMemcpyHostToDevice);

The four directions are HostToDevice, DeviceToHost, DeviceToDevice, and HostToHost. cudaMemcpy is synchronous with respect to the host — it blocks until the copy finishes. For overlapping compute and transfers you'd use cudaMemcpyAsync with streams, but for the basics this is fine.

Pick a block size (256 is a safe default — multiple of 32, balances occupancy on most GPUs). Then compute the number of blocks with ceiling division so you launch at least N threads:

const int threads = 256;
const int blocks  = (N + threads - 1) / threads;
vecAdd<<<blocks, threads>>>(d_a, d_b, d_c, N);

For N = 1,048,576 and threads = 256 you get 4096 blocks × 256 = 1,048,576 threads — exact. For N = 1,000,000 you'd get 3907 blocks × 256 = 1,000,192 threads and the last 192 threads would fall into the if (i < n) guard. The launch is asynchronous — control returns to the host immediately. You sync when you need the result.

Kernel launches don't return errors directly — they can't, they're async. You ask for them with two API calls:

cudaError_t launchErr = cudaGetLastError();   // catches launch-time errors
if (launchErr != cudaSuccess) {
  fprintf(stderr, "launch: %s\n", cudaGetErrorString(launchErr));
  exit(1);
}

cudaError_t runErr = cudaDeviceSynchronize();  // catches runtime errors
if (runErr != cudaSuccess) {
  fprintf(stderr, "runtime: %s\n", cudaGetErrorString(runErr));
  exit(1);
}

cudaGetLastError catches problems at launch (bad config, missing symbol). cudaDeviceSynchronize waits for the kernel to finish and surfaces problems that happened during execution (out-of-bounds, illegal address). In production code, wrap both in a CUDA_CHECK macro and use it everywhere.

The six-step ritual.

Putting it all together in one runnable file. Compile with nvcc vec_add.cu -o vec_add.

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

__global__ void vecAdd(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];
}

int main() {
  const int N = 1 << 20;
  size_t bytes = N * sizeof(float);

  float *h_a = (float*)malloc(bytes);
  float *h_b = (float*)malloc(bytes);
  float *h_c = (float*)malloc(bytes);
  for (int i = 0; i < N; ++i) { h_a[i] = i; h_b[i] = 2.0f * i; }

  float *d_a, *d_b, *d_c;
  cudaMalloc(&d_a, bytes);
  cudaMalloc(&d_b, bytes);
  cudaMalloc(&d_c, bytes);

  cudaMemcpy(d_a, h_a, bytes, cudaMemcpyHostToDevice);
  cudaMemcpy(d_b, h_b, bytes, cudaMemcpyHostToDevice);

  int threads = 256;
  int blocks  = (N + threads - 1) / threads;
  vecAdd<<<blocks, threads>>>(d_a, d_b, d_c, N);

  if (cudaGetLastError() != cudaSuccess) return 1;
  if (cudaDeviceSynchronize() != cudaSuccess) return 2;

  cudaMemcpy(h_c, d_c, bytes, cudaMemcpyDeviceToHost);
  printf("c[42] = %f\n", h_c[42]);  // expect 126.0

  cudaFree(d_a); cudaFree(d_b); cudaFree(d_c);
  free(h_a); free(h_b); free(h_c);
}

Vector add does one add per two loads + one store — three memory ops per FLOP. It's the textbook memory-bound kernel: even with perfect coalescing the GPU spends most of its time waiting on VRAM. Don't expect 80 TFLOPS here; expect a number close to your card's memory bandwidth divided by 12 bytes.

Also, cudaMemcpy is dwarfing the kernel. For N=1M floats:

• Kernel: ~50 µs
• 8 MB host→device: ~250 µs
• 4 MB device→host: ~125 µs

This is why nobody runs a single GPU kernel in isolation. Real CUDA code uploads once and runs hundreds of kernels back-to-back. Vector add is the right thing to write to learn the ritual, and the wrong thing to write if you want to demonstrate GPU speed.