AI / Transformer Kernel Programming with CUDA

Modern AI systems are not powered by neural-network mathematics alone. The real performance of a Transformer depends heavily on how those mathematical operations are mapped onto GPU hardware.

A Transformer may contain operations such as:

text
1Matrix Multiplication
2       ↓
3Bias
4       ↓
5Activation
6       ↓
7Normalization
8       ↓
9Q / K / V projections
10       ↓
11Attention
12       ↓
13MLP / SwiGLU
14       ↓
15KV Cache

At the framework level, these operations may look like simple Python or PyTorch expressions.

For example:

python
1Q = x @ Wq
2K = x @ Wk
3V = x @ Wv
4
5scores = Q @ K.transpose(-2, -1)
6scores = softmax(scores)
7
8output = scores @ V

But underneath these operations, the GPU must perform enormous amounts of:

memory loads
memory stores
floating-point operations
thread synchronization
shared-memory communication
register operations
tensor-core operations
kernel launches

This is where GPU kernel programming becomes important.

The goal of this course is to understand how Transformer operations can be implemented and optimized at the CUDA-kernel level.

1. What Is a Neural-Network Kernel?

A kernel is a GPU function executed by many CUDA threads in parallel.

A simple CUDA kernel might look like this:

cpp
1__global__ void add_vectors(
2    const float* a,
3    const float* b,
4    float* c,
5    int n
6) {
7    int i = blockIdx.x * blockDim.x + threadIdx.x;
8
9    if (i < n) {
10        c[i] = a[i] + b[i];
11    }
12}

The kernel is executed using many GPU threads:

cpp
1int threads = 256;
2int blocks = (n + threads - 1) / threads;
3
4add_vectors<<<blocks, threads>>>(a, b, c, n);

Conceptually:

text
1CPU
2 │
3 │ Launch
4 ▼
5CUDA Kernel
6 │
7 ├── Thread 0
8 ├── Thread 1
9 ├── Thread 2
10 ├── Thread 3
11 ├── ...
12 └── Thread N

Each thread processes part of the data.

This simple model is the foundation for understanding much more complicated Transformer kernels.

Module 1 — Neural-Network Kernels

2. Vector Operations

Vector operations are among the simplest GPU workloads.

For example:

text
1C[i] = A[i] + B[i]

A CUDA implementation:

cpp
1__global__ void vector_add(
2    const float* A,
3    const float* B,
4    float* C,
5    int N
6) {
7    int idx = blockIdx.x * blockDim.x + threadIdx.x;
8
9    if (idx < N) {
10        C[idx] = A[idx] + B[idx];
11    }
12}

The important idea is that each thread owns one or more elements.

text
1Thread 0 → C[0]
2Thread 1 → C[1]
3Thread 2 → C[2]
4Thread 3 → C[3]
5...

This is called element-wise parallelism.

3. Matrix Multiplication

Matrix multiplication is one of the most important operations in AI.

Given:

text
1C = A × B

where:

text
1A = M × K
2B = K × N
3C = M × N

each output element is:

text
1C[i][j] = Σ A[i][k] × B[k][j]

A basic CUDA implementation:

cpp
1__global__ void matmul(
2    const float* A,
3    const float* B,
4    float* C,
5    int M,
6    int K,
7    int N
8) {
9    int row = blockIdx.y * blockDim.y + threadIdx.y;
10    int col = blockIdx.x * blockDim.x + threadIdx.x;
11
12    if (row < M && col < N) {
13        float sum = 0.0f;
14
15        for (int k = 0; k < K; ++k) {
16            sum += A[row * K + k] *
17                   B[k * N + col];
18        }
19
20        C[row * N + col] = sum;
21    }
22}

Launch:

cpp
1dim3 block(16, 16);
2
3dim3 grid(
4    (N + block.x - 1) / block.x,
5    (M + block.y - 1) / block.y
6);
7
8matmul<<<grid, block>>>(A, B, C, M, K, N);

This is a useful learning implementation, but it is not how high-performance AI libraries normally implement GEMM.

The naive implementation repeatedly reads data from global memory.

A faster design uses tiling and shared memory.

4. Tiled Matrix Multiplication

Suppose a block needs a small portion of matrices A and B.

Instead of repeatedly loading values from global memory, threads cooperatively load tiles into shared memory.

text
1Global Memory
2      │
3      ▼
4┌───────────────┐
5│ Tile A        │
6│ Tile B        │
7└───────────────┘
8      │
9      ▼
10Shared Memory
11      │
12      ▼
13CUDA Threads
14      │
15      ▼
16Register Accumulation

