The CUDA Model: Host and Device

CUDA is C with a few extras that let the CPU hand work to the GPU. Two worlds, connected by a bus:

saxpy<<< blocks, threadsPerBlock >>>(n, a, d_x, d_y);
//     ^^^^^^^^^^^^^^^^^^^^^^^^^^^^  the launch configuration
//     "run this kernel on (blocks x threadsPerBlock) threads"

Thread → Block → Grid

Threads are organized in a three-level hierarchy. Every thread knows its own coordinates and uses them to find its slice of the data:

GRID  (all the threads a kernel launches)
 +-----------------------------------------------+
 |  BLOCK 0            BLOCK 1            BLOCK 2  |
 | [t0 t1 t2 t3]      [t0 t1 t2 t3]      [t0 ...] |   threadIdx.x = 0..3
 +-----------------------------------------------+
   blockIdx.x=0        blockIdx.x=1       blockIdx.x=2      blockDim.x = 4

 global index of a thread:
   i = blockIdx.x * blockDim.x + threadIdx.x
   e.g. block 2, thread 1  ->  2 * 4 + 1 = 9

The demo 04-cuda-threads.html lets you click a thread and watch this exact formula light up.

Memory: Two Separate Worlds

Every CUDA program follows the same five steps

1. cudaMalloc(&d_x, bytes);                  // allocate on the DEVICE
2. cudaMemcpy(d_x, h_x, bytes, H2D);         // copy inputs  HOST -> DEVICE
3. kernel<<<blocks, threads>>>(d_x, ...);      // launch: one thread per element
4. cudaMemcpy(h_y, d_y, bytes, D2H);         // copy result  DEVICE -> HOST
5. cudaFree(d_x);                            // free device memory

Vector Add / SAXPY — CUDA's "Hello World"

y[i] = a*x[i] + y[i] for a million elements. Every element is independent, so we give one thread to one element — the CPU's loop disappears:

__global__ void saxpy(int n, float a, const float *x, float *y) {
    // this thread's unique global index
    int i = blockIdx.x * blockDim.x + threadIdx.x;

    if (i < n) {              // bounds check (next slide)
        y[i] = a * x[i] + y[i];
    }
}

// launch: enough blocks of 256 threads to cover all n elements
int threads = 256;
int blocks  = (n + threads - 1) / threads;   // round UP
saxpy<<<blocks, threads>>>(n, 2.0f, d_x, d_y);

Worked Example: Which Thread Gets Which Element?

This one line is the whole of CUDA indexing, and it's where nearly every first bug lives. N = 1000, block size 256.

__global__ void vecAdd(float *a, float *b, float *c, int n) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < n)                      // <-- do NOT skip this
        c[i] = a[i] + b[i];
}

int blocks = (n + 255) / 256;       // = (1000+255)/256 = 4  (round UP)
vecAdd<<<blocks, 256>>>(a, b, c, n);   // launches 4 x 256 = 1024 threads

Why if (i < n) Matters

Threads come in whole blocks. If n = 1000 and a block holds 256 threads, you need 4 blocks — but 4 × 256 = 1024 threads for 1000 elements. The extra 24 must do nothing.

blocks = (1000 + 256 - 1) / 256 = 4      // ceiling division, rounds UP
threads launched = 4 * 256 = 1024
elements         = 1000
                   ^^^^^^^^  threads 1000..1023 run off the end!

    if (i < n)   // <-  this skips them. Without it, they write
                 //     past the array and corrupt GPU memory.

Watch: CUDA in a Nutshell

A fast, high-energy tour of what CUDA is and why the GPU changed computing. Great warm-up before the code.

Inside the GPU: Three Kinds of Memory

Getting speed out of a GPU is mostly about which memory you use. Closer = faster = smaller.

Worked Example: Coalescing Is Worth 30×

A warp's 32 threads issue their memory request at the same instant. If those 32 addresses are consecutive, the hardware fetches them in one transaction. If they're scattered, it needs 32.

c[i] = a[i] * 2;

warp reads a[0..31]
= 128 consecutive bytes
= ONE transaction
100% of it used

c[i] = a[i*32] * 2;

