GPU & CUDA

CUDA organizes parallel execution into a three-level hierarchy. Threads are the finest unit; each runs one instance of a kernel function. Threads group into Blocks (1D, 2D, or 3D arrays up to 1024 threads per block). Blocks group into a Grid. On hardware, groups of 32 threads form a Warp — the true unit of GPU scheduling. A Streaming Multiprocessor (SM) executes multiple warps concurrently (up to thousands of threads). Blocks are assigned to SMs at launch; all threads in a block run on the same SM and can share its shared memory and synchronize with __syncthreads(). Choosing block size as a multiple of 32 ensures full warp utilization with no idle threads.

# CUDA kernel launch
# Grid of 256 blocks, each with 128 threads (4 warps)
my_kernel<<<256, 128>>>(d_input, d_output, N)

Warp divergence occurs when threads within the same warp (32 threads) take different execution paths due to conditional branches (if/else, switch). Since all 32 threads in a warp execute in lock-step (Single Instruction Multiple Threads — SIMT), both branches must be executed serially with inactive threads masked. If half the warp takes the if-branch and half takes the else-branch, effective GPU throughput halves. Severe divergence in tight loops dramatically reduces utilization. Minimize divergence by: sorting input data so similar cases are adjacent (coalescing branch decisions within warps), restructuring algorithms to use branchless code (predicated instructions), or using separate kernels for different cases. Modern GPUs (Volta+) have Independent Thread Scheduling that reduces some divergence costs.

// Avoid divergence
// Bad: half warp takes if, half takes else
if (threadIdx.x % 2 == 0) { ... } else { ... }
// Better: separate kernels or branchless arithmetic

Memory coalescing occurs when threads in a warp access consecutive memory addresses simultaneously, allowing the GPU memory controller to service all requests in a single or few memory transactions. If threads access scattered non-consecutive addresses, the controller issues multiple separate transactions — a coalescing failure that can reduce memory bandwidth utilization by 32x. For optimal coalescing, thread N in a warp should access element N of a contiguous array (stride-1 access). Common pitfalls: accessing a 2D array column-by-column (stride = num_columns) rather than row-by-row, or using struct-of-arrays incorrectly. Structure data as Structure of Arrays (SoA) rather than Array of Structures (AoS) for coalesced access in compute kernels.

// Coalesced: thread i accesses data[i]
int idx = threadIdx.x + blockIdx.x * blockDim.x;
float val = data[idx];  // stride-1: coalesced

// Non-coalesced: thread i accesses col i of row-major matrix
float val = matrix[row * N + threadIdx.x];  // ok if row is fixed

Shared memory is a fast on-chip SRAM (typically 48-164 KB per SM) shared among all threads in a block, with latency ~100x lower than global memory. Use it to cache frequently accessed data — each thread loads one element from global memory to shared memory, then all threads in the block read from the shared cache. In matrix multiplication (tiled GEMM), each block loads 16x16 or 32x32 tiles of A and B matrices into shared memory, then computes partial dot products without additional global memory reads. This reduces global memory bandwidth consumption by a factor equal to the tile size. Declare with __shared__ float tile[32][32]. Use __syncthreads() to ensure all threads have finished loading before computing.

__global__ void matmul(float* A, float* B, float* C, int N) {
    __shared__ float tileA[16][16];
    __shared__ float tileB[16][16];
    // load tiles, __syncthreads(), compute partial sum
}

Occupancy is the ratio of active warps to the maximum possible warps on a Streaming Multiprocessor. Higher occupancy hides memory latency by switching to other ready warps while some wait for data. Occupancy is limited by three SM resources: registers per thread (register file shared across warps), shared memory per block, and max threads per SM. Reduce register usage by setting --maxrregcount or using launch_bounds to guide the compiler. Reduce shared memory if possible. Choose block sizes that are multiples of 32 and that allow 2-4 active blocks per SM. Use the CUDA Occupancy Calculator or cudaOccupancyMaxActiveBlocksPerMultiprocessor() API. Note: 100% occupancy is not always optimal — sometimes fewer warps with more registers perform better than high-occupancy low-register kernels.

import pycuda.autoinit
from pycuda.compiler import SourceModule
# Use NVIDIA Nsight Compute for occupancy analysis
# Or: cuda.Device(0).get_attributes() for SM limits

A CUDA stream is a sequence of operations (kernel launches, memory copies) executed in order on the GPU, but different streams can overlap. The default stream (stream 0) serializes all operations. Create non-default streams with cudaStreamCreate to enable concurrency. Overlap strategies: launch compute kernels on one stream while transferring data on another stream using cudaMemcpyAsync. Multiple kernels can execute concurrently if the GPU has free SMs (requires compute capability 2.0+, Hyper-Q for up to 32 concurrent kernels on Maxwell+). Use CUDA events (cudaEventRecord, cudaStreamWaitEvent) for inter-stream synchronization. Stream-based overlapping is essential for maximizing GPU utilization during inference — pre-fetch next batch while current batch computes.

cudaStream_t s1, s2;
cudaStreamCreate(&s1);
cudaStreamCreate(&s2);
// Overlap: transfer on s1, compute on s2
cudaMemcpyAsync(d_in, h_in, bytes, cudaMemcpyHostToDevice, s1);
my_kernel<<<grid, block, 0, s2>>>(d_in_prev, d_out);

Pinned memory (page-locked host memory) is allocated with cudaMallocHost() and is guaranteed to never be swapped to disk by the OS virtual memory system. DMA (Direct Memory Access) engines on the GPU can directly transfer data to/from pinned memory at peak PCI-e bandwidth without the CPU copying data through an intermediary staging buffer. Pageable host memory requires an extra intermediate copy to a temporary page-locked buffer before DMA transfer, halving effective bandwidth. Pinned memory is also required for cudaMemcpyAsync to be truly asynchronous — without pinned memory, async copies may block. Downside: pinned memory is a scarce system resource; over-allocating it starves the OS and other processes. Use for data-loading buffers, not for general application memory.

