NVIDIA GPU and CUDA Architecture Deep Dive: Why GPUs Dominate AI

CPU Architecture (latency-optimized):
┌────────────────────────────────────────────┐
│  Core 0  │  Core 1  │  Core 2  │  Core 3  │
│ (powerful)│(powerful)│(powerful)│(powerful)│
│ OOO exec │ OOO exec │ OOO exec │ OOO exec │
├──────────┴──────────┴──────────┴───────────┤
│           L3 Cache (32MB+)                 │
├────────────────────────────────────────────┤
│           Main Memory (DDR5, ~50 GB/s)     │
└────────────────────────────────────────────┘

GPU Architecture (throughput-optimized) — H100:
┌──┬──┬──┬──┬──┬──┬──┬──┬──┬──┬──┐
│SM│SM│SM│SM│SM│SM│SM│SM│SM│SM│...│  132 SMs
│  │  │  │  │  │  │  │  │  │  │   │  each SM = 128 CUDA cores
└──┴──┴──┴──┴──┴──┴──┴──┴──┴──┴───┘
       L2 Cache (50MB)
       HBM3 (80GB, 3.35 TB/s)

GPU Device (H100 SXM5)
│
├── SM 0 (Streaming Multiprocessor)
│   ├── 4x Warp Schedulers
│   ├── 4x Dispatch Units
│   ├── 128x CUDA Cores (FP32)
│   ├── 64x FP64 Cores
│   ├── 4x Tensor Cores (4th gen: FP8/FP16/BF16/TF32)
│   ├── 1x Special Function Unit (SFU: sin, cos, sqrt)
│   ├── Load/Store Units (LD/ST)
│   ├── L1 Cache + Shared Memory = 228KB (ratio configurable)
│   └── Register File (65536 x 32-bit registers)
│
├── SM 1 ...
├── SM 2 ...
│   ...
├── SM 131 ...
│
├── L2 Cache (50MB)
└── HBM3 Memory (80GB, 3.35 TB/s bandwidth)

Thread Block (example: 256 threads)
├── Warp 0: Threads  0-31  → execute same instruction in lockstep
├── Warp 1: Threads 32-63  → execute same instruction in lockstep
├── Warp 2: Threads 64-95  → execute same instruction in lockstep
├── ...
└── Warp 7: Threads 224-255

// This code causes warp divergence
__global__ void divergent_kernel(float* data, int n) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx % 2 == 0) {
        // Even-indexed threads in the warp execute this
        data[idx] = data[idx] * 2.0f;  // Step 1
    } else {
        // Odd-indexed threads in the warp execute this
        data[idx] = data[idx] + 1.0f;  // Step 2
    }
    // Result: During Step 1, odd threads are masked off (idle)
    //         During Step 2, even threads are masked off (idle)
    // Effective throughput: 50% of theoretical
}

// No divergence: all threads in a warp take the same path
__global__ void coherent_kernel(float* data, int n) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    // All threads in a block share the same blockIdx.x,
    // so all 32 threads in a warp take the same branch
    if (blockIdx.x % 2 == 0) {
        data[idx] = data[idx] * 2.0f;
    } else {
        data[idx] = data[idx] + 1.0f;
    }
}

Memory Hierarchy (fast → slow):

┌─────────────────────────────────────────────────────┐
│  Registers (32-bit, 65536 per SM)                   │
│  Latency: ~1 cycle     BW: enormous (local access)  │
│  Scope: single thread only                          │
├─────────────────────────────────────────────────────┤
│  L1 Cache / Shared Memory (228KB/SM on H100)        │
│  Latency: ~32 cycles   BW: ~19 TB/s                 │
│  Scope: all threads in a Thread Block               │
├─────────────────────────────────────────────────────┤
│  L2 Cache (50MB on H100)                            │
│  Latency: ~200 cycles  BW: ~12 TB/s                 │
│  Scope: all SMs share                               │
├─────────────────────────────────────────────────────┤
│  Global Memory HBM3 (80GB on H100)                  │
│  Latency: ~700 cycles  BW: 3.35 TB/s                │
│  Scope: all threads on the device                   │
└─────────────────────────────────────────────────────┘

// Naive matrix multiply kernel — understand WHY it's slow
__global__ void matmul_naive(float* A, float* B, float* C, int N) {
    // Each thread computes one element of the output matrix
    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++) {
            // Two global memory reads per iteration → ~700 cycles latency each
            sum += A[row * N + k] * B[k * N + col];
        }
        C[row * N + col] = sum;
    }
}

// Launch configuration
dim3 blockDim(16, 16);  // 256 threads per block
dim3 gridDim((N + 15) / 16, (N + 15) / 16);
matmul_naive<<<gridDim, blockDim>>>(A, B, C, N);

#define TILE_SIZE 16