Example:

cpp
1__global__ void tiled_matmul(
2    const float* A,
3    const float* B,
4    float* C,
5    int M,
6    int K,
7    int N
8) {
9    __shared__ float tileA[16][16];
10    __shared__ float tileB[16][16];
11
12    int row = blockIdx.y * 16 + threadIdx.y;
13    int col = blockIdx.x * 16 + threadIdx.x;
14
15    float sum = 0.0f;
16
17    for (int t = 0; t < (K + 15) / 16; ++t) {
18
19        int a_col = t * 16 + threadIdx.x;
20        int b_row = t * 16 + threadIdx.y;
21
22        tileA[threadIdx.y][threadIdx.x] =
23            (row < M && a_col < K)
24                ? A[row * K + a_col]
25                : 0.0f;
26
27        tileB[threadIdx.y][threadIdx.x] =
28            (b_row < K && col < N)
29                ? B[b_row * N + col]
30                : 0.0f;
31
32        __syncthreads();
33
34        for (int k = 0; k < 16; ++k) {
35            sum += tileA[threadIdx.y][k] *
36                   tileB[k][threadIdx.x];
37        }
38
39        __syncthreads();
40    }
41
42    if (row < M && col < N) {
43        C[row * N + col] = sum;
44    }
45}

The key optimization is data reuse.

Instead of:

text
1Load → Calculate → Load → Calculate

we try to do:

text
1Load tile
2    ↓
3Reuse tile many times
4    ↓
5Calculate

This principle appears repeatedly throughout optimized Transformer kernels.

5. Bias Kernels

A bias operation is usually element-wise:

text
1Y = X + bias

CUDA:

cpp
1__global__ void add_bias(
2    float* x,
3    const float* bias,
4    int rows,
5    int cols
6) {
7    int idx = blockIdx.x * blockDim.x + threadIdx.x;
8
9    int total = rows * cols;
10
11    if (idx < total) {
12        int col = idx % cols;
13        x[idx] += bias[col];
14    }
15}

Although mathematically simple, kernel fusion can make this operation more efficient.

Instead of:

text
1MatMul
2 ↓
3Global Memory
4 ↓
5Bias Kernel
6 ↓
7Global Memory

a fused kernel can potentially perform:

text
1MatMul
2 ↓
3Bias
4 ↓
5Store

This reduces intermediate memory traffic and kernel launches.

6. ReLU

ReLU is:

text
1ReLU(x) = max(0, x)

CUDA:

cpp
1__global__ void relu(
2    const float* x,
3    float* y,
4    int n
5) {
6    int i = blockIdx.x * blockDim.x + threadIdx.x;
7
8    if (i < n) {
9        y[i] = fmaxf(0.0f, x[i]);
10    }
11}

ReLU is easy to parallelize because every element is independent.

7. GELU

Transformers commonly use smoother activation functions such as GELU.

An approximate form is:

text
1GELU(x) ≈ 0.5x(1 + tanh(...))

CUDA:

cpp
1__global__ void gelu(
2    const float* x,
3    float* y,
4    int n
5) {
6    int i = blockIdx.x * blockDim.x + threadIdx.x;
7
8    if (i < n) {
9        float v = x[i];
10
11        float c =
12            0.044715f * v * v * v;
13
14        y[i] =
15            0.5f * v *
16            (1.0f + tanhf(
17                0.79788456f * (v + c)
18            ));
19    }
20}

The important lesson is that activation kernels are generally memory-bound rather than compute-heavy.

That makes memory access and fusion important optimization targets.

8. SiLU

SiLU is defined as:

text
1SiLU(x) = x × sigmoid(x)

where:

text
1sigmoid(x) = 1 / (1 + e^-x)

CUDA:

cpp
1__global__ void silu(
2    const float* x,
3    float* y,
4    int n
5) {
6    int i = blockIdx.x * blockDim.x + threadIdx.x;
7
8    if (i < n) {
9        float v = x[i];
10        float sigmoid =
11            1.0f / (1.0f + expf(-v));
12
13        y[i] = v * sigmoid;
14    }
15}

SiLU is especially important because modern Transformer architectures frequently use it as part of SwiGLU.

9. SwiGLU

SwiGLU combines a gating projection with an activation.

A simplified representation is:

text
1SwiGLU(x) = SiLU(xW_gate) × (xW_up)

Conceptually:

text
1                ┌── W_gate ── SiLU ──┐
2Input ──────────┤                    × ── Output
3                └── W_up ────────────┘

A simple CUDA element-wise kernel:

cpp
1__global__ void swiglu(
2    const float* gate,
3    const float* up,
4    float* output,
5    int n
6) {
7    int i = blockIdx.x * blockDim.x + threadIdx.x;
8
9    if (i < n) {
10        float g = gate[i];
11
12        float sigmoid =
13            1.0f / (1.0f + expf(-g));
14
15        float silu = g * sigmoid;
16
17        output[i] = silu * up[i];
18    }
19}

In production systems, the surrounding projections and activation may be fused to reduce memory traffic.

10. Softmax

Softmax converts values into normalized probabilities:

text
1softmax(xᵢ) = exp(xᵢ) / Σ exp(xⱼ)

A numerically stable implementation subtracts the maximum value:

text
1m = max(x)
2
3softmax(xᵢ) =
4    exp(xᵢ - m)
5    ----------------
6    Σ exp(xⱼ - m)

This prevents large exponentials from overflowing.

A simplified CUDA implementation:

cpp
1__global__ void softmax(
2    const float* input,
3    float* output,
4    int rows,
5    int cols
6) {
7    int row = blockIdx.x;
8
9    if (row >= rows)
10        return;
11
12    float max_value = -INFINITY;
13
14    for (int col = threadIdx.x;
15         col < cols;
16         col += blockDim.x) {
17
18        max_value =
19            fmaxf(max_value,
20                  input[row * cols + col]);
21    }
22
23    // A production implementation should
24    // perform a block-wide reduction here.
25}

The difficult part is not calculating exp().

The difficult part is performing the parallel reductions efficiently.

Softmax typically requires:

text
11. Find maximum
22. Subtract maximum
33. Calculate exponentials
44. Sum exponentials
55. Divide

That makes it a good example of why GPU programming requires synchronization and reduction algorithms.

11. LayerNorm

LayerNorm normalizes activations across a feature dimension.

Conceptually:

text
1mean = average(x)
2
3variance = average((x - mean)²)
4
5normalized =
6    (x - mean) / sqrt(variance + epsilon)
7
8output =
9    normalized × gamma + beta

A simplified CPU-style equation:

cpp
1float mean = 0.0f;
2
3for (int i = 0; i < hidden; ++i) {
4    mean += x[i];
5}
6
7mean /= hidden;

Then:

cpp
1float variance = 0.0f;
2
3for (int i = 0; i < hidden; ++i) {
4    float d = x[i] - mean;
5    variance += d * d;
6}
7
8variance /= hidden;

Finally:

cpp
1x[i] =
2    (x[i] - mean) /
3    sqrtf(variance + 1e-5f);

The GPU version must parallelize the reduction operations.

12. RMSNorm

Modern LLMs frequently use RMSNorm.

Unlike LayerNorm, RMSNorm does not require calculating the mean.

The RMS value is:

text
1RMS(x) =
2sqrt(
3    mean(x²) + ε
4)

Then:

text
1yᵢ = xᵢ / RMS(x) × weightᵢ

A simplified implementation:

cpp
1__global__ void rmsnorm(
2    const float* x,
3    const float* weight,
4    float* y,
5    int hidden,
6    float eps
7) {
8    int row = blockIdx.x;
9
10    float sum = 0.0f;
11
12    for (int i = threadIdx.x;
13         i < hidden;
14         i += blockDim.x) {
15
16        float v = x[row * hidden + i];
17        sum += v * v;
18    }
19
20    // Production implementation:
21    // reduce sum across the block.
22
23    float rms =
24        sqrtf(sum / hidden + eps);
25
26    for (int i = threadIdx.x;
27         i < hidden;
28         i += blockDim.x) {
29
30        y[row * hidden + i] =
31            (x[row * hidden + i] / rms)
32            * weight[i];
33    }
34}

This kernel illustrates a recurring pattern:

text
1Load
2 ↓
3Reduction
4 ↓
5Compute normalization factor
6 ↓
7Load/compute
8 ↓
9Store

Module 2 — Transformer Kernels

13. Q, K and V Projections

The attention mechanism begins by transforming hidden states into:

text
1Q = XWq
2
3K = XWk
4
5V = XWv

These are matrix multiplications.

text
1Input X
2 │
3 ├────────── Wq → Q
4 │
5 ├────────── Wk → K
6 │
7 └────────── Wv → V

For a Transformer:

text
1X:  [batch, sequence, hidden]
2
3Wq: [hidden, heads × head_dim]
4
5Q:  [batch, sequence, heads × head_dim]

The computational core is still GEMM.

This is why highly optimized matrix multiplication is fundamental to Transformer performance.

14. QKᵀ

Attention scores are calculated using:

text
1Scores = QKᵀ

For one attention head:

text
1Q: [sequence, head_dim]
2
3K: [sequence, head_dim]
4
5Kᵀ: [head_dim, sequence]
6
7Scores:
8[sequence, sequence]

Therefore:

text
1Attention matrix size = sequence²

This quadratic relationship is one reason long-context inference is computationally expensive.

For example:

text
1Sequence = 1024
2
3Attention matrix =
41024 × 1024

while:

text
1Sequence = 8192
2
3Attention matrix =
48192 × 8192

The number of attention-score elements grows quadratically.

15. Scaling

Attention scores are scaled:

text
1Scores =
2(QKᵀ) / sqrt(head_dim)

CUDA:

cpp
1float scale =
2    rsqrtf((float)head_dim);
3
4score *= scale;

Scaling helps maintain stable numerical behavior as the dimensionality of the query and key vectors increases.

16. Causal Masking

Autoregressive language models cannot allow a token to see future tokens.

For example:

text
1Token 1 → Token 1
2
3Token 2 → Token 1, Token 2
4
5Token 3 → Token 1, Token 2, Token 3

The attention matrix therefore looks conceptually like:

Masked positions can be assigned a very negative value:

cpp
1if (key_position > query_position) {
2    score = -INFINITY;
3}

After softmax, these positions effectively receive probability zero.

17. Attention × V

After softmax:

text
1P = softmax(QKᵀ / sqrt(d))

the output is:

text
1Output = P × V

The complete attention calculation becomes:

text
1Q
2 │
3 ├───────┐
4         ▼
5        QKᵀ
6         │
7         ▼
8      Scaling
9         │
10         ▼
11       Masking
12         │
13         ▼
14       Softmax
15         │
16         ▼
17       × V
18         │
19         ▼
20     Attention
21       Output

A naive implementation materializes the entire score matrix:

text
1Q
2 ↓
3QKᵀ
4 ↓
5Global Memory
6 ↓
7Softmax
8 ↓
9Global Memory
10 ↓
11× V
12 ↓
13Output

This is straightforward but can create substantial memory traffic.

That leads directly to attention optimization.

Module 3 — Modern Transformer Kernels

18. Rotary Positional Embeddings — RoPE

RoPE encodes positional information by rotating query and key vector pairs.

For a pair:

text
1(x₀, x₁)

the rotation can be expressed as:

text
1x₀' = x₀ cos(θ) - x₁ sin(θ)
2
3x₁' = x₀ sin(θ) + x₁ cos(θ)

CUDA-style implementation:

cpp
1__global__ void rope(
2    float* q,
3    float* k,
4    const float* cos_table,
5    const float* sin_table,
6    int n
7) {
8    int i = blockIdx.x * blockDim.x + threadIdx.x;
9
10    int pair = i * 2;
11
12    if (pair + 1 < n) {
13
14        float q0 = q[pair];
15        float q1 = q[pair + 1];
16
17        float k0 = k[pair];
18        float k1 = k[pair + 1];
19
20        float c = cos_table[i];
21        float s = sin_table[i];
22
23        q[pair]     = q0 * c - q1 * s;
24        q[pair + 1] = q0 * s + q1 * c;
25
26        k[pair]     = k0 * c - k1 * s;
27        k[pair + 1] = k0 * s + k1 * c;
28    }
29}

Real implementations must account for tensor layout, head dimensions, position indices, data types, and batching.

19. Grouped Query Attention — GQA

Traditional multi-head attention can have:

text
1Number of Q heads =
2Number of K heads =
3Number of V heads

Grouped Query Attention changes this relationship.

For example:

text
1Q heads = 32
2
3K heads = 8
4
5V heads = 8

Each KV head serves multiple query heads.

text
1Q0 ─┐
2Q1 ─┼── KV0
3Q2 ─┤
4Q3 ─┘
5
6Q4 ─┐
7Q5 ─┼── KV1
8Q6 ─┤
9Q7 ─┘

This reduces KV memory requirements compared with standard multi-head attention.

A kernel must correctly map:

text
1query_head
2      ↓
3kv_head

For example:

cpp
1int group_size = num_q_heads / num_kv_heads;
2
3int kv_head =
4    query_head / group_size;

That mapping is simple mathematically but becomes important when designing efficient memory access.

20. Multi-Query Attention — MQA

MQA takes the idea further:

text
1Many Q heads
2      ↓
3One K head
4      ↓
5One V head

For example:

text
1Q = 32 heads
2
3K = 1 head
4
5V = 1 head

This can significantly reduce KV-cache memory.