float* h_data;
cudaMallocHost(&h_data, N * sizeof(float));  // pinned
// Now async copy is truly async:
cudaMemcpyAsync(d_data, h_data, bytes, cudaMemcpyHostToDevice, stream);

Unified Memory (cudaMallocManaged) creates a single memory pool accessible by both CPU and GPU. The CUDA runtime and GPU hardware automatically migrate pages between host and device as needed, managed by the unified memory driver. Benefits: simplified programming (no explicit cudaMemcpy), ability to process datasets larger than GPU memory, and deep pointer support for complex data structures. Drawbacks: page migration overhead (first-access faults trigger page transfers) can be severe for random access patterns; memory may migrate repeatedly if CPU and GPU access the same data alternately. Best for: datasets that don't fit in GPU memory (oversubscription), complex pointer-based structures, or rapid prototyping. Optimize with cudaMemAdvise() (set preferred location, read-mostly flags) and cudaMemPrefetchAsync() to pre-migrate pages.

float* data;
cudaMallocManaged(&data, N * sizeof(float));
// Pre-migrate to GPU before kernel
cudaMemPrefetchAsync(data, N*sizeof(float), device_id, stream);
my_kernel<<<grid, block>>>(data, N);

Tensor Cores are specialized hardware units introduced in NVIDIA Volta GPUs (V100) that perform matrix-multiply-accumulate (MMA) operations on 4x4 submatrices in a single clock cycle. They operate on reduced-precision formats (FP16, BF16, TF32, INT8, INT4, FP8 on Hopper) and accumulate in higher precision (FP32). The primary beneficiaries are GEMM operations (General Matrix Multiply) that dominate deep learning: linear layers, attention Q*K^T products, and convolutions. Tensor Cores require matrix dimensions aligned to multiples of 8 (FP16) or 16 (TF32). Libraries like cuDNN, cuBLAS, and Transformer Engine use Tensor Cores automatically when data is in a compatible format. Mixed-precision training with torch.cuda.amp leverages Tensor Cores transparently for 2-8x throughput gains.

import torch
# Enable TF32 for Tensor Core usage on Ampere+
torch.backends.cuda.matmul.allow_tf32 = True
torch.backends.cudnn.allow_tf32 = True

Mixed-precision training uses FP16 (or BF16 on Ampere+) for forward and backward passes (exploiting Tensor Cores) while maintaining FP32 master weights and optimizer states to preserve numerical stability. PyTorch's Automatic Mixed Precision (AMP) automates this via torch.amp.autocast(device_type='cuda') which casts operations to FP16/BF16 where safe, and torch.cuda.amp.GradScaler which scales loss before backward() to prevent FP16 gradient underflow (values that round to zero), then unscales and clips gradients before optimizer.step(). BF16 (bfloat16) is preferred on Ampere/Hopper as it has the same exponent range as FP32 (avoiding underflow) with less precision, eliminating the need for loss scaling in most cases.

scaler = torch.cuda.amp.GradScaler()
with torch.autocast(device_type='cuda', dtype=torch.float16):
    output = model(batch)
    loss = criterion(output, target)
scaler.scale(loss).backward()
scaler.step(optimizer)
scaler.update()

Gradient checkpointing (also called activation checkpointing) reduces memory usage during backpropagation by not storing all intermediate activations. In standard backprop, all layer activations are kept in memory during the forward pass for use in the backward pass. With checkpointing, activations are discarded during the forward pass and recomputed on-demand during backpropagation. This reduces activation memory from O(n) to O(sqrt(n)) for uniform networks. The trade-off is approximately 33% more compute (recomputing one full forward pass per backward pass). Essential for training very deep or large models that would otherwise OOM. In PyTorch, use torch.utils.checkpoint.checkpoint(fn, *args) around blocks of layers, or use gradient_checkpointing_enable() in HuggingFace Transformers.

import torch.utils.checkpoint as cp

class TransformerBlock(nn.Module):
    def forward(self, x):
        # Recompute attention during backward
        return cp.checkpoint(self._forward, x, use_reentrant=False)
    def _forward(self, x): ...

NCCL (NVIDIA Collective Communications Library) provides highly optimized collective operations (AllReduce, AllGather, Reduce, Broadcast, ReduceScatter) for GPU-to-GPU communication. In data-parallel training, AllReduce aggregates gradients across all GPUs: each GPU contributes its locally computed gradient, and all GPUs receive the averaged gradient. NCCL automatically uses the fastest available transport: NVLink (80-600 GB/s on same-node GPUs), InfiniBand (multi-node), or PCIe (fallback). Ring AllReduce is the standard algorithm: data flows around a ring of GPUs in 2n steps (scatter-reduce + all-gather) where each link carries 2*(n-1)/n of total data volume, maximizing bandwidth utilization. PyTorch's DistributedDataParallel uses NCCL as its default backend.

import torch.distributed as dist
dist.init_process_group('nccl', rank=rank, world_size=world_size)
# AllReduce happens automatically in DDP
model = DDP(model, device_ids=[local_rank])

DataParallel (DP) uses Python threads to replicate the model across GPUs on one machine. The GIL limits parallelism; the parameter server pattern centralizes gradients on GPU 0, creating a bandwidth bottleneck proportional to number of GPUs. DistributedDataParallel (DDP) uses separate processes (one per GPU), avoiding the GIL. Gradients are synchronized via NCCL AllReduce, distributing communication load across all links equally. DDP scales linearly and is faster even on single machines. DDP also supports multi-node training. Use DP only for quick prototyping on a single machine. Always use DDP for production training. Launch DDP with torchrun or torch.distributed.launch, providing RANK, WORLD_SIZE, and LOCAL_RANK environment variables per process.

# DDP launch with torchrun
# torchrun --nproc-per-node=8 train.py
dist.init_process_group('nccl')
local_rank = int(os.environ['LOCAL_RANK'])
model = model.to(local_rank)
model = DDP(model, device_ids=[local_rank])