warp reads a[0], a[32],
a[64] ... a[992]
= 32 transactions
3% of each one used

// thread per ROW — feels natural
for (j=0;j<N;j++)
    sum += a[tid*N + j];

// thread per COLUMN
for (j=0;j<N;j++)
    sum += a[j*N + tid];

Lab 5 — Parallel Reduction on the GPU

"Sum a million numbers into one" is not embarrassingly parallel — the threads must cooperate. The classic answer is the tree reduction: halve the active threads each step.

step 0:  [3][1][7][0][4][1][6][3]     8 values, 4 threads active
step 1:  [7][2][13][3] . . . .        add partner, halve  (stride 4)
step 2:  [20][5] . . . . . . .        halve again         (stride 2)
step 3:  [25] . . . . . . . .         one value = block's sum (stride 1)

  for (int stride = blockDim.x/2; stride > 0; stride >>= 1) {
      if (tid < stride) sdata[tid] += sdata[tid + stride];
      __syncthreads();   // all threads finish this level first
  }

Worked Example: Trace the Reduction Tree

Sum [3, 1, 7, 0, 4, 1, 6, 3] — total 25 — with 8 threads. Serially that's 7 additions in 7 steps. In parallel it's 7 additions in 3.

for (s = blockDim.x/2; s > 0; s >>= 1) {
    if (tid < s) data[tid] += data[tid + s];
    __syncthreads();                          // everyone finishes the step
}

start   [ 3   1   7   0 | 4   1   6   3 ]
                          |   |   |   |
s=4       t0: 3+4 = 7     |   |   |   |
          t1: 1+1 = 2 <---+   |   |   |
          t2: 7+6 = 13 <------+---+   |
          t3: 0+3 = 3  <--------------+
        [ 7   2  13   3 | .   .   .   . ]    4 threads active

s=2       t0: 7 + 13 = 20
          t1: 2 +  3 =  5
        [20   5 | .   . | .   .   .   . ]    2 threads active

s=1       t0: 20 + 5 = 25
        [25 | . | .   . | .   .   .   . ]    1 thread active

DONE in 3 = log2(8) steps.

// the "obvious" version
for (s=1; s<blockDim.x; s*=2) {
  if (tid % (2*s) == 0)
    data[tid] += data[tid+s];
  __syncthreads();
}

// the version above
for (s=blockDim.x/2; s>0; s>>=1) {
  if (tid < s)
    data[tid] += data[tid+s];
  __syncthreads();
}

Lab 3 — Image Processing: One Thread per Pixel

Every output pixel is computed the same way and is independent — the ultimate data-parallel shard. A 512×512 image becomes 262,144 threads. Because images are 2D, we use a 2D grid:

__global__ void to_grayscale(const unsigned char *rgb,
                             unsigned char *gray, int w, int h) {
    int col = blockIdx.x * blockDim.x + threadIdx.x;   // x direction
    int row = blockIdx.y * blockDim.y + threadIdx.y;   // y direction

    if (col < w && row < h) {
        int p = row * w + col;
        gray[p] = 0.299f*rgb[p*3] + 0.587f*rgb[p*3+1] + 0.114f*rgb[p*3+2];
    }
}
dim3 threads(16, 16);                     // 16x16 = 256 threads per block
dim3 blocks((w+15)/16, (h+15)/16);        // enough blocks to cover the image

The Modern GPU Stack: From Kernel to Data Center

You just wrote a kernel by hand. In production, layers on top of CUDA turn one GPU into thousands serving AI at scale:

Unit IV Wrap-Up: Three Tools, One Idea

Unit IV gave you the three programming models that run essentially all parallel code today:

Recap & What's Next

Key Takeaways

• Host vs device: the CPU orchestrates, the GPU runs kernels across thousands of threads.
• Thread → block → grid: every thread finds its data with i = blockIdx.x*blockDim.x + threadIdx.x, guarded by if (i < n).
• Five steps: malloc → copy up → launch → copy down → free. Data movement is the cost.
• Reduction (Lab 5) cooperates via shared memory; image processing (Lab 3) is one thread per pixel.
• The production stack — Kubernetes, TensorRT, CUDA-X — is built on exactly these kernels.