The tradeoff is that sharing K/V representations can change model quality compared with a full multi-head design.

21. KV Cache

During autoregressive generation, the model repeatedly generates one token at a time.

Without caching:

text
1Token 1
2 ↓
3Recompute everything
4
5Token 2
6 ↓
7Recompute previous information
8
9Token 3
10 ↓
11Recompute previous information

With KV caching:

text
1Previous K ──────┐
2                 ├── Attention
3Previous V ──────┤
4                 │
5Current Q ───────┘

The cache stores previous keys and values.

Conceptually:

text
1KV Cache
2
3Layer
4 ├── Head
5 │    ├── K[positions]
6 │    └── V[positions]

A simplified append kernel:

cpp
1__global__ void append_kv_cache(
2    const float* new_k,
3    const float* new_v,
4    float* cache_k,
5    float* cache_v,
6    int position,
7    int head_dim
8) {
9    int i =
10        blockIdx.x * blockDim.x +
11        threadIdx.x;
12
13    if (i < head_dim) {
14        cache_k[position * head_dim + i] =
15            new_k[i];
16
17        cache_v[position * head_dim + i] =
18            new_v[i];
19    }
20}

Real inference engines require more sophisticated cache layouts because memory capacity, batching, paging, and fragmentation become important.

22. Mixture of Experts — MoE

Mixture-of-Experts models contain multiple expert networks.

Instead of sending every token through every expert:

text
1Token
2  ↓
3Router
4  ↓
5Select experts
6  ↓
7Expert 3
8Expert 7
9  ↓
10Combine

The router determines which experts process each token.

A conceptual flow:

text
1Input Tokens
2      │
3      ▼
4    Router
5      │
6 ┌────┴────┐
7 ▼         ▼
8Expert A  Expert B
9 │         │
10 └────┬────┘
11      ▼
12   Combine

Kernel optimization becomes challenging because tokens may need to be rearranged and dispatched to different experts.

This introduces operations such as:

routing
sorting
token permutation
expert batching
GEMM
output permutation

MoE therefore combines compute optimization with memory-layout optimization.

Module 4 — Attention Optimization

23. Naive Attention

A straightforward implementation computes:

text
1S = QKᵀ
2P = softmax(S)
3O = PV

The problem is that S and P can be very large.

For sequence length N:

text
1Attention matrix = N × N

Therefore memory consumption grows quadratically.

The naive approach:

text
1Q ────────┐
2          ▼
3        QKᵀ
4          │
5          ▼
6    Store attention
7          │
8          ▼
9       Softmax
10          │
11          ▼
12    Store probabilities
13          │
14          ▼
15        × V

creates large intermediate tensors.

24. Tiled Attention

Instead of processing the entire attention matrix at once, process it in tiles.

text
1Q
2 │
3 ├── Q tile
4 │
5 ▼
6┌──────────────┐
7│ K tile       │
8├──────────────┤
9│ Compute      │
10│ scores       │
11└──────────────┘
12      │
13      ▼
14  Softmax tile
15      │
16      ▼
17   × V tile

This improves data reuse.

The kernel repeatedly loads small portions of:

text
1Q
2K
3V

into faster GPU memory.

25. Memory-Efficient Attention

The next step is avoiding unnecessary materialization.

Instead of:

text
1QKᵀ
2 ↓
3Store huge matrix
4 ↓
5Softmax
6 ↓
7Store probability matrix
8 ↓
9× V

we compute blocks and accumulate the output.

Conceptually:

text
1Q tile
2  ↓
3K tile
4  ↓
5Score tile
6  ↓
7Online softmax
8  ↓
9V tile
10  ↓
11Output accumulator

The attention matrix does not need to exist completely in global memory.

This dramatically changes the memory behavior of the algorithm.

26. FlashAttention-Style Kernels

FlashAttention-style algorithms combine:

tiling
shared memory
register accumulation
online softmax
reduced global-memory traffic
careful synchronization

The central idea is:

Compute attention in blocks while avoiding unnecessary reads and writes of the full attention matrix.

Conceptually:

text
1                 ┌──────────────┐
2Q ──────────────►│              │
3                 │ Tiled        │
4K ──────────────►│ Attention    │
5                 │              │
6V ──────────────►│              │
7                 └──────┬───────┘
8                        │
9                        ▼
10                     Output

The algorithm maintains running softmax statistics.

For a row, conceptually maintain:

text
1m = running maximum
2
3l = running normalization factor
4
5O = running output

When a new tile arrives, the previous accumulated values are rescaled and combined with the new tile.