FSDP (Fully Sharded Data Parallel, PyTorch 1.11+) shards model parameters, gradients, and optimizer states across all GPUs, unlike DDP which replicates the full model on every GPU. FSDP implements ZeRO-3 (Zero Redundancy Optimizer Stage 3) from DeepSpeed. Each GPU holds only 1/N of each parameter. Before a forward or backward pass on a layer, FSDP performs an AllGather to temporarily reconstruct the full parameters, then discards them after the computation. Gradients are immediately reduced via ReduceScatter after each layer's backward pass. FSDP enables training models too large for a single GPU (100B+ parameters) at the cost of higher communication volume vs. DDP. Critical config: sharding_strategy (FULL_SHARD, SHARD_GRAD_OP, NO_SHARD), mixed_precision, and CPU offloading.

from torch.distributed.fsdp import FullyShardedDataParallel as FSDP
from torch.distributed.fsdp import ShardingStrategy
model = FSDP(
    model,
    sharding_strategy=ShardingStrategy.FULL_SHARD,
    device_id=local_rank
)

Standard self-attention materializes the full N×N attention matrix in GPU global memory, requiring O(N²) memory and O(N²) reads/writes — becoming the bottleneck for long sequences. Flash Attention (Dao et al., 2022) implements attention using tiling: it reads blocks of Q, K, V into shared memory, computes partial softmax results (using the online softmax algorithm to maintain correct normalization), and accumulates outputs without ever materializing the full N×N matrix. This reduces global memory reads/writes from O(N²) to O(N) for the attention matrix, making attention IO-bound operations significantly faster. Flash Attention 2 and 3 add further optimizations: parallel over sequence length, reduced warp synchronization, and FP8 support (FA3 on Hopper). Enabled in PyTorch 2.0+ via scaled_dot_product_attention.

# PyTorch 2.0+ - Flash Attention via SDPA
with torch.backends.cuda.sdp_kernel(enable_flash=True):
    output = F.scaled_dot_product_attention(q, k, v, is_causal=True)

Nsight Systems (nsys) provides system-level profiling: timeline visualization of CPU, GPU kernel activity, CUDA API calls, memory transfers, and streams — ideal for identifying pipeline bottlenecks, idle GPU time, and PCIe transfer overhead. Run with nsys profile python train.py and open the .nsys-rep file in the Nsight Systems GUI. Nsight Compute (ncu) provides deep per-kernel analysis: occupancy, memory throughput, cache hit rates, instruction mix, roofline model position, and SASS-level metrics. Run with ncu --set full -k kernel_name python train.py. Use Nsight Systems first to identify the critical kernel, then use Nsight Compute for deep micro-architectural analysis. PyTorch Profiler (torch.profiler) provides Python-level integration with Chrome trace export.

# System-level profiling
nsys profile --trace=cuda,nvtx python train.py

# Kernel-level deep dive
ncu --set full --kernel-name my_kernel python train.py

# PyTorch Profiler
with torch.profiler.profile(activities=[torch.profiler.ProfilerActivity.CUDA]) as p:
    train_step()

GPU shared memory is divided into 32 banks (on modern GPUs) for simultaneous multi-thread access. A bank conflict occurs when two or more threads in the same warp access different addresses that map to the same bank, forcing the accesses to be serialized. The number of conflicts equals the worst-case number of threads mapping to one bank within a warp. Common cause: accessing a 2D shared memory array with a stride equal to the number of banks (e.g., float tile[32][32] accessed as tile[col][row]). Avoidance: add 1-element padding to rows (float tile[32][33]) so row-column swaps break the bank alignment. Accessing the same address by multiple threads (broadcast) does not cause conflicts — it is handled by a hardware broadcast mechanism.

// Avoid bank conflicts with padding
__shared__ float tile[32][33];  // +1 padding breaks bank alignment
// vs conflict-prone:
__shared__ float tile[32][32];  // col-major access causes conflicts

nvidia-smi is NVIDIA's command-line management interface for monitoring GPU state. Key production metrics: gpu_utilization (percent of time at least one warp is active — low means compute starved, not necessarily under-used); memory_utilization (% of time memory is being read/written — high indicates memory bandwidth bottleneck); memory_used/memory_total (VRAM usage — monitor for OOM risk); power_draw vs power_limit (thermal throttling indicator); temperature (exceeding 80-85°C triggers throttling); sm_clock (check for throttling below boost clock); ecc_errors (double-bit errors signal hardware failure). Use nvidia-smi dmon for continuous streaming metrics. Integrate with Prometheus via dcgm-exporter (DCGM) for Grafana dashboards in production ML serving infrastructure.

nvidia-smi --query-gpu=name,utilization.gpu,utilization.memory,memory.used,memory.total,temperature.gpu,power.draw --format=csv -l 1

Move tensors to GPU with tensor.to('cuda') or tensor.cuda(). For specific GPU indices: tensor.to('cuda:1'). Move the entire model with model.to(device). Common mistakes: forgetting to move the model (model on CPU, inputs on GPU → RuntimeError); creating new tensors during the forward pass that default to CPU (use torch.zeros(..., device=x.device)); calling .numpy() on a CUDA tensor (must call .cpu().numpy() first); model output on GPU but loss function on CPU; and mixing devices in DDP (each rank must use its LOCAL_RANK device). Best practice: define device = torch.device('cuda' if torch.cuda.is_available() else 'cpu') at script start and use it consistently. Use model.parameters().__next__().device to check model's current device.

device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
model = model.to(device)
for batch, labels in loader:
    batch = batch.to(device)
    labels = labels.to(device)
    output = model(batch)

The Roofline model visualizes whether a kernel is compute-bound or memory-bandwidth-bound. The x-axis is arithmetic intensity (FLOP/byte of memory accessed); the y-axis is attained FLOP/s. Two ceilings bound performance: the compute roof (peak FLOP/s for the operation type, e.g., FP16 Tensor Core peak) and the memory bandwidth roof (peak bandwidth × arithmetic intensity). If arithmetic intensity places the kernel left of the ridge point (intersection), it is bandwidth-limited; right of the ridge point, it is compute-limited. For memory-bound kernels, optimize data access patterns and reduce memory transactions. For compute-bound kernels, increase parallelism and use Tensor Cores. Use Nsight Compute's Roofline view to plot actual achieved performance against hardware ceilings.

