Compute bound: vLLM

.gitignore

@@ -143,6 +143,8 @@ ENV/
 env.bak/
 venv.bak/
 
+*.nsys-rep
+
 # Spyder project settings
 .spyderproject
 .spyproject
@@ -185,4 +187,4 @@ outputs
 external
 
 # results and analysis
-results
+*.deb

README.md

@@ -16,3 +16,4 @@ A repository to hold my examples and projects in the blog [agramunt.me](https://
 - [CUDA Utils](https://agramunt.me/posts/cuda-utils/) - [Code](https://github.com/SebastiaAgramunt/blogging-code/tree/main/cuda-utils)
 - [CUDA Performance](https://agramunt.me/posts/cuda-performance/) - [Code](https://github.com/SebastiaAgramunt/blogging-code/tree/main/cuda-performance)
 - [BLAS and LAPACKE](https://agramunt.me/posts/blas-lapack/) - [Code](https://github.com/SebastiaAgramunt/blogging-code/tree/main/blas-lapacke)
+- [CUDA GMMA and Roofline model](https://agramunt.me/posts/posts/cuda-mma/) - [Code](https://github.com/SebastiaAgramunt/blogging-code/tree/main/cuda-mma)

cuda-mma/.gitignore

@@ -0,0 +1,3 @@
+build/
+output/
+vendor/

cuda-mma/Makefile

@@ -0,0 +1,60 @@
+# ---------------------------------------------------------------------------
+# Build system for cuda-mma
+# ---------------------------------------------------------------------------
+
+NVCC      := nvcc
+CXX_FLAGS := -std=c++17 -O3
+
+# sm_80 = Ampere (A100 / RTX 30xx)   sm_86 = RTX 30xx consumer
+# sm_70 = Volta (V100)                sm_75 = Turing (RTX 20xx / T4)
+# sm_89 = Ada (RTX 40xx)
+# Override on the command line: make ARCH=sm_75
+ARCH      ?= sm_80
+
+CUTLASS_DIR := vendor/cutlass
+
+NVCC_FLAGS := $(CXX_FLAGS) -arch=$(ARCH) \
+              -rdc=true \
+              --generate-line-info \
+              -Xcompiler -Wall \
+              -diag-suppress 20013,20015 \
+              -I. -Iinclude -I$(CUTLASS_DIR)/include
+
+BIN_DIR   := build/bin
+TARGET    := $(BIN_DIR)/cuda_mma
+TEST      := $(BIN_DIR)/test_correctness
+SRCS      := main.cu src/kernels.cu src/utils.cu src/benchmarks.cu src/benchmark_cutlass.cu
+TEST_SRCS := test/test_correctness.cu src/kernels.cu src/utils.cu
+HDRS      := include/cuda_check.cuh include/kernels.cuh \
+             include/timer.h include/utils.h include/benchmarks.h
+
+.PHONY: all clean run test profile
+
+all: $(TARGET) $(TEST)
+
+$(BIN_DIR):
+	mkdir -p $(BIN_DIR)
+
+$(TARGET): $(SRCS) $(HDRS) | $(BIN_DIR)
+	$(NVCC) $(NVCC_FLAGS) -o $@ $(SRCS) -lcublas -lopenblas -Xlinker -rpath,/usr/local/cuda/targets/x86_64-linux/lib
+
+$(TEST): $(TEST_SRCS) $(HDRS) | $(BIN_DIR)
+	$(NVCC) $(NVCC_FLAGS) -o $@ $(TEST_SRCS)
+
+run: $(TARGET)
+	./$(TARGET)
+
+test: $(TEST)
+	./$(TEST)
+
+# Requires Nsight Compute (ncu).
+profile: $(TARGET)
+	ncu --metrics \
+	    l1tex__t_bytes_pipe_lsu_mem_global_op_ld.sum,\
+	    l1tex__t_bytes_pipe_lsu_mem_global_op_st.sum,\
+	    sm__warps_active.avg.pct_of_peak_sustained_active,\
+	    smsp__sass_thread_inst_executed_op_ffma_pred_on.sum \
+	    ./$(TARGET)
+
+clean:
+	rm -rf build

cuda-mma/include/benchmarks.h

@@ -0,0 +1,12 @@
+#pragma once
+
+float benchmark_naive(int S, const float* h_A, const float* h_B);
+float benchmark_tiled(int S, const float* h_A, const float* h_B);
+float benchmark_coalesced(int S, const float* h_A, const float* h_B);
+float benchmark_cublas(int S, const float* h_A, const float* h_B);
+float benchmark_cblas(int S, const float* h_A, const float* h_B);
+float benchmark_coarsened(int S, const float* h_A, const float* h_B);
+float benchmark_cutlass_fp16(int S, const float* h_A, const float* h_B);
+float benchmark_cutlass_fp32(int S, const float* h_A, const float* h_B);
+float benchmark_cutlass_tf32(int S, const float* h_A, const float* h_B);
+

cuda-mma/include/cuda_check.cuh

@@ -0,0 +1,40 @@
+// cuda_check.cuh
+#pragma once
+
+#include <cuda_runtime.h>
+#include <cstdio>
+#include <cstdlib>
+
+// ── Host-side API checks ──────────────────────────────────────────────
+#define CUDA_CHECK(call)                                               \
+  do {                                                                 \
+    cudaError_t _e = (call);                                           \
+    if (_e != cudaSuccess) {                                           \
+      fprintf(stderr, "[CUDA ERROR] %s:%d  %s\n  → %s\n",              \
+              __FILE__, __LINE__, #call,                               \
+              cudaGetErrorString(_e));                                 \
+      std::exit(EXIT_FAILURE);                                         \
+    }                                                                  \
+  } while (0)
+
+// ── Kernel launch checks ──────────────────────────────────────────────
+#define CHECK_LAST_ERROR()                                             \
+  do {                                                                 \
+    cudaError_t _e = cudaGetLastError();                               \
+    if (_e != cudaSuccess) {                                           \
+      fprintf(stderr, "[KERNEL LAUNCH ERROR] %s:%d  → %s\n",           \
+              __FILE__, __LINE__, cudaGetErrorString(_e));             \
+      std::exit(EXIT_FAILURE);                                         \
+    }                                                                  \
+  } while (0)
+
+// ── Post-kernel execution checks ──────────────────────────────────────
+#define CHECK_SYNC()                                                   \
+  do {                                                                 \
+    cudaError_t _e = cudaDeviceSynchronize();                          \
+    if (_e != cudaSuccess) {                                           \
+      fprintf(stderr, "[KERNEL EXEC ERROR] %s:%d  → %s\n",             \
+              __FILE__, __LINE__, cudaGetErrorString(_e));             \
+      std::exit(EXIT_FAILURE);                                         \
+    }                                                                  \
+  } while (0)

cuda-mma/include/kernels.cuh

@@ -0,0 +1,52 @@
+#pragma once
+
+#include <cuda_runtime.h>
+
+
+__global__ void sgemm_naive(
+    size_t M,             // Number of rows in A and C
+    size_t N,             // Number of columns in B and C
+    size_t K,             // Number of columns in A and rows in B
+    float alpha,       // Scaling factor for the product of A and B
+    const float *A,    // [M x K] row-major
+    const float *B,    // [K x N] row-major
+    float beta,        // Scaling factor for C
+    float *C);         // [M x N] row-major  (in-out: C = alpha*A*B + beta*C)
+
+
+template <int BLOCKSIZE>
+__global__ void sgemm_coalesced(
+    size_t M,             // Number of rows in A and C
+    size_t N,             // Number of columns in B and C
+    size_t K,             // Number of columns in A and rows in B
+    float alpha,       // Scaling factor for the product of A and B
+    const float *A,    // [M x K] row-major
+    const float *B,    // [K x N] row-major
+    float beta,        // Scaling factor for C
+    float *C);         // [M x N] row-major  (in-out: C = alpha*A*B + beta*C)
+
+
+template <int TILE_SIZE>
+__global__ void sgemm_tiled(
+    size_t M,             // Number of rows in A and C
+    size_t N,             // Number of columns in B and C
+    size_t K,             // Number of columns in A and rows in B
+    float alpha,       // Scaling factor for the product of A and B
+    const float *A,    // [M x K] row-major
+    const float *B,    // [K x N] row-major
+    float beta,        // Scaling factor for C
+    float *C);         // [M x N] row-major  (in-out: C = alpha*A*B + beta*C)
+
+
+// Thread-coarsened tiling: each thread computes TM consecutive rows of one C column.
+// Block tile BM×BN, K-tile BK. Threads per block: BM*BN/TM.
+template <int BM, int BN, int BK, int TM>
+__global__ void sgemm_coarsened(
+    size_t M,
+    size_t N,
+    size_t K,
+    float alpha,
+    const float *A,    // [M x K] row-major
+    const float *B,    // [K x N] row-major
+    float beta,
+    float *C);         // [M x N] row-major  (in-out: C = alpha*A*B + beta*C)

cuda-mma/include/timer.h

@@ -0,0 +1,29 @@
+#pragma once
+
+#include <cuda_runtime.h>
+
+// ---------------------------------------------------------------------------
+// CUDA event-based timer.  Usage:
+//
+//   GpuTimer t;
+//   t.start();
+//   kernel<<<...>>>(...);
+//   float ms = t.stop();   // blocks until kernel finishes
+// ---------------------------------------------------------------------------
+struct GpuTimer {
+    cudaEvent_t _start, _stop;
+
+    GpuTimer()  { cudaEventCreate(&_start); cudaEventCreate(&_stop); }
+    ~GpuTimer() { cudaEventDestroy(_start); cudaEventDestroy(_stop);  }
+
+    void start() { cudaEventRecord(_start); }
+
+    // Returns elapsed milliseconds.
+    float stop() {
+        cudaEventRecord(_stop);
+        cudaEventSynchronize(_stop);
+        float ms = 0.0f;
+        cudaEventElapsedTime(&ms, _start, _stop);
+        return ms;
+    }
+};

cuda-mma/include/utils.h

@@ -0,0 +1,17 @@
+#pragma once
+
+#include <cstdlib>
+
+// Utility macros and functions for the CUDA MMA example.
+#define CEIL_DIV(M, N) (((M) + (N)-1) / (N))
+
+// Fill an array with uniform random floats in [-1, 1].
+void fill_random(float* data, int n);
+
+// CPU reference: C += A * B  (row-major, m×k × k×n → m×n).
+void matmul_cpu(const float* A, const float* B, float* C, int m, int n, int k);
+
+// Element-wise comparison with absolute + relative tolerance.
+// Prints the first few mismatches to stderr; returns true if all match.
+bool verify(const float* ref, const float* gpu, int total_elements,
+            float atol = 1e-5f, float rtol = 1e-5f);

cuda-mma/main.cu

@@ -0,0 +1,85 @@
+#include <cstdio>
+#include <cstdlib>
+#include <filesystem>
+
+#include "utils.h"
+#include "benchmarks.h"
+
+#define WARMUP 3
+#define ITERS  5
+
+using BenchmarkFn = float (*)(int, const float*, const float*);
+
+static void run_sweep(const char* label, BenchmarkFn fn,
+                      const int* sizes, int n_sizes,
+                      const float* h_A, const float* h_B,
+                      const char* csv_path,
+                      size_t elem_bytes_ab = sizeof(float),
+                      size_t elem_bytes_cd = sizeof(float))
+{
+    printf("%s roofline sweep (%d iters, %d warmup)\n", label, ITERS, WARMUP);
+    printf("%-6s  %9s  %9s  %11s  %9s\n",
+           "Size", "Time(ms)", "GFLOP/s", "BW(GB/s)", "AI(F/B)");
+    printf("------  ---------  ---------  -----------  ---------\n");
+
+    FILE* csv = fopen(csv_path, "w");
+    fprintf(csv, "size,time_ms,gflops,bandwidth_gbs,arithmetic_intensity\n");
+
+    for (int i = 0; i < n_sizes; ++i) {
+        int S = sizes[i];
+
+        float avg_ms = fn(S, h_A, h_B);
+
+        double flops = 2.0 * S * S * S;
+        double bytes = 2.0 * S * S * elem_bytes_ab    // A + B reads
+                      + 2.0 * S * S * elem_bytes_cd;   // C read + D write
+        double ai        = flops / bytes;
+        double gflops    = flops / (avg_ms * 1e-3) / 1e9;
+        double bandwidth = bytes / (avg_ms * 1e-3) / 1e9;
+
+        printf("%-6d  %9.3f  %9.1f  %11.1f  %9.2f\n",
+               S, avg_ms, gflops, bandwidth, ai);
+        fprintf(csv, "%d,%.3f,%.3f,%.3f,%.4f\n",
+                S, avg_ms, gflops, bandwidth, ai);
+    }
+    printf("\n");
+    fclose(csv);
+}
+
+int main() {
+    srand(42);
+
+    std::filesystem::create_directories("output");
+
+    const int sizes[] = {128, 256, 512, 1024, 2048, 4096, 8192, 16384, 32768};
+    const int N_SIZES = sizeof(sizes) / sizeof(sizes[0]);
+
+    const int S_max = sizes[N_SIZES - 1];
+    float* h_A = new float[(size_t)S_max * S_max];
+    float* h_B = new float[(size_t)S_max * S_max];
+    fill_random(h_A, S_max * S_max);
+    fill_random(h_B, S_max * S_max);
+
+    run_sweep("Naive",     benchmark_naive,     sizes, N_SIZES, h_A, h_B, "output/naive.csv");
+    run_sweep("Tiled",     benchmark_tiled,     sizes, N_SIZES, h_A, h_B, "output/tiled.csv");
+    run_sweep("Coalesced", benchmark_coalesced, sizes, N_SIZES, h_A, h_B, "output/coalesced.csv");
+    run_sweep("cuBLAS",    benchmark_cublas,    sizes, N_SIZES, h_A, h_B, "output/cublas.csv");
+    run_sweep("Coarsened", benchmark_coarsened, sizes, N_SIZES, h_A, h_B, "output/coarsened.csv");
+    run_sweep("CUTLASS_fp32",   benchmark_cutlass_fp32,   sizes, N_SIZES, h_A, h_B, "output/cutlass_fp32.csv");
+
+    // CUTLASS version using tf32 tensor cores; A/B are bitcast (still 4 bytes), not narrowed.
+    run_sweep("CUTLASS_tf32",   benchmark_cutlass_tf32,   sizes, N_SIZES, h_A, h_B, "output/cutlass_tf32.csv");
+
+    // CUTLASS version for half_t (2-byte) operands for A/B; C/D stay float (4 bytes).
+    run_sweep("CUTLASS_fp16",   benchmark_cutlass_fp16,   sizes, N_SIZES, h_A, h_B, "output/cutlass_fp16.csv", 2, sizeof(float));
+
+    // CBLAS runs on CPU; cap at 4096 to keep runtime reasonable
+    const int cblas_sizes[] = {128, 256, 512, 1024, 2048, 4096, 8192};
+    const int N_CBLAS = sizeof(cblas_sizes) / sizeof(cblas_sizes[0]);
+    run_sweep("CBLAS",     benchmark_cblas,     cblas_sizes, N_CBLAS, h_A, h_B, "output/cblas.csv");
+
+    delete[] h_A;
+    delete[] h_B;
+
+    return 0;
+}

cuda-mma/scripts/download_vendor.sh

@@ -0,0 +1,28 @@
+#!/usr/bin/env bash
+# Downloads vendored dependencies into vendor/.
+# Re-running is safe: existing checkouts are skipped.
+set -euo pipefail
+
+REPO_ROOT="$(cd "$(dirname "$0")/.." && pwd)"
+VENDOR_DIR="$REPO_ROOT/vendor"
+
+CUTLASS_TAG="v4.4.2"
+CUTLASS_DIR="$VENDOR_DIR/cutlass"
+
+echo "==> Downloading CUTLASS $CUTLASS_TAG into $CUTLASS_DIR ..."
+
+if [ -d "$CUTLASS_DIR/.git" ]; then
+    echo "    Already present, skipping."
+else
+    mkdir -p "$VENDOR_DIR"
+    git clone \
+        --filter=blob:none \
+        --sparse \
+        --depth=1 \
+        --branch "$CUTLASS_TAG" \
+        https://github.com/NVIDIA/cutlass.git \
+        "$CUTLASS_DIR"
+    git -C "$CUTLASS_DIR" sparse-checkout set include
+fi
+
+echo "Done. vendor/ is ready."

cuda-mma/scripts/install_env.sh

@@ -0,0 +1,8 @@
+#!/bin/bash
+
+THIS_DIR=$(dirname "$(realpath "$0")")
+
+rm -rf "$THIS_DIR"/.venv
+uv venv "$THIS_DIR"/.venv --python 3.12
+
+cd "$THIS_DIR" && uv pip install matplotlib path pandas