This makes it possible to compute mathematically equivalent attention without storing the entire attention probability matrix.

27. Why Memory Matters So Much

GPU performance is not only about arithmetic throughput.

Consider:

text
1Global Memory
2     ↓
3Shared Memory
4     ↓
5Registers
6     ↓
7Arithmetic

The closer data is to the computation, the faster it can generally be accessed.

A useful mental model is:

text
1Registers
2   ↑
3Shared Memory
4   ↑
5L2 Cache
6   ↑
7Global Memory
8   ↑
9Host Memory

A poorly designed kernel might repeatedly perform:

text
1Global Load
2     ↓
3Compute
4     ↓
5Global Store
6     ↓
7Global Load
8     ↓
9Compute

An optimized kernel tries to reuse data:

text
1Global Load
2     ↓
3Shared Memory
4     ↓
5Registers
6     ↓
7Many calculations
8     ↓
9Global Store

This is why Transformer kernel optimization is fundamentally a data-movement problem as well as a computation problem.

28. Kernel Fusion

Suppose a model performs:

text
1RMSNorm
2 ↓
3RoPE
4 ↓
5Attention preparation

A naive implementation may launch multiple kernels:

text
1Kernel 1 → RMSNorm
2Kernel 2 → RoPE
3Kernel 3 → Preparation

Each kernel can require reading and writing intermediate data.

Fusion combines operations:

text
1Fused Kernel
2 ├── RMSNorm
3 ├── RoPE
4 └── Preparation

Benefits can include:

fewer kernel launches
fewer global-memory operations
better cache reuse
lower intermediate-memory usage

However, fusion is not automatically better.

A fused kernel can also:

increase register usage
reduce occupancy
make synchronization harder
increase compilation complexity

Therefore optimization requires measurement.

29. Warp-Level Programming

CUDA threads are organized into warps.

A warp commonly contains:

text
132 threads

Threads in a warp can cooperate using warp-level primitives.

For example:

cpp
1float value = ...;
2
3float sum =
4    __shfl_down_sync(
5        0xffffffff,
6        value,
7        16
8    );

Repeated shuffle operations can perform reductions without using shared memory for every step.

A conceptual warp reduction:

text
132 values
2   ↓
316 values
4   ↓
58 values
6   ↓
74 values
8   ↓
92 values
10   ↓
111 value

Warp-level programming is particularly useful for:

reductions
softmax
normalization
attention
small matrix operations

30. Shared Memory

Shared memory is visible to threads within a CUDA block.

Example:

cpp
1__shared__ float buffer[256];

Threads can load data into it:

cpp
1buffer[threadIdx.x] = input[idx];
2
3__syncthreads();

Then reuse the data:

cpp
1float value = buffer[threadIdx.x];

The synchronization point:

cpp
1__syncthreads();

ensures that threads reach the required memory state before continuing.

This is essential when multiple threads cooperate on shared data.

31. Registers

Registers are extremely fast and private to each thread.

For example:

cpp
1float accumulator = 0.0f;

may be stored in a register.

Matrix multiplication kernels often keep partial sums in registers:

text
1A tile
2 ↓
3B tile
4 ↓
5Multiply
6 ↓
7Register accumulator
8 ↓
9Repeat
10 ↓
11Store result

But there is a tradeoff.

Excessive register usage can reduce the number of active warps.

Therefore:

text
1More registers per thread
2        ↓
3Potentially fewer active warps

Optimization requires balancing register usage with occupancy and instruction-level parallelism.

32. Tensor Cores

Modern NVIDIA GPUs provide specialized hardware for matrix operations.

Tensor Cores can accelerate operations such as:

text
1FP16
2BF16
3TF32
4FP8

depending on GPU architecture and supported instructions.

This is extremely important for Transformer workloads because matrix multiplication dominates much of the compute.

Conceptually:

text
1CUDA Threads
2     ↓
3Matrix Fragments
4     ↓
5Tensor Core
6     ↓
7Matrix Accumulation

High-performance Transformer libraries therefore often use specialized matrix-multiply paths rather than implementing everything with scalar floating-point operations.

33. Data Types and Precision

Transformer kernels may operate with different numerical formats:

text
1FP32
2FP16
3BF16
4FP8
5INT8

For example:

cpp
1__half

can represent FP16 values.

Reduced precision can improve:

memory bandwidth
cache efficiency
tensor-core throughput

But numerical stability must be considered.

For example, softmax is sensitive to overflow and underflow.

A common strategy is:

text
1Input:
2FP16 / BF16
3
4Accumulation:
5FP32
6
7Output:
8FP16 / BF16