# View in Nsight Compute
ncu --set roofline -k my_kernel python script.py
# Arithmetic intensity = FLOPs / bytes_transferred
# ai = (2 * M * N * K) / ((M*K + K*N + M*N) * 2)  # FP16 matmul

Cooperative Groups (CUDA 9+) provides a structured, flexible API for expressing synchronization between arbitrary subsets of threads beyond the block boundary. Traditional __syncthreads() only synchronizes all threads in a block. Cooperative Groups enables: intra-warp synchronization (cg::coalesced_threads(), cg::tiled_partition<32>(block)), intra-block synchronization (cg::this_thread_block()), and grid-wide synchronization (cg::this_grid(), requires cooperative launch). Grid-wide sync enables split algorithms like multi-pass reductions and iterative solvers without kernel re-launches. Warp-level operations (cg::reduce(), cg::scan()) replace legacy shuffle intrinsics (__shfl_down_sync) with a cleaner API. Essential for implementing efficient reductions, prefix sums, and graph algorithms on GPU.

#include <cooperative_groups.h>
namespace cg = cooperative_groups;
__global__ void reduce(float* data, float* result) {
    auto block = cg::this_thread_block();
    auto warp = cg::tiled_partition<32>(block);
    float val = cg::reduce(warp, data[threadIdx.x], cg::plus<float>());
}

cuDNN (CUDA Deep Neural Network library) is NVIDIA's library of GPU-accelerated primitives for deep learning. It provides highly optimized implementations of: convolution forward and backward, pooling (max, average), activation functions (ReLU, sigmoid, tanh, GELU), normalization (batch, layer, instance norm), dropout, LSTM and GRU cells, attention mechanisms, and softmax. cuDNN uses Tensor Cores automatically when inputs meet alignment requirements. It auto-selects the fastest algorithm for each operation given the input shape via benchmarking (torch.backends.cudnn.benchmark=True runs a one-time search). cuDNN algorithms differ by workspace memory requirements, enabling memory vs. speed trade-offs. PyTorch, TensorFlow, and MXNet all use cuDNN as the GPU backend for neural network operations.

import torch
# Enable cuDNN auto-tuner for fixed input shapes
torch.backends.cudnn.benchmark = True
# Deterministic mode (slower, for reproducibility)
torch.backends.cudnn.deterministic = True

NVIDIA Ampere (A100, A30, RTX 3000 series, 2020) introduced several major improvements. Third-generation Tensor Cores add TF32 (19-bit mantissa of FP32 with FP16 exponent range) for automatic FP32 precision boost, and structured sparsity support (2:4 sparsity pattern) doubling effective Tensor Core FLOP/s for sparse models. The A100 has 80 GB HBM2e memory (vs. V100's 32 GB) at 2 TB/s bandwidth. NVLink 3.0 enables 600 GB/s GPU-to-GPU bandwidth. Multi-Instance GPU (MIG) partitions one A100 into up to 7 independent GPU instances for multi-tenant serving. SXM4 form factor and NVSwitch enable all-to-all NVLink topology for 16-GPU pods. The H100 (Hopper, 2022) further adds FP8 Tensor Cores and the Transformer Engine for adaptive precision.

# Check Tensor Core support
import torch
print(torch.cuda.get_device_capability())  # (8,0) = Ampere A100
# Enable TF32 for Ampere Tensor Core matrix ops
torch.backends.cuda.matmul.allow_tf32 = True

MIG (Multi-Instance GPU) is a hardware partitioning feature introduced on NVIDIA A100 and H100 GPUs that divides a single physical GPU into up to 7 completely isolated GPU instances. Each instance has a dedicated portion of the SM count, memory, L2 cache, bandwidth, and memory controllers — providing hardware-level isolation, not just software-level scheduling. MIG instances cannot interfere with each other's performance or access each other's memory. Production use cases: multi-tenant ML serving platforms where each tenant gets a guaranteed performance slice; running multiple different models simultaneously with predictable latency; cost-efficient inference serving where a full A100 would be over-provisioned for each individual model. Kubernetes supports MIG via the NVIDIA Device Plugin with time-slicing or MIG strategies.

# Configure MIG (requires root)
nvidia-smi mig -cgi 9,9,9,9,9,9,9 -C  # 7x 1g.10gb instances on A100
# Kubernetes: nvidia.com/mig-1g.10gb: 1 in pod resource request

NVLink is NVIDIA's proprietary high-bandwidth GPU interconnect. NVLink 3.0 (A100) provides 600 GB/s bidirectional bandwidth between GPU pairs — 10x more than PCIe 4.0's 64 GB/s. Higher NVLink bandwidth directly benefits AllReduce communication in distributed training: larger gradient tensors can be synchronized faster. NVSwitch is a chip that creates a high-bandwidth non-blocking crossbar connecting all GPUs in a node, enabling any GPU to communicate with any other at full NVLink bandwidth simultaneously (all-to-all communication). An 8-GPU DGX A100 system with NVSwitch provides 2.4 TB/s aggregate NVLink bandwidth. Multi-node training still uses InfiniBand between nodes. The NVLink fabric makes FSDP and tensor parallelism practical within a node.

# Check NVLink topology
nvidia-smi topo -m
# P2P bandwidth test
python -c "import torch; a=torch.randn(1024,device='cuda:0'); b=a.to('cuda:1')"

CUDA atomic operations perform read-modify-write operations on global or shared memory locations guaranteed to be indivisible — no other thread can observe a partial update. Supported operations: atomicAdd, atomicSub, atomicMax, atomicMin, atomicAnd, atomicOr, atomicXor, atomicExch, atomicCAS (compare-and-swap). Necessary when multiple threads update the same memory location (e.g., histogramming where multiple threads increment the same bin, or sparse gradient accumulation). Atomics serialize conflicting accesses and degrade performance when many threads contend for the same location (atomic bottleneck). Mitigation: use shared memory atomics (much faster than global) for intra-block aggregation, then one global atomic per block; or use warp reduction to produce one value per warp before the global atomic.

__global__ void histogram(int* data, int* hist, int N) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < N) {
        atomicAdd(&hist[data[idx]], 1);  // thread-safe increment
    }
}

