CUDA and GPU languages

Deep dive

Decision table

Execution and memory model

• Grid → blocks (up to 1024 threads, run on one SM) → warps of 32 lanes. A block's threads share fast shared memory and __syncthreads().
• Memory hierarchy (per H100-class GPU, order of): registers (64K × 32-bit per SM), shared memory/L1 (up to about 228 KB per SM), L2 (tens of MB), HBM (80-192 GB).
• Coalescing: a warp's 32 loads should hit consecutive addresses (one 128-byte transaction); strided access multiplies traffic.
• Divergence: if lanes in a warp take different branches, both paths run serially.
• Occupancy: resident warps per SM limited by registers and shared memory per block; enough to hide latency, not maximal.
• Asynchrony: kernels and copies are queued on streams; overlap compute with transfers using pinned host memory and multiple streams.

Code

CUDA C++: grid-stride loop, warp-shuffle reduction, one atomic per block (nvcc -O3 -arch=sm_80 saxpy.cu).

#include <cstdio>
#include <cstdlib>
#include <vector>
#include <cuda_runtime.h>

#define CUDA_CHECK(x) do { cudaError_t e = (x); if (e != cudaSuccess) { \
  std::fprintf(stderr, "%s:%d %s\n", __FILE__, __LINE__, cudaGetErrorString(e)); std::exit(1); } } while (0)

__device__ float warp_sum(float v) {
  for (int off = 16; off > 0; off >>= 1)
    v += __shfl_down_sync(0xffffffffu, v, off);  // register exchange inside the warp
  return v;
}

// y = a*x + y, and accumulate sum(y) in the same pass (fusion: one read, one write).
__global__ void saxpy_sum(int n, float a, const float* __restrict__ x,
                          float* __restrict__ y, float* total) {
  __shared__ float partial[32];
  float acc = 0.f;
  for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < n; i += gridDim.x * blockDim.x) {
    float v = a * x[i] + y[i];  // neighbouring lanes touch neighbouring addresses: coalesced
    y[i] = v;
    acc += v;
  }
  acc = warp_sum(acc);
  const int lane = threadIdx.x % 32, warp = threadIdx.x / 32;
  if (lane == 0) partial[warp] = acc;
  __syncthreads();
  if (warp == 0) {
    acc = lane < blockDim.x / 32 ? partial[lane] : 0.f;
    acc = warp_sum(acc);
    if (lane == 0) atomicAdd(total, acc);  // float atomics: result order is nondeterministic
  }
}

int main() {
  const int n = 1 << 24;
  std::vector<float> hx(n, 1.f), hy(n, 2.f);
  float *x, *y, *total;
  CUDA_CHECK(cudaMalloc(&x, n * sizeof(float)));
  CUDA_CHECK(cudaMalloc(&y, n * sizeof(float)));
  CUDA_CHECK(cudaMalloc(&total, sizeof(float)));
  CUDA_CHECK(cudaMemcpy(x, hx.data(), n * sizeof(float), cudaMemcpyHostToDevice));
  CUDA_CHECK(cudaMemcpy(y, hy.data(), n * sizeof(float), cudaMemcpyHostToDevice));
  CUDA_CHECK(cudaMemset(total, 0, sizeof(float)));

  int sms = 0;
  CUDA_CHECK(cudaDeviceGetAttribute(&sms, cudaDevAttrMultiProcessorCount, 0));
  saxpy_sum<<<sms * 4, 256>>>(n, 3.f, x, y, total);
  CUDA_CHECK(cudaGetLastError());                       // launch configuration errors
  float h_total = 0.f;
  CUDA_CHECK(cudaMemcpy(&h_total, total, sizeof(float), cudaMemcpyDeviceToHost)); // waits, surfaces faults
  std::printf("sum = %.0f (expected %.0f)\n", h_total, 5.0 * n);
  cudaFree(x); cudaFree(y); cudaFree(total);
}

Triton: the same blocked thinking from Python (pip install triton torch).

import torch
import triton
import triton.language as tl

@triton.jit
def add_relu_kernel(x_ptr, y_ptr, out_ptr, n, BLOCK: tl.constexpr):
    pid = tl.program_id(0)                       # one program instance per block of elements
    offs = pid * BLOCK + tl.arange(0, BLOCK)
    mask = offs < n                              # tail handling
    x = tl.load(x_ptr + offs, mask=mask)
    y = tl.load(y_ptr + offs, mask=mask)
    tl.store(out_ptr + offs, tl.maximum(x + y, 0.0), mask=mask)  # fused: 2 reads, 1 write