The exact strategy depends on the kernel and model.

34. Building the Transformer CUDA Kernel Library

The major project for this course is a complete CUDA kernel library.

Recommended structure:

text
1transformer-kernels/
2│
3├── CMakeLists.txt
4├── include/
5│   └── kernels.cuh
6│
7├── kernels/
8│   ├── matmul.cu
9│   ├── softmax.cu
10│   ├── rmsnorm.cu
11│   ├── rope.cu
12│   ├── silu.cu
13│   ├── swiglu.cu
14│   ├── attention.cu
15│   ├── gqa.cu
16│   ├── kv_cache.cu
17│   └── fused_transformer.cu
18│
19├── tests/
20│   ├── test_matmul.cu
21│   ├── test_softmax.cu
22│   ├── test_rmsnorm.cu
23│   └── test_attention.cu
24│
25├── benchmarks/
26│   ├── benchmark_matmul.cu
27│   ├── benchmark_attention.cu
28│   └── benchmark_kv_cache.cu
29│
30└── README.md

35. Example Kernel API

A common design is to expose simple C++ interfaces:

cpp
1void launch_matmul(
2    const float* A,
3    const float* B,
4    float* C,
5    int M,
6    int K,
7    int N
8);

Then the implementation remains inside the .cu file:

cpp
1void launch_matmul(
2    const float* A,
3    const float* B,
4    float* C,
5    int M,
6    int K,
7    int N
8) {
9    dim3 block(16, 16);
10
11    dim3 grid(
12        (N + 15) / 16,
13        (M + 15) / 16
14    );
15
16    matmul<<<grid, block>>>(
17        A, B, C, M, K, N
18    );
19}

This separation makes the library easier to test and integrate.

36. Connecting CUDA Kernels to PyTorch

A useful next step is exposing CUDA kernels through a Python interface.

Conceptually:

text
1Python
2  │
3  ▼
4PyTorch Extension
5  │
6  ▼
7C++ Binding
8  │
9  ▼
10CUDA Kernel
11  │
12  ▼
13GPU

A Python API might eventually look like:

python
1import transformer_kernels
2
3output = transformer_kernels.rmsnorm(
4    x,
5    weight
6)

The Python user does not need to know how the GPU kernel works internally.

This creates a practical bridge between:

text
1AI Framework Programming

and:

text
1GPU Kernel Programming

37. Benchmarking Kernels

A kernel should never be called "optimized" simply because the code looks complicated.

Measure it.

Important metrics include:

text
1Latency
2Throughput
3Memory bandwidth
4FLOPS
5GPU utilization
6Occupancy
7Register usage
8Shared-memory usage
9Kernel launch overhead

A simple CUDA timing pattern uses CUDA events:

cpp
1cudaEvent_t start, stop;
2
3cudaEventCreate(&start);
4cudaEventCreate(&stop);
5
6cudaEventRecord(start);
7
8my_kernel<<<grid, block>>>(...);
9
10cudaEventRecord(stop);
11cudaEventSynchronize(stop);
12
13float milliseconds = 0.0f;
14
15cudaEventElapsedTime(
16    &milliseconds,
17    start,
18    stop
19);

Then calculate throughput where appropriate.

For matrix multiplication:

text
1FLOPs ≈ 2 × M × N × K

and:

text
1Performance =
2FLOPs / execution_time

This allows different kernel implementations to be compared objectively.

38. Correctness Testing

Performance is meaningless if the kernel produces incorrect results.

A good testing strategy compares the CUDA implementation with a trusted reference.

For example:

python
1reference = torch.nn.functional.softmax(
2    x,
3    dim=-1
4)
5
6custom = transformer_kernels.softmax(x)

Then compare:

python
1torch.testing.assert_close(
2    custom,
3    reference,
4    rtol=1e-4,
5    atol=1e-5
6)

Different numerical precisions require appropriate tolerances.

Test:

text
1Small tensors
2Large tensors
3Odd dimensions
4Different batch sizes
5Different sequence lengths
6Different data types
7Boundary conditions

39. Profiling

CUDA profiling helps identify where the GPU spends time.

A kernel may be:

text
1Compute-bound

or:

text
1Memory-bound

or limited by:

text
1Latency
2Synchronization
3Register pressure
4Occupancy
5Memory access patterns

The important workflow is:

text
1Implement
2   ↓
3Validate
4   ↓
5Benchmark
6   ↓
7Profile
8   ↓
9Identify bottleneck
10   ↓
11Optimize
12   ↓
13Benchmark again

Do not optimize based only on intuition.

40. Complete Transformer Kernel Pipeline