Custom CUDA kernels can be called from Python without writing C++ using CuPy's RawKernel or PyTorch's torch.utils.cpp_extension. CuPy RawKernel compiles a CUDA C string at runtime using NVRTC and launches it with Python; input/output arrays are CuPy ndarrays. PyTorch's load_inline() compiles a C++ extension with pybind11 bindings, enabling persistent compiled modules. For simpler element-wise or reduction kernels, use torch.compile() (TorchInductor) or Triton, a Python-native GPU kernel DSL. Triton's @triton.jit decorator compiles Python-like kernel code to PTX, with automatic tiling and swizzling optimizations. Triton is now the preferred path for custom high-performance kernels: easier than CUDA C, produces near-cuBLAS-quality code for GEMM-like patterns.

import triton
import triton.language as tl

@triton.jit
def add_kernel(x_ptr, y_ptr, out_ptr, N, BLOCK: tl.constexpr):
    pid = tl.program_id(0)
    offs = pid * BLOCK + tl.arange(0, BLOCK)
    mask = offs < N
    x = tl.load(x_ptr + offs, mask=mask)
    y = tl.load(y_ptr + offs, mask=mask)
    tl.store(out_ptr + offs, x + y, mask=mask)

Paged Attention (vLLM, 2023) solves the KV cache memory fragmentation problem in LLM inference. Traditional attention allocates contiguous GPU memory for each request's KV cache sized to maximum sequence length, wasting 60-80% of memory to worst-case reservation. Paged Attention manages KV cache in fixed-size pages (blocks), similar to OS virtual memory. Each request's KV cache is stored in non-contiguous pages linked by a page table. Pages are allocated on demand, enabling near-zero internal fragmentation and allowing multiple requests to share GPU memory more efficiently. This enables 24x higher throughput than naive implementations by fitting more concurrent requests in GPU memory. All major LLM serving frameworks (vLLM, TGI, TensorRT-LLM) now implement paged or similar blocked KV caching.

from vllm import LLM, SamplingParams
llm = LLM(model='meta-llama/Llama-3-8B', gpu_memory_utilization=0.9)
outputs = llm.generate(['Tell me about paged attention'], SamplingParams(max_tokens=200))

INT8 quantization represents weights and activations as 8-bit integers instead of FP32 or FP16, reducing memory by 4x or 2x respectively and enabling INT8 CUDA Core operations (not Tensor Cores for standard GEMM). On Ampere+, INT8 Tensor Core GEMM provides 2x the FLOP/s of FP16 Tensor Cores. Two quantization modes: post-training quantization (PTQ) calibrates scale factors from a representative dataset without retraining — fast but may degrade accuracy 0.5-2%; quantization-aware training (QAT) simulates quantization during training via fake quantization, recovering accuracy at the cost of retraining. The key challenge is activating outlier values that exceed INT8 range, causing clipping errors — mitigated by LLM.int8() (mixed-precision with FP16 outlier processing) or GPTQ/AWQ for weight-only quantization.

from transformers import BitsAndBytesConfig
quant_config = BitsAndBytesConfig(load_in_8bit=True)
model = AutoModelForCausalLM.from_pretrained(
    'meta-llama/Llama-3-8B',
    quantization_config=quant_config
)

PyTorch uses its own CUDA memory allocator (a caching allocator) that retains freed GPU memory for reuse rather than returning it to the OS immediately. Fragmentation occurs when many small allocations leave gaps between large blocks, preventing large contiguous allocations even if total free memory is sufficient. Mitigation: torch.cuda.empty_cache() returns cached memory to the OS (useful between independent workloads). torch.cuda.memory_stats() and torch.cuda.memory_summary() diagnose fragmentation. Use memory pools sparingly: avoid creating many small intermediate tensors (prefer in-place operations). Enable PYTORCH_CUDA_ALLOC_CONF=expandable_segments:True (PyTorch 2.1+) which allows the allocator to grow and shrink segments dynamically, significantly reducing fragmentation in training with variable sequence lengths.

import torch
import os
# Enable expandable segments (PyTorch 2.1+)
os.environ['PYTORCH_CUDA_ALLOC_CONF'] = 'expandable_segments:True'
# Diagnose fragmentation
print(torch.cuda.memory_summary())
torch.cuda.empty_cache()

Speculative decoding uses a small draft model (7B or smaller) to quickly generate multiple candidate tokens, then verifies them in parallel with the large target model in a single forward pass. If the target model accepts all draft tokens (they match its distribution), multiple tokens are produced in the time of one large model forward pass. Rejected tokens cause a rollback to the divergence point. GPU resource usage: both draft and target models must fit in GPU memory simultaneously. The draft model runs on the same GPU (or a separate smaller GPU in dedicated setups). The key insight is that verification of K tokens by the large model costs roughly the same as generating 1 token (attention is dominated by the K×K matrix for K tokens), making multi-token output nearly free when drafts are accepted.

# vLLM speculative decoding config
llm = LLM(
    model='meta-llama/Llama-3-70B',
    speculative_model='meta-llama/Llama-3-7B',
    num_speculative_tokens=5
)

CUDA events are GPU-side timestamps recorded into a stream. Measure kernel execution time accurately with: cudaEvent_t start, stop; cudaEventRecord(start, stream) before the kernel; cudaEventRecord(stop, stream) after; cudaEventSynchronize(stop) to wait for completion; cudaEventElapsedTime(&ms, start, stop) to retrieve milliseconds. This is more accurate than CPU-side timing (which includes kernel launch overhead, async gaps, and driver latency) because events are timestamped on the GPU itself at stream execution time. In Python: torch.cuda.Event(enable_timing=True) provides the equivalent API. Always warm up the GPU with a few dry runs before timing to avoid CUDA context initialization overhead contaminating measurements.

