A collection of CUDA kernels implementing common deep-learning primitives, with benchmarking, correctness checks, and PyTorch bindings.
Three implementations of element-wise float addition on a 1D array:
- Naive — one thread per element, simple index mapping
- Grid-stride — threads loop over the array; saturates the GPU with fewer blocks
- Float4 — vectorized 128-bit loads/stores for higher memory bandwidth
- Naive — each thread computes one output element with uncached global memory reads
- Tiled-16 — 16×16 shared-memory tiling; reduces DRAM traffic by 16×
- Tiled-32 — 32×32 variant for higher occupancy on larger matrices
Tiling reduces loads from
2Nto2N/TILEper output element.
Numerically stable, single-pass row softmax using warp shuffle intrinsics (__shfl_down_sync):
- Parallel reduction → row max (for numerical stability)
exp(x - max)written to output in registers- Parallel reduction → sum of exponentials
- In-place divide by sum
Warp shuffles replace shared memory for intra-warp reductions (~4 cycles vs ~20 for shared memory round-trips).
Fused layer normalization using Welford's online algorithm for numerically stable, single-pass mean and variance:
y = gamma * (x - mean) / sqrt(var + eps) + beta
Warp-level reduction with shared memory fallback for cross-warp accumulation.
Tanh-approximation GELU:
GELU(x) = 0.5 * x * (1 + tanh(sqrt(2/π) * (x + 0.044715 * x³)))
- Naive — scalar, one element per thread
- Vec4 — float4 vectorized loads for 4× memory throughput
.
├── kernels/
│ ├── vector_add.cu
│ ├── matmul_naive.cu
│ ├── matmul_tiled.cu
│ ├── softmax.cu
│ ├── layernorm.cu
│ └── gelu.cu
├── benchmarks/
│ └── benchmark.cu # Unified benchmark harness
├── bindings/
│ └── torch.cpp # PyTorch C++ extension bindings
├── include/
│ └── utils.h # CUDA_CHECK, CUDATimer, fill_random, arrays_close
└── analysis/
├── roofline.md
├── bank_conflicts.md
└── memory_coalescing.md
| Utility | Description |
|---|---|
CUDA_CHECK(call) |
Wraps any CUDA API call; prints file/line and exits on error |
CUDATimer |
RAII wrapper around cudaEvent_t; call .begin() / .end() |
fill_random(ptr, n, lo, hi) |
Fills a host float array with uniform random values |
arrays_close(a, b, n, tol) |
Element-wise correctness check with configurable tolerance |
optimal_block_size(kernel) |
Wraps cudaOccupancyMaxPotentialBlockSize |
Each kernel is exposed as a Python function via pybind11:
import cuda_kernel_suite as ck
y = ck.Softmax(x) # [B, N] → [B, N]
y = ck.LayerNorm(x, gamma, beta, eps=1e-5) # [B, N] → [B, N]
y = ck.GeLU(x) # [N] → [N]
c = ck.Matmul_Tiled(a, b) # [N, N] × [N, N] → [N, N]All inputs must be contiguous float32 CUDA tensors.
# Individual kernels
nvcc -O3 -arch=sm_80 kernels/softmax.cu -o softmax
nvcc -O3 -arch=sm_80 kernels/layernorm.cu -o layernorm
nvcc -O3 -arch=sm_80 kernels/gelu.cu -o gelu
nvcc -O3 -arch=sm_80 kernels/matmul_tiled.cu -o matmul_tiled
# PyTorch extension
python setup.py build_ext --inplaceEach kernel prints a one-liner with timing and throughput:
[Softmax] | B = 512 | N = 1024 | Time = 0.041 ms | BW = 102.3 GB/s
[GELU Naive] | N = 4194304 | Time = 0.312 ms | BW = 214.7 GB/s
[Tiled-16] | N = 1024 | Time = 1.23 ms | GFLOP/s = 1745.2
[Layer Norm] | B = 512 | N = 768 | Time = 0.038 ms | BW = 210.1 GB/s
Correctness is verified either by row-sum checks (softmax), mean/variance checks (layernorm), or direct comparison against a CPU reference (arrays_close).