CUDA Programming Basics

CUDA is NVIDIA’s platform for running your own code on the GPU. You write a kernel, a function that thousands of GPU threads run at once, and launch it from ordinary C++ code on the CPU. Every program below ran on a Tesla T4 with CUDA 13.0, built with plain nvcc.

The Thread Hierarchy

A launch creates a grid of blocks, and each block holds threads. Inside the kernel, every thread reads its own position from three built-in variables, blockIdx.x, blockDim.x and threadIdx.x, and combines them into a global index:

global index = blockIdx.x * blockDim.x + threadIdx.x

The CPU side (the host) owns ordinary memory from malloc. The GPU side (the device) has its own memory, which you allocate with cudaMalloc and fill with cudaMemcpy.

Vector Addition

The classic first kernel adds two arrays element by element. Each of about a million threads handles one element, and the program checks every result on the CPU.

// vector_add.cu - add two vectors on the GPU and check every element
#include <cstdio>
#include <cstdlib>
#include <cuda_runtime.h>

__global__ void vector_add(const float* A, const float* B, float* C, int N) {
  int i = blockIdx.x * blockDim.x + threadIdx.x;  // this thread's element
  if (i < N) C[i] = A[i] + B[i];                  // the last block may have spare threads
}

int main() {
  const int N = 1 << 20;
  const 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] = (float)i; h_B[i] = 2.0f * i; }

  float *d_A, *d_B, *d_C;                         // device (GPU) memory
  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;       // round up
  vector_add<<<blocks, threads>>>(d_A, d_B, d_C, N);
  cudaMemcpy(h_C, d_C, bytes, cudaMemcpyDeviceToHost);  // waits for the kernel

  int wrong = 0;
  for (int i = 0; i < N; i++) if (h_C[i] != 3.0f * i) wrong++;
  printf("blocks=%d threads=%d -> %d threads for %d elements\n", blocks, threads, blocks * threads, N);
  printf("C[1]=%.1f  C[N-1]=%.1f  wrong=%d\n", h_C[1], h_C[N - 1], wrong);

  cudaFree(d_A); cudaFree(d_B); cudaFree(d_C);
  free(h_A); free(h_B); free(h_C);
  return 0;
}

$ nvcc -O2 -arch=sm_75 -o vector_add vector_add.cu && ./vector_add
blocks=4096 threads=256 -> 1048576 threads for 1048576 elements
C[1]=3.0  C[N-1]=3145725.0  wrong=0

N is 2^20, an exact multiple of 256, so 4096 blocks give exactly one thread per element. -arch=sm_75 builds for the T4’s compute capability 7.5. All 1,048,576 sums are right, and C[N-1] is 3 × 1,048,575.

Nothing in the program waits for the kernel explicitly. The launch returns at once, and the cudaMemcpy back to the host is what waits: it can’t copy C until the kernel that writes C has finished.

Which Thread Gets Which Element

To see the index formula at work, this kernel records, for each element, the block and thread that handled it. Ten elements with four threads per block need three blocks.

// index_map.cu - which block and thread handled each element?
#include <cstdio>
#include <cuda_runtime.h>

__global__ void record(int* who_block, int* who_thread, int N) {
  int i = blockIdx.x * blockDim.x + threadIdx.x;
  if (i < N) { who_block[i] = blockIdx.x; who_thread[i] = threadIdx.x; }
}

int main() {
  const int N = 10, threads = 4;
  int blocks = (N + threads - 1) / threads;       // 3 blocks of 4 = 12 threads for 10 elements
  int *d_b, *d_t, h_b[N], h_t[N];
  cudaMalloc(&d_b, N * sizeof(int));
  cudaMalloc(&d_t, N * sizeof(int));
  record<<<blocks, threads>>>(d_b, d_t, N);
  cudaMemcpy(h_b, d_b, sizeof(h_b), cudaMemcpyDeviceToHost);
  cudaMemcpy(h_t, d_t, sizeof(h_t), cudaMemcpyDeviceToHost);
  for (int i = 0; i < N; i++)
    printf("element %d <- block %d, thread %d  (%d*%d + %d)\n", i, h_b[i], h_t[i], h_b[i], threads, h_t[i]);
  cudaFree(d_b); cudaFree(d_t);
  return 0;
}

$ nvcc -O2 -arch=sm_75 -o index_map index_map.cu && ./index_map
element 0 <- block 0, thread 0  (0*4 + 0)
element 1 <- block 0, thread 1  (0*4 + 1)
element 2 <- block 0, thread 2  (0*4 + 2)
element 3 <- block 0, thread 3  (0*4 + 3)
element 4 <- block 1, thread 0  (1*4 + 0)
element 5 <- block 1, thread 1  (1*4 + 1)
element 6 <- block 1, thread 2  (1*4 + 2)
element 7 <- block 1, thread 3  (1*4 + 3)
element 8 <- block 2, thread 0  (2*4 + 0)
element 9 <- block 2, thread 1  (2*4 + 1)