start = torch.cuda.Event(enable_timing=True)
end = torch.cuda.Event(enable_timing=True)
start.record()
output = model(batch)
end.record()
torch.cuda.synchronize()
print(f'{start.elapsed_time(end):.2f} ms')

torch.compile() (PyTorch 2.0+) compiles a model or function using TorchDynamo (Python bytecode capture), Torch FX graph transformations, and TorchInductor (a code generation backend targeting Triton for GPU). The result is optimized GPU kernel fusion: multiple element-wise operations (add + relu + dropout) are fused into a single GPU kernel, eliminating intermediate memory allocations and global memory round-trips. This is particularly impactful for transformer architectures where attention computation and MLP layers chain many small ops. Expected speedups: 30-100% on training throughput for typical transformer models. Modes: default (balanced), reduce-overhead (reduces Python overhead via graph capture), max-autotune (runs expensive kernel auto-tuning for maximum performance). Requires PyTorch 2.0+ and CUDA, and falls back gracefully for unsupported operations.

model = torch.compile(model, mode='max-autotune')
# Then train normally
for batch in loader:
    loss = model(batch)  # compiled kernels execute

Data loading is CPU-bound when the GPU waits for batches to arrive from the CPU data pipeline. Symptoms: GPU utilization fluctuates (drops to 0% between batches), nvidia-smi shows low sustained GPU utilization despite training being active. Solutions: (1) Use PyTorch DataLoader with num_workers > 0 (multiprocessing data loading, typically 4-8 workers); (2) Use pin_memory=True in DataLoader to allocate batches in pinned memory for faster GPU transfer; (3) prefetch_factor=2 pre-loads 2 batches ahead; (4) NVIDIA DALI (Data Loading Library) performs image decode and augmentation on GPU, eliminating CPU as the bottleneck entirely; (5) pre-tokenize and cache datasets to avoid repeated tokenization; (6) store data in a format optimal for sequential reads (sharded WebDataset, FFCV, or Arrow files).

loader = DataLoader(
    dataset,
    batch_size=256,
    num_workers=8,
    pin_memory=True,
    prefetch_factor=2,
    persistent_workers=True
)

Tensor parallelism (TP) splits individual weight matrices across GPUs within a layer. For a linear layer with weight W, TP column-splits W across N GPUs so each holds W/N columns; outputs are gathered via AllGather. This requires N-GPU synchronization per layer but exposes intra-layer parallelism. Typical TP degree is 4-8 GPUs (within a node to exploit NVLink bandwidth). Pipeline parallelism (PP) assigns different layers to different GPUs (GPU 0 runs layers 1-8, GPU 1 runs layers 9-16). Forward pass data flows sequentially between stages. Micro-batching (GPipe) or interleaved schedules (PipeDream) reduce pipeline bubbles. PP requires less communication than TP (only activations at stage boundaries) but introduces pipeline bubbles. For 100B+ models, 3D parallelism combines data, tensor, and pipeline parallelism simultaneously.

# Megatron-style tensor parallel linear
# Each GPU holds 1/N of the columns
# ColParallelLinear + RowParallelLinear with AllReduce

CUDA graphs record a sequence of GPU operations (kernel launches, memory copies, memsets) into a graph during a capture phase, then replay the entire graph in a single API call. This eliminates the CPU-side kernel launch overhead: each CUDA API call has ~3-5 μs of latency; for models with thousands of small operations per step, this overhead accumulates to milliseconds. CUDA graphs are beneficial when the computation graph is static (same operations, same tensor shapes, same addresses). Enable in PyTorch with torch.cuda.CUDAGraph() or use torch.cuda.make_graphed_callables(). The graph must be recaptured if batch size, sequence length, or other dynamic dimensions change. Particularly impactful for inference serving where low latency is critical and input shapes are fixed.

g = torch.cuda.CUDAGraph()
with torch.cuda.graph(g):
    output = model(static_input)  # capture
# Replay on each inference call
g.replay()

cuBLAS is NVIDIA's CUDA-accelerated BLAS (Basic Linear Algebra Subprograms) library, providing highly optimized implementations of dense matrix operations (GEMM, GEMV, SYMM, TRSM). In deep learning, cuBLAS is called by cuDNN and frameworks for linear layer forward/backward passes, attention Q*K^T products, and embedding lookups. It automatically uses Tensor Cores when the data type is FP16, BF16, or TF32 and matrix dimensions are aligned to multiples of 8 or 16. cuBLAS-Lt (a lightweight version) exposes additional controls for Tensor Core algorithm selection and workspace management. torch.backends.cuda.matmul.allow_tf32=True enables TF32 in cuBLAS for linear operations. cuBLAS uses heuristics and a precomputed algorithm database to select the best algorithm for each matrix shape.

import torch
# Enable TF32 in cuBLAS for Ampere+ GPUs
torch.backends.cuda.matmul.allow_tf32 = True
# cuBLAS handles this automatically:
result = torch.matmul(A_fp16, B_fp16)  # uses Tensor Cores

Debug CUDA OOM with a systematic approach. First, enable memory snapshot: torch.cuda.memory._record_memory_history() captures allocation stack traces; on OOM, call torch.cuda.memory._snapshot() and save it. Visualize with PyTorch's memory viewer (memory_viz.py). Check the error message: it reports the failed allocation size and total/free memory. Common causes: batch size too large (reduce), gradient accumulation missing (add optimizer.zero_grad()); unintended computation graph retention (detach tensors before storing metrics: loss.item() not loss); DataLoader workers holding tensors from previous batches; model inputs not freed between validation steps. Use torch.cuda.max_memory_allocated() to track peak usage. Enable expandable_segments to reduce fragmentation-related OOM.

torch.cuda.memory._record_memory_history()
try:
    train_step()