// Tiled matrix multiply using shared memory
__global__ void matmul_tiled(float* A, float* B, float* C, int N) {
    // Declare tiles in shared memory (~32-cycle latency vs 700 for HBM)
    __shared__ float tile_A[TILE_SIZE][TILE_SIZE];
    __shared__ float tile_B[TILE_SIZE][TILE_SIZE];

    int tx = threadIdx.x, ty = threadIdx.y;
    int row = blockIdx.y * TILE_SIZE + ty;
    int col = blockIdx.x * TILE_SIZE + tx;
    float sum = 0.0f;

    // Iterate over the K dimension in tile-sized chunks
    for (int t = 0; t < N / TILE_SIZE; t++) {
        // Step 1: Each thread loads one element into shared memory
        tile_A[ty][tx] = A[row * N + t * TILE_SIZE + tx];
        tile_B[ty][tx] = B[(t * TILE_SIZE + ty) * N + col];

        // Step 2: Wait for ALL threads to finish loading
        __syncthreads();

        // Step 3: Compute partial dot product from this tile
        // All accesses here hit shared memory (~32 cycles)
        for (int k = 0; k < TILE_SIZE; k++) {
            sum += tile_A[ty][k] * tile_B[k][tx];
        }

        // Step 4: Sync before loading the next tile
        __syncthreads();
    }

    C[row * N + col] = sum;
}

Tiling principle:

Matrix A:         Matrix B:
┌──┬──┬──┬──┐  ┌──┬──┬──┬──┐
│T0│T1│T2│T3│  │T0│T1│T2│T3│
├──┼──┼──┼──┤  ├──┼──┼──┼──┤
│  │  │  │  │  │  │  │  │  │
└──┴──┴──┴──┘  └──┴──┴──┴──┘

Each Thread Block owns one 16×16 output tile of C.
It iterates over A's row tiles and B's column tiles,
loading each pair into shared memory before computing.
Result: global memory traffic reduced by TILE_SIZE factor.

// Direct Tensor Core access via WMMA API
#include <mma.h>
using namespace nvcuda::wmma;

__global__ void tensor_core_matmul(half* a_ptr, half* b_ptr,
                                    float* c_ptr, int M, int N, int K) {
    // Each warp handles one 16x16 output tile
    // fragment = a matrix tile that the Tensor Core operates on
    fragment<matrix_a, 16, 16, 16, half, row_major> a_frag;
    fragment<matrix_b, 16, 16, 16, half, col_major> b_frag;
    fragment<accumulator, 16, 16, 16, float> c_frag;

    // Initialize accumulator to zero
    fill_fragment(c_frag, 0.0f);

    // Iterate over K dimension
    for (int k = 0; k < K; k += 16) {
        // Load fragments from global memory
        load_matrix_sync(a_frag, a_ptr + /* offset */, K);
        load_matrix_sync(b_frag, b_ptr + /* offset */, N);

        // ONE instruction: 16x16x16 matrix multiply-accumulate
        // Internally dispatches to Tensor Core hardware
        mma_sync(c_frag, a_frag, b_frag, c_frag);
    }

    // Store result back to global memory
    store_matrix_sync(c_ptr + /* offset */, c_frag, N, mem_row_major);
}

CUDA Execution Hierarchy:

Grid (entire kernel launch)
├── Block (0,0): 256 threads
│   ├── Warp 0: Threads  0-31 ── execute in lockstep on SM
│   ├── Warp 1: Threads 32-63
│   └── ...
├── Block (1,0): 256 threads
├── Block (2,0): 256 threads
└── ...

Each Block is assigned to exactly one SM.
One SM can host multiple Blocks simultaneously (occupancy).

import torch

A = torch.randn(4096, 4096, device='cuda', dtype=torch.float16)
B = torch.randn(4096, 4096, device='cuda', dtype=torch.float16)

# What happens inside this one line:
# 1. cuBLAS or CUTLASS selects the optimal kernel
# 2. Execution config: e.g., Grid(256,256), Block(16,16)
# 3. Tiled GEMM using Tensor Cores
# 4. ~4096 Thread Blocks, each with 256 threads
result = torch.matmul(A, B)  # internally calls cublasSgemmEx

SM Resource Limits (H100):
- Max concurrent Thread Blocks: 32
- Max concurrent threads: 2048
- Registers: 65536 total (divided among all threads in resident warps)
- Shared Memory: 228KB (divided among resident blocks)

Occupancy = Active Warps / Maximum Possible Warps
Goal: 50%+ (needed to hide memory latency through warp switching)

import torch
import time

# 1. Use FP16 or BF16 to activate Tensor Cores
A = torch.randn(8192, 8192, device='cuda', dtype=torch.bfloat16)
B = torch.randn(8192, 8192, device='cuda', dtype=torch.bfloat16)

# 2. Warmup: first runs include kernel compilation/selection
for _ in range(5):
    _ = torch.matmul(A, B)
torch.cuda.synchronize()

# 3. Measure
start = time.perf_counter()
for _ in range(100):
    C = torch.matmul(A, B)
torch.cuda.synchronize()
end = time.perf_counter()

# 4. Compute TFLOPS
# 8192x8192 matmul: 2 * N^3 = 2 * 8192^3 = 1.1e12 FLOPs per call
flops = 2 * 8192**3 * 100
elapsed = end - start
tflops = flops / elapsed / 1e12
print(f"Achieved: {tflops:.1f} TFLOPS")
# Target: 700+ TFLOPS (70%+ of 989 TFLOPS theoretical is achievable)