Block 1, thread 2 is element 1 × 4 + 2 = 6. The launch started 12 threads for 10 elements, so threads 2 and 3 of block 2 would have been elements 10 and 11. The if (i < N) guard stopped them writing past the end of the arrays, which is why there are only ten lines.

Unified Memory: One Pointer for Both Sides

cudaMallocManaged returns memory that the CPU and the GPU can both use, so the copies disappear. It does not remove the need to wait for the GPU.

// managed.cu - the same vector add with unified memory: one pointer, no cudaMemcpy
#include <cstdio>
#include <cuda_runtime.h>

__global__ void vector_add(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;
  float *A, *B, *C;
  cudaMallocManaged(&A, N * sizeof(float));       // one address, usable on CPU and GPU
  cudaMallocManaged(&B, N * sizeof(float));
  cudaMallocManaged(&C, N * sizeof(float));
  for (int i = 0; i < N; i++) { A[i] = (float)i; B[i] = 2.0f * i; C[i] = -1.0f; }

  vector_add<<<(N + 255) / 256, 256>>>(A, B, C, N);
  printf("before cudaDeviceSynchronize: C[N-1] = %.1f\n", C[N - 1]);
  cudaDeviceSynchronize();                        // wait for the GPU to finish
  printf("after  cudaDeviceSynchronize: C[N-1] = %.1f\n", C[N - 1]);

  cudaFree(A); cudaFree(B); cudaFree(C);
  return 0;
}

$ nvcc -O2 -arch=sm_75 -o managed managed.cu && ./managed
before cudaDeviceSynchronize: C[N-1] = -1.0
after  cudaDeviceSynchronize: C[N-1] = 3145725.0

The first printf ran while the kernel was still in flight, and it read the -1.0 the CPU had written. That was the result on all four runs here, but it’s a race, not a guarantee: the CPU might also see the new value. After cudaDeviceSynchronize() the answer is the right one, 3,145,725. With cudaMemcpy the wait was hidden inside the copy; with unified memory you add it yourself.

Checking for Launch Errors

A block can hold at most 1024 threads on this GPU. Here the first launch asks for 2048, and the program finds out only because it asks.

// launch_check.cu - a launch that fails silently, until you ask
#include <cstdio>
#include <cuda_runtime.h>

__global__ void fill(int* out, int v) { out[blockIdx.x * blockDim.x + threadIdx.x] = v; }

int main() {
  int *d, h = 0;
  cudaMalloc(&d, 2048 * sizeof(int));
  cudaMemset(d, 0, 2048 * sizeof(int));

  fill<<<1, 2048>>>(d, 7);                        // 2048 threads in one block: over the limit
  cudaMemcpy(&h, d, sizeof(int), cudaMemcpyDeviceToHost);
  printf("no check:   out[0] = %d\n", h);
  printf("with check: %s\n", cudaGetErrorString(cudaGetLastError()));

  fill<<<2, 1024>>>(d, 7);                        // the same 2048 threads, a legal shape
  cudaError_t err = cudaGetLastError();           // errors in the launch itself
  if (err == cudaSuccess) err = cudaDeviceSynchronize();  // errors while the kernel ran
  cudaMemcpy(&h, d, sizeof(int), cudaMemcpyDeviceToHost);
  printf("fixed:      %s, out[0] = %d\n", cudaGetErrorString(err), h);
  cudaFree(d);
  return 0;
}

$ nvcc -O2 -arch=sm_75 -o launch_check launch_check.cu && ./launch_check
no check:   out[0] = 0
with check: invalid argument
fixed:      no error, out[0] = 7

The bad launch never ran, so out[0] is still the 0 from cudaMemset, and the program carried on as if nothing had happened. cudaGetLastError() reports the problem as invalid argument. The fix keeps the same 2048 threads but splits them into two blocks of 1024.

Checking has two parts. cudaGetLastError() right after the launch catches a launch that couldn’t start, such as a wrong shape. cudaDeviceSynchronize() waits and then reports errors raised while the kernel ran, such as an illegal memory access. Real projects wrap both in a small macro.

A 2D Grid for Images and Matrices

Image-like data maps naturally onto two dimensions: x for columns and y for rows. Both the grid and the blocks can be 2D.

// add_2d.cu - a 2D grid of 2D blocks over a 1920 x 1080 image
#include <cstdio>
#include <cstdlib>
#include <cuda_runtime.h>

__global__ void add_2d(const float* A, const float* B, float* C, int W, int H) {
  int x = blockIdx.x * blockDim.x + threadIdx.x;  // column
  int y = blockIdx.y * blockDim.y + threadIdx.y;  // row
  if (x < W && y < H) { int idx = y * W + x; C[idx] = A[idx] + B[idx]; }
}