except torch.cuda.OutOfMemoryError:
    snap = torch.cuda.memory._snapshot()
    torch.save(snap, 'oom_snapshot.pt')

NVIDIA's Transformer Engine (TE) is a library that enables FP8 (8-bit floating point) mixed-precision training on H100 GPUs (Hopper architecture), providing 2x the FLOP/s of FP16 at the same power. FP8 exists in two formats: E4M3 (4-bit exponent, 3-bit mantissa — for forward pass activations and weights) and E5M2 (5-bit exponent, 2-bit mantissa — for backward pass gradients, needing more range). TE's key innovation is automatic per-tensor scale factor management (amax tracking): it maintains a delayed scaling factor per tensor based on historical absolute maxima, converting FP32 master weights to FP8 just-in-time for GEMM operations. TE integrates with PyTorch via transformer_engine.pytorch.Linear and supports TransformerLayer with fused multi-head attention and LayerNorm.

import transformer_engine.pytorch as te
with te.fp8_autocast(enabled=True):
    output = te_linear(input)  # uses FP8 Tensor Cores on H100

Gradient compression reduces the volume of data transmitted during AllReduce, addressing communication bottlenecks in distributed training — especially over slow interconnects (Ethernet vs. InfiniBand). Techniques: (1) FP16 gradient communication: cast FP32 gradients to FP16 before AllReduce, halving bandwidth. (2) 1-bit SGD / PowerSGD: compress gradients to 1-bit representations with error feedback to correct compression bias. (3) TopK sparsification: transmit only the top-K% largest gradients, zeroing the rest; use error accumulation to account for skipped updates. (4) Gradient quantization: quantize to INT8. (5) All-reduce with reduced precision: NCCL supports FP16 all-reduce natively. These techniques trade accuracy and implementation complexity for reduced communication time, most beneficial when the compute-to-communication ratio is low (small batches, slow network).

from torch.distributed.algorithms.ddp_comm_hooks import powerSGD_hook
model = DDP(model)
state = powerSGD_hook.PowerSGDState(process_group=None, matrix_approximation_rank=1)
model.register_comm_hook(state, powerSGD_hook.powerSGD_hook)

A GPU kernel is compute-bound when the bottleneck is the number of arithmetic operations per second (FLOP/s) — more data could be fed but the compute units are saturated. This occurs in large GEMM operations (dense attention, big linear layers) where Tensor Cores are fully utilized. A kernel is memory-bandwidth-bound when the bottleneck is how fast data is read/written from global memory — compute units idle waiting for data. This occurs in element-wise operations (activation functions, layer norm, dropout) and small GEMMs. Diagnosis: Roofline model — kernels left of the ridge point are bandwidth-bound. Optimization: memory-bound kernels benefit from kernel fusion (combine multiple ops to reduce memory round-trips); compute-bound kernels benefit from Tensor Cores, mixed precision, and increased arithmetic intensity via algorithmic changes.

# Check if kernel is memory or compute bound via Nsight Compute
ncu --metrics l1tex__t_bytes_pipe_lsu_mem_global_op_ld.sum,\
smsp__sass_thread_inst_executed_op_fadd_pred_on.sum \
-k my_kernel python script.py