After implementing the individual kernels, the system begins to look like a real Transformer execution engine.

text
1Input
2  │
3  ▼
4Embedding
5  │
6  ▼
7RMSNorm
8  │
9  ▼
10Q / K / V Projection
11  │
12  ▼
13RoPE
14  │
15  ▼
16GQA / MQA
17  │
18  ▼
19Attention
20  │
21  ├── QKᵀ
22  ├── Scaling
23  ├── Masking
24  ├── Softmax
25  └── Attention × V
26  │
27  ▼
28Output Projection
29  │
30  ▼
31Residual
32  │
33  ▼
34RMSNorm
35  │
36  ▼
37SwiGLU
38  │
39  ├── Gate Projection
40  ├── SiLU
41  ├── Up Projection
42  └── Element-wise Multiply
43  │
44  ▼
45Down Projection
46  │
47  ▼
48Residual
49  │
50  ▼
51Next Transformer Layer

During autoregressive inference, the KV cache participates in the attention stage:

text
1Current Q
2   │
3   ├──────────────┐
4   │              │
5   ▼              ▼
6Current K      KV Cache
7                  │
8                  ▼
9               Attention
10                  ▲
11                  │
12               Current V

41. Putting the Project Together

The final project should contain at least these kernels:

text
1kernels/
2├── matmul.cu
3├── softmax.cu
4├── rmsnorm.cu
5├── rope.cu
6├── silu.cu
7├── swiglu.cu
8├── attention.cu
9├── gqa.cu
10├── kv_cache.cu
11└── fused_transformer.cu

A useful development order is:

text
11. Vector operations
2       ↓
32. Matrix multiplication
4       ↓
53. Tiled matrix multiplication
6       ↓
74. Softmax
8       ↓
95. RMSNorm
10       ↓
116. SiLU
12       ↓
137. SwiGLU
14       ↓
158. RoPE
16       ↓
179. QKV projection
18       ↓
1910. Attention
20       ↓
2111. GQA
22       ↓
2312. KV Cache
24       ↓
2513. Kernel fusion
26       ↓
2714. Attention optimization
28       ↓
2915. Benchmarking

This progression moves from simple CUDA programming toward real Transformer inference workloads.

42. What You Should Understand After This Course

By the end of this specialization, you should understand the connection between:

text
1Transformer Mathematics
2        ↓
3Tensor Operations
4        ↓
5CUDA Kernels
6        ↓
7GPU Memory
8        ↓
9Warps
10        ↓
11Shared Memory
12        ↓
13Registers
14        ↓
15Tensor Cores
16        ↓
17Kernel Fusion
18        ↓
19Performance

You should be able to explain why a Transformer operation is expensive, identify whether it is primarily compute-bound or memory-bound, and design a CUDA implementation around the GPU's execution and memory model.

More importantly, you should understand that a high-performance AI system is not simply:

text
1Model + GPU

It is closer to:

text
1Model Architecture
2       +
3Algorithms
4       +
5Memory Layout
6       +
7CUDA Kernels
8       +
9GPU Hardware
10       +
11Kernel Scheduling
12       +
13Numerical Precision
14       +
15Profiling
16       =
17AI Inference Performance

43. Final Architecture of the Transformer CUDA Kernel Library

The completed project can evolve into:

text
1                    Python / C++ API
2                           │
3                           ▼
4                Transformer Kernel Library
5                           │
6          ┌────────────────┼────────────────┐
7          │                │                │
8          ▼                ▼                ▼
9      Math Kernels    Attention Kernels   Fusion
10          │                │                │
11          ├─ MatMul        ├─ QKᵀ           ├─ RMSNorm
12          ├─ Softmax       ├─ Softmax        ├─ RoPE
13          ├─ SiLU          ├─ GQA            └─ Attention
14          ├─ RMSNorm       ├─ MQA
15          └─ SwiGLU        └─ KV Cache
16                           │
17                           ▼
18                    CUDA Execution
19                           │
20             ┌─────────────┼─────────────┐
21             ▼             ▼             ▼
22          Registers    Shared Memory   Global Memory
23             │             │             │
24             └─────────────┼─────────────┘
25                           ▼
26                         GPU
27                           │
28                           ▼
29                  Transformer Inference

The ultimate goal is not simply to write CUDA code.

The goal is to learn how to translate Transformer mathematics into efficient GPU execution.

Once you understand that translation, concepts such as FlashAttention, fused kernels, GQA, KV caching, tensor-core GEMM, and modern LLM inference optimization become much easier to understand.

That is the core of AI / Transformer Kernel Programming.