def add_relu(x: torch.Tensor, y: torch.Tensor) -> torch.Tensor:
    out = torch.empty_like(x)
    n = out.numel()
    grid = lambda meta: (triton.cdiv(n, meta["BLOCK"]),)
    add_relu_kernel[grid](x, y, out, n, BLOCK=1024)
    return out

x = torch.randn(1 << 24, device="cuda"); y = torch.randn_like(x)
torch.testing.assert_close(add_relu(x, y), torch.relu(x + y))

Efficiency

• Cost: H100 HBM bandwidth is 3.35 TB/s (SXM) to 3.9 TB/s (NVL) (NVIDIA H100 page); PCIe to the host is order of 32-64 GB/s per direction, so keep data on the GPU. Launching a kernel is not free: on a Tesla V100 a 2.9 µs kernel cost 9.6 µs per launch with per-kernel synchronization, 3.8 µs with overlapped launches and 3.4 µs with a captured CUDA Graph (NVIDIA blog). Most kernels are bound by memory bandwidth, not FLOPs.
• Levers:
  • CUDA Graphs — 9.6 → 3.4 µs per short launch on V100 (NVIDIA blog) — loops of many tiny kernels
  • fuse chains of elementwise ops (torch.compile, Triton) — one read and write of HBM instead of one per op (order of) — memory-bound op chains
  • coalesced access and shared-memory tiling — consecutive lanes hit consecutive addresses, fewer transactions (CUDA Best Practices Guide) — any global-memory kernel
  • pinned host memory and streams — copies overlap compute (CUDA Best Practices Guide) — host-device transfer loops
  • FP16/BF16/FP8 on tensor cores — up to 4x faster GPT-3 training on H100 vs the prior generation (NVIDIA claim) — matmul-heavy deep learning
• Naive → standard → state of the art:
  • one tiny kernel per op, pageable copies, FP32 · launch- and bandwidth-bound
  • cuBLAS/cuDNN and torch.compile fusion, pinned memory, BF16 · library-tuned
  • Triton or CUTLASS kernels fused for the model, CUDA Graphs, FP8 where the model tolerates it · near the bandwidth or tensor-core roofline

Build and deploy

• CUDA: nvcc or Clang; compile for the target architectures (-gencode for each sm_XX, plus PTX for forward compatibility). Driver must be at least as new as the toolkit requires.
• Containers: NVIDIA base images (nvidia/cuda:*-runtime) with the NVIDIA Container Toolkit on the host.
• Profile with Nsight Systems (timeline: launches, copies, gaps) then Nsight Compute (per-kernel: bandwidth, occupancy, stalls).
• Triton kernels JIT-compile on first call and cache; use @triton.autotune over block sizes and warps.

Pitfalls

• Forgetting that launches are asynchronous: timing without cudaDeviceSynchronize()/events measures nothing; errors appear at the next sync.
• Host-device copies inside loops; pageable host memory is slower and can't overlap with compute (use pinned memory).
• Uncoalesced or strided global memory access; array-of-structs layouts (use struct-of-arrays).
• Shared-memory bank conflicts (32 banks of 4 bytes); pad tiles by one element.
• Too many tiny kernels (launch-bound); fuse or capture with CUDA Graphs.
• Float atomics and parallel reductions are nondeterministic; tests need tolerances.
• Assuming FP32 tensor-core matmuls are exact FP32 (TF32 is often enabled by default in frameworks).
• Writing a custom kernel before checking whether the time is actually in that op (profile first).

Implementations

Sources

• NVIDIA, CUDA C++ Programming Guide docs and CUDA C++ Best Practices Guide docs.
• Tillet, Kung, Cox, Triton: An Intermediate Language and Compiler for Tiled Neural Network Computations, MAPL 2019. paper · code
• Williams, Waterman, Patterson, Roofline: An Insightful Visual Performance Model for Multicore Architectures, CACM 2009. paper
• Kirk, Hwu, Programming Massively Parallel Processors, Morgan Kaufmann (4th ed. 2022, with El Hajj).
• AMD, HIP Programming Guide (ROCm documentation). docs
• NVIDIA, NVIDIA H100 Tensor Core GPU. site · Gray, Getting Started with CUDA Graphs, NVIDIA Developer Blog. post