ZeRO (Microsoft DeepSpeed) eliminates memory redundancy in data-parallel training by partitioning optimizer states, gradients, and parameters across GPUs instead of replicating them. Stage 1 (ZeRO-1): partitions optimizer states (Adam's momentum, variance) across GPUs — 4x memory reduction for FP32 Adam. Stage 2 (ZeRO-2): additionally partitions gradients — 8x reduction. Stage 3 (ZeRO-3): partitions all three (optimizer states + gradients + parameters) — linear memory reduction with GPU count (64x with 64 GPUs). PyTorch FSDP implements ZeRO-3 natively. Each stage increases communication overhead: Stage 3 performs AllGather before each forward and ReduceScatter after each backward. ZeRO-Infinity extends Stage 3 with CPU/NVMe offloading, enabling trillion-parameter model training.

# DeepSpeed ZeRO Stage 3
ds_config = {
    'zero_optimization': {
        'stage': 3,
        'offload_optimizer': {'device': 'cpu'},
        'offload_param': {'device': 'cpu'}
    }
}

When a CUDA kernel is launched from the host, the CUDA driver queues the launch command in the associated stream. The GPU's GigaThread Engine reads from the stream and dispatches thread blocks to Streaming Multiprocessors (SMs) that have available slots. Each SM runs its assigned blocks concurrently via warp-level time-multiplexing — warps waiting for memory become inactive while other ready warps execute, hiding memory latency. The pipeline has four main stages: instruction fetch (I-cache), decode, issue/dispatch, and execution. Warps are selected for issue using a zero-overhead warp scheduler — no context switch cost. A kernel completes when all thread blocks finish; the GPU signals completion via a CUDA event or cudaDeviceSynchronize(). Multiple kernels from different streams can overlap on the GPU if both have sufficient occupancy headroom.

// Kernel launch from host
my_kernel<<<numBlocks, threadsPerBlock, sharedMem, stream>>>(args);
// Async: returns immediately; GPU executes in parallel
// Sync when needed:
cudaStreamSynchronize(stream);

PCIe bandwidth (PCIe 4.0: 32 GB/s unidirectional per GPU) is the bottleneck for host-to-device transfers — critical during data loading, model initialization, and multi-node gradient communication over Ethernet (when GPUs are connected via CPU+Ethernet rather than NVLink+InfiniBand). Impact areas: data loading (minimize with pin_memory, async transfers, DALI); gradient all-reduce in multi-node training (minimize by using InfiniBand directly between GPUs bypassing PCIe/CPU with GPUDirect RDMA); CPU offloading (ZeRO-Infinity parameter transfers are PCIe-limited). Within a single node, use NVLink for GPU-GPU communication and keep GPU-to-GPU data on the GPU without CPU round-trips. For inference, pre-load models to GPU and avoid per-request CPU-GPU model weight transfers.

# GPUDirect RDMA for peer-to-peer transfers
# (Requires MLNX OFED + MOFED + compatible NIC)
export NCCL_IB_DISABLE=0
export NCCL_P2P_LEVEL=NVL  # prefer NVLink for P2P

For inference (not training), DDP is overkill. Use tensor parallelism libraries like accelerate (HuggingFace) which distributes model layers across GPUs via device_map='auto', automatically placing each layer on the GPU with most free memory. For batch-parallel inference, use torch.multiprocessing to run independent model replicas on separate GPUs and aggregate results. vLLM and TensorRT-LLM implement optimized multi-GPU inference with tensor parallelism built in. For simple cases: create one model per GPU using model.to(f'cuda:{gpu_id}'), process different batches on each GPU in separate threads, then collect results. Ensure each GPU thread has its own CUDA context (separate processes or careful context management). Flash Attention and PagedAttention further maximize throughput per GPU.

from accelerate import Accelerator
accelerator = Accelerator()
model = AutoModelForCausalLM.from_pretrained(
    'meta-llama/Llama-3-70B',
    device_map='auto',  # auto-shards across GPUs
    torch_dtype=torch.bfloat16
)

Constant memory is a 64 KB read-only cache in CUDA designed for data that all threads read uniformly (the same address at the same time). When all threads in a warp access the same constant memory address, it is served in a single broadcast from the constant cache (extremely fast — similar to L1 cache). If threads access different addresses, accesses are serialized. Declare with __constant__ qualifier; copy from host with cudaMemcpyToSymbol(). Ideal use cases: mathematical constants, lookup tables, model parameters accessed by all threads (e.g., a small weight matrix used by every thread for a dot product). Avoid for large or irregularly accessed data — use shared memory or texture memory instead. Exceeding 64 KB spills to global memory.

__constant__ float weights[256];
// Host: copy once before kernel launch
cudaMemcpyToSymbol(weights, h_weights, 256*sizeof(float));
__global__ void kernel(float* out) {
    out[threadIdx.x] = weights[threadIdx.x % 256]; // broadcast if same idx
}

Texture memory is a read-only global memory region accessed via the texture cache, optimized for spatial locality in 2D/3D access patterns. Key advantages: automatic caching for 2D spatial locality (fetching pixel [x,y] also caches nearby [x±1,y±1] pixels); hardware bilinear and trilinear interpolation at no extra cost; automatic out-of-bounds handling (clamp, wrap, or mirror border modes) without branch code; unnormalized or normalized coordinate addressing. For image processing kernels where each thread accesses a small neighborhood (convolution, bilateral filter), texture cache hit rates are much higher than L1 cache for global memory. In CUDA C++, declare with cudaTextureObject_t; in PyCUDA, use TextureObject. On Ampere+ GPUs, global memory L2 cache improvements have reduced the relative advantage of texture memory for newer architectures.

// Create texture object
cudaTextureObject_t tex;
cudaCreateTextureObject(&tex, &resDesc, &texDesc, nullptr);
// Access in kernel
float val = tex2D<float>(tex, x, y);  // hardware interpolation

NVIDIA GPU memory has multiple levels with dramatically different latencies. Registers: per-thread, zero latency (1 cycle), limited to ~255 per thread (varies by architecture). Shared memory: per-SM, ~20-40 cycles for non-bank-conflict access, 48-164 KB per SM. L1 cache: unified with shared memory on modern GPUs, ~30-50 cycles. L2 cache: across all SMs, ~200 cycles, 40-80 MB on Ampere A100. Global memory (HBM2e/HBM3): ~600-800 cycles, 40-80 GB on A100/H100, bandwidth 2-3.35 TB/s. Constant memory cache: fastest for uniform broadcasts, falls back to global for non-uniform. Texture cache: optimized for 2D spatial access, similar latency to L1. The fundamental optimization insight: move data up the hierarchy (global → shared → registers) to maximize reuse before eviction.

# Memory latency comparison (approximate cycles)
# Registers: 1
# Shared memory: 20-40
# L1 cache: 30-50
# L2 cache: ~200
# Global memory (HBM): 600-800
# Bandwidth: HBM3 = 3.35 TB/s on H100

Warp-level primitives provide fast communication between the 32 threads of a warp without shared memory, leveraging the GPU's register file directly. Key intrinsics: __shfl_down_sync(mask, var, delta) passes a register value from thread T to thread T-delta within the warp; __shfl_xor_sync implements butterfly reductions; __shfl_sync broadcasts from a source lane; __ballot_sync returns a bitmask of which lanes satisfy a condition; __reduce_add_sync, __reduce_max_sync provide hardware-accelerated warp reductions on Hopper. Use for: parallel reduction (sum all 32 values in log2(32)=5 steps), prefix scan, warp-level softmax. Warp primitives are 2-5x faster than shared memory reductions for warp-wide operations because they avoid bank conflicts and __syncthreads() overhead.

// Warp reduction using __shfl_down_sync
__device__ float warp_reduce_sum(float val) {
    for (int offset = 16; offset > 0; offset /= 2)
        val += __shfl_down_sync(0xffffffff, val, offset);
    return val;  // lane 0 holds the sum
}

Kernel fusion combines multiple sequential GPU kernels into a single kernel, eliminating intermediate global memory round-trips. Example: applying LayerNorm + GELU + Dropout sequentially requires three kernel launches and three reads/writes of the full activation tensor from global memory. A fused kernel loads data once, applies all three operations in registers/shared memory, and writes results once — reducing global memory bandwidth by 3x. Modern frameworks automate fusion: PyTorch's torch.compile (TorchInductor) fuses element-wise ops and pointwise compositions automatically via Triton kernels. cuDNN performs conv-BN-ReLU fusion natively. Flash Attention is the extreme example of fusion: computing the entire attention mechanism (Q*K softmax*V) in shared memory without materializing the O(N²) attention matrix in global memory.

# torch.compile automatically fuses these:
@torch.compile
def fused_mlp(x, w1, w2):
    return F.relu(x @ w1.T) @ w2.T  # fused into 2 kernels instead of 3