int main() {
  const int W = 1920, H = 1080;
  const size_t bytes = (size_t)W * H * sizeof(float);
  float* h_A = (float*)malloc(bytes);
  float* h_B = (float*)malloc(bytes);
  float* h_C = (float*)malloc(bytes);
  for (int y = 0; y < H; y++)
    for (int x = 0; x < W; x++) { h_A[y * W + x] = (float)x; h_B[y * W + x] = 10000.0f * y; }

  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);

  dim3 threads(16, 16);                           // 256 threads per block
  dim3 blocks((W + threads.x - 1) / threads.x, (H + threads.y - 1) / threads.y);
  add_2d<<<blocks, threads>>>(d_A, d_B, d_C, W, H);
  cudaMemcpy(h_C, d_C, bytes, cudaMemcpyDeviceToHost);

  printf("grid %d x %d blocks of %d x %d threads\n", blocks.x, blocks.y, threads.x, threads.y);
  printf("C(x=7, y=3) = %.0f   C(x=1919, y=1079) = %.0f\n", h_C[3 * W + 7], h_C[1079 * W + 1919]);
  cudaFree(d_A); cudaFree(d_B); cudaFree(d_C);
  free(h_A); free(h_B); free(h_C);
  return 0;
}

$ nvcc -O2 -arch=sm_75 -o add_2d add_2d.cu && ./add_2d
grid 120 x 68 blocks of 16 x 16 threads
C(x=7, y=3) = 30007   C(x=1919, y=1079) = 10791919

1920 / 16 is exactly 120, but 1080 / 16 is 67.5, so the grid rounds up to 68 rows of blocks. The bottom row of blocks has 8 spare rows of threads, and the x < W && y < H guard keeps them from writing. Each pixel holds x + 10,000 × y, so C(7, 3) is 30,007 and the last pixel is 1919 + 10,790,000.

Common Gotchas

• Out-of-bounds threads. Rounding the block count up is right, but then you need the if (i < N) guard. Practice 2 shows that leaving it out can look harmless.
• Reading results too early. Kernel launches are asynchronous. cudaMemcpy back to the host waits for you; unified memory and timing code don’t.
• Silent launch failures. A launch returns nothing. Check cudaGetLastError() after it and cudaDeviceSynchronize() before you trust the results.
• Block size limits. 1024 threads per block is the maximum on this GPU. 256 is a common default that fits well in every SM.

Practice

1. Change N in vector_add.cu to 1,000,000, which isn’t a multiple of 256. How many threads does the launch start, and are all the results still right?

$ nvcc -O2 -arch=sm_75 -o p1_vector_add p1_vector_add.cu && ./p1_vector_add
blocks=3907 threads=256 -> 1000192 threads for 1000000 elements
C[1]=3.0  C[N-1]=2999997.0  wrong=0

2. Remove the guard from the kernel (C[i] = A[i] + B[i]; with no if) and keep N = 1000000. Does the program still print wrong=0? Then run it under compute-sanitizer.

$ nvcc -O2 -arch=sm_75 -o p2_noguard p2_noguard.cu && ./p2_noguard
blocks=3907 threads=256 -> 1000192 threads for 1000000 elements
C[1]=3.0  C[N-1]=2999997.0  wrong=0
$ compute-sanitizer ./p2_noguard | grep -E -m 1 -o 'Invalid __global__ [a-z]+ of size [0-9]+ bytes'; compute-sanitizer ./p2_noguard | grep -E '^C\['
Invalid __global__ read of size 4 bytes
C[1]=0.0  C[N-1]=0.0  wrong=999999

3. In index_map.cu, use 3 threads per block instead of 4. Predict the block count and the last line before you run it.

$ nvcc -O2 -arch=sm_75 -o p3_index_map p3_index_map.cu && ./p3_index_map
element 0 <- block 0, thread 0  (0*3 + 0)
element 1 <- block 0, thread 1  (0*3 + 1)
element 2 <- block 0, thread 2  (0*3 + 2)
element 3 <- block 1, thread 0  (1*3 + 0)
element 4 <- block 1, thread 1  (1*3 + 1)
element 5 <- block 1, thread 2  (1*3 + 2)
element 6 <- block 2, thread 0  (2*3 + 0)
element 7 <- block 2, thread 1  (2*3 + 1)
element 8 <- block 2, thread 2  (2*3 + 2)
element 9 <- block 3, thread 0  (3*3 + 0)

4. Find the block-size limit by experiment: launch an empty kernel with 512, 1024, 1025 and 2048 threads per block, check each launch, and compare with what the device reports.

// p4_limit.cu - find the largest legal block size by launching and checking
#include <cstdio>
#include <cuda_runtime.h>

__global__ void noop() {}

int main() {
  for (int t : {512, 1024, 1025, 2048}) {
    noop<<<1, t>>>();
    cudaError_t err = cudaGetLastError();
    if (err == cudaSuccess) err = cudaDeviceSynchronize();
    printf("%4d threads per block: %s\n", t, cudaGetErrorString(err));
  }
  cudaDeviceProp p;
  cudaGetDeviceProperties(&p, 0);
  printf("maxThreadsPerBlock reported by the device: %d\n", p.maxThreadsPerBlock);
  return 0;
}

$ nvcc -O2 -arch=sm_75 -o p4_limit p4_limit.cu && ./p4_limit
 512 threads per block: no error
1024 threads per block: no error
1025 threads per block: invalid argument
2048 threads per block: invalid argument
maxThreadsPerBlock reported by the device: 1024