From 09cd3b4a2dbc37c07ab2eb3c892855a8cc4e126a Mon Sep 17 00:00:00 2001 From: Sebas <12949040+SebastiaAgramunt@users.noreply.github.com> Date: Tue, 31 Mar 2026 20:26:54 -0700 Subject: [PATCH 01/21] compute bound --- cuda-memory-compute-bound/include/kernels.h | 43 +++ cuda-memory-compute-bound/include/utils.h | 29 ++ cuda-memory-compute-bound/scripts/analyze.py | 326 +++++++++++++++++ cuda-memory-compute-bound/scripts/compile.sh | 41 +++ cuda-memory-compute-bound/scripts/execute.sh | 50 +++ cuda-memory-compute-bound/src/kernels.cu | 102 ++++++ cuda-memory-compute-bound/src/main.cu | 356 +++++++++++++++++++ 7 files changed, 947 insertions(+) create mode 100644 cuda-memory-compute-bound/include/kernels.h create mode 100644 cuda-memory-compute-bound/include/utils.h create mode 100755 cuda-memory-compute-bound/scripts/analyze.py create mode 100755 cuda-memory-compute-bound/scripts/compile.sh create mode 100755 cuda-memory-compute-bound/scripts/execute.sh create mode 100644 cuda-memory-compute-bound/src/kernels.cu create mode 100644 cuda-memory-compute-bound/src/main.cu diff --git a/cuda-memory-compute-bound/include/kernels.h b/cuda-memory-compute-bound/include/kernels.h new file mode 100644 index 0000000..a1b7f71 --- /dev/null +++ b/cuda-memory-compute-bound/include/kernels.h @@ -0,0 +1,43 @@ +#ifndef KERNELS_H +#define KERNELS_H + +#include "utils.h" +#include + +// ─── Benchmark kernel ──────────────────────────────────────────────────────── +// +// Each thread: +// 1. Loads a[i] and b[i] → 2 × sizeof(float) = 8 bytes read +// 2. Computes val = a[i] + b[i] (1 FLOP) +// 3. Runs m FMA iterations: val = val * alpha + beta (2 FLOPs each) +// 4. Stores c[i] → 1 × sizeof(float) = 4 bytes written +// +// Totals per element +// Bytes: 3 × 4 = 12 +// FLOPs: 1 + 2·m +// Arithmetic intensity: (1 + 2·m) / 12 [FLOPs / byte] +// +// The GPU hides FMA latency (~4 cycles) through inter-warp parallelism, so for +// large N the kernel will approach peak FP32 throughput once AI exceeds the +// ridge point. +void vector_add_fma( + const float* d_a, const float* d_b, float* d_c, + int N, int m, int threads_per_block, + cudaEvent_t start, cudaEvent_t stop, float& elapsed_ms); + +// ─── Peak-bandwidth kernel ──────────────────────────────────────────────────── +// +// Pure memory copy: c[i] = a[i] +// Bytes: 2 × 4 = 8 per element +// FLOPs: 0 → arithmetic intensity ≈ 0 +// +// Used to measure empirical peak HBM / GDDR bandwidth. +void bandwidth_test( + const float* d_a, float* d_b, + int N, int threads_per_block, + cudaEvent_t start, cudaEvent_t stop, float& elapsed_ms); + +// GPU warm-up (first kernel launch incurs driver-init overhead) +void warmup(); + +#endif diff --git a/cuda-memory-compute-bound/include/utils.h b/cuda-memory-compute-bound/include/utils.h new file mode 100644 index 0000000..3e70ca5 --- /dev/null +++ b/cuda-memory-compute-bound/include/utils.h @@ -0,0 +1,29 @@ +#ifndef UTILS_H +#define UTILS_H + +#include +#include + +inline void check(cudaError_t err, const char* const func, const char* const file, const int line) +{ + if (err != cudaSuccess) { + std::cerr << "CUDA Runtime Error at: " << file << ":" << line << "\n"; + std::cerr << cudaGetErrorString(err) << " " << func << "\n"; + std::exit(EXIT_FAILURE); + } +} + +inline void checkLast(const char* const file, const int line) +{ + cudaError_t const err{cudaGetLastError()}; + if (err != cudaSuccess) { + std::cerr << "CUDA Runtime Error at: " << file << ":" << line << "\n"; + std::cerr << cudaGetErrorString(err) << "\n"; + std::exit(EXIT_FAILURE); + } +} + +#define CHECK_CUDA_ERROR(expr) (check((expr), #expr, __FILE__, __LINE__)) +#define CHECK_LAST_CUDA_ERROR() (checkLast(__FILE__, __LINE__)) + +#endif diff --git a/cuda-memory-compute-bound/scripts/analyze.py b/cuda-memory-compute-bound/scripts/analyze.py new file mode 100755 index 0000000..f0a01ee --- /dev/null +++ b/cuda-memory-compute-bound/scripts/analyze.py @@ -0,0 +1,326 @@ +#!/usr/bin/env python3 +""" +analyze.py — Visualise memory-bound vs compute-bound GPU behaviour. + +Reads: + results/roofline.csv + results/bw_saturation.csv + results/device_info.csv + +Produces: + results/roofline.png — classic roofline model + results/bw_vs_m.png — bandwidth & GFLOP/s vs m + results/bw_saturation.png — bandwidth vs vector size (N) + results/time_breakdown.png — compute time vs m with regime labels +""" + +import pathlib +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +import matplotlib.ticker as ticker + +THIS_DIR = pathlib.Path(__file__).parent.resolve() +ROOT_DIR = THIS_DIR.parent.resolve() +RESULTS = ROOT_DIR / "results" + +plt.rcParams.update({ + "font.size": 13, + "axes.titlesize": 15, + "axes.labelsize": 13, + "legend.fontsize": 11, + "figure.dpi": 120, +}) + +# ───────────────────────────────────────────────────────────────────────────── +# Load data +# ───────────────────────────────────────────────────────────────────────────── + +def load_device_info(path: pathlib.Path) -> dict: + df = pd.read_csv(path, index_col="key") + return {k: float(v) if v.replace(".", "", 1).isdigit() else v + for k, v in df["value"].items()} + + +def load_csv(name: str) -> pd.DataFrame: + path = RESULTS / name + if not path.exists(): + raise FileNotFoundError(f"Missing {path} — run execute.sh first") + return pd.read_csv(path) + + +# ───────────────────────────────────────────────────────────────────────────── +# Plot 1 — Roofline model +# ───────────────────────────────────────────────────────────────────────────── + +def plot_roofline(df: pd.DataFrame, info: dict): + """ + Classic roofline: GFLOP/s on y-axis, arithmetic intensity on x-axis. + + Two hardware ceilings: + • Memory-bandwidth ceiling: GFLOP/s = BW_GB/s × AI + • Compute ceiling: GFLOP/s = peak_FP32 + + Points below the ridge are memory-bound; above are compute-bound. + """ + peak_bw = info["theoretical_bw_gbs"] # GB/s + peak_flops = info["theoretical_fp32_gflops"] # GFLOP/s + emp_bw = info.get("empirical_bw_gbs", peak_bw) + ridge_pt = info["ridge_point_flop_per_byte"] + + ai_range = np.logspace(-2, 4, 500) + + # Roofline ceiling + roof_mem = peak_bw * ai_range # memory-bandwidth ceiling (GFLOP/s) + roof_compute = np.full_like(ai_range, peak_flops) + roofline = np.minimum(roof_mem, roof_compute) + + fig, ax = plt.subplots(figsize=(10, 6), layout="constrained") + + ax.loglog(ai_range, roofline, "k-", lw=2.5, label="Roofline (theoretical)", zorder=3) + ax.loglog(ai_range, emp_bw * ai_range, "k--", lw=1.5, + label=f"Empirical BW ({emp_bw:.0f} GB/s)", zorder=3) + ax.axhline(peak_flops, color="steelblue", lw=1.5, ls=":", + label=f"Peak FP32 ({peak_flops:.0f} GFLOP/s)") + + # Mark ridge point + ax.axvline(ridge_pt, color="gray", lw=1.2, ls="--") + ax.text(ridge_pt * 1.08, peak_flops * 0.55, f"Ridge\n{ridge_pt:.1f} FLOP/byte", + fontsize=10, color="gray", va="center") + + # Shade regimes + ax.axvspan(ai_range[0], ridge_pt, alpha=0.06, color="royalblue", + label="Memory-bound region") + ax.axvspan(ridge_pt, ai_range[-1], alpha=0.06, color="tomato", + label="Compute-bound region") + + # Measured points, coloured by m + m_vals = df["m"].values + sc = ax.scatter(df["arithmetic_intensity"], df["achieved_gflops"], + c=np.log2(m_vals.astype(float)), cmap="plasma", + s=70, zorder=5, edgecolors="k", linewidths=0.4) + + # Annotate a few interesting m values + for m_label in [1, 32, 128, 1024, 8192]: + row = df[df["m"] == m_label] + if row.empty: + continue + ax.annotate(f"m={m_label}", + xy=(row["arithmetic_intensity"].values[0], + row["achieved_gflops"].values[0]), + xytext=(6, 4), textcoords="offset points", + fontsize=8, color="black") + + cbar = fig.colorbar(sc, ax=ax, pad=0.02) + cbar.set_label("log₂(m) — FMA iterations per element") + + device = info.get("device_name", "GPU") + ax.set_xlabel("Arithmetic Intensity [FLOP / byte]") + ax.set_ylabel("Achieved Throughput [GFLOP/s]") + ax.set_title(f"Roofline Model — {device}\nvector_add + m FMAs per element") + ax.legend(loc="lower right") + ax.grid(True, which="both", ls="--", lw=0.6, alpha=0.5) + ax.set_xlim(0.05, 5000) + ax.set_ylim(1, peak_flops * 3) + + fig.savefig(RESULTS / "roofline.png", bbox_inches="tight") + plt.close(fig) + print("Saved: roofline.png") + + +# ───────────────────────────────────────────────────────────────────────────── +# Plot 2 — Bandwidth and GFLOP/s vs m +# ───────────────────────────────────────────────────────────────────────────── + +def plot_bw_vs_m(df: pd.DataFrame, info: dict): + """ + Two y-axes: left = achieved GB/s, right = achieved GFLOP/s. + Dashed horizontal lines show hardware ceilings. + The cross-over from BW-saturation to FLOP-saturation is the ridge point. + """ + peak_bw = info["theoretical_bw_gbs"] + peak_flops = info["theoretical_fp32_gflops"] + emp_bw = info.get("empirical_bw_gbs", peak_bw) + + fig, ax1 = plt.subplots(figsize=(11, 6), layout="constrained") + ax2 = ax1.twinx() + + m = df["m"].values + + ln1, = ax1.semilogx(m, df["achieved_bw_gbs"], "o-", color="royalblue", + lw=2, ms=5, label="Achieved bandwidth") + ln2, = ax2.semilogx(m, df["achieved_gflops"], "s-", color="tomato", + lw=2, ms=5, label="Achieved GFLOP/s") + + ax1.axhline(emp_bw, color="royalblue", ls="--", lw=1.5, + label=f"Empirical peak BW ({emp_bw:.0f} GB/s)") + ax2.axhline(peak_flops, color="tomato", ls="--", lw=1.5, + label=f"Peak FP32 ({peak_flops:.0f} GFLOP/s)") + + # Regime annotations + ax1.axvspan(m[0], m[-1] // 4, alpha=0.05, color="royalblue") + ax1.axvspan(m[-1] // 4, m[-1], alpha=0.05, color="tomato") + ax1.text(1.5, emp_bw * 0.6, "← Memory-\nbound", fontsize=10, + color="royalblue", alpha=0.8) + ax1.text(m[-1] * 0.3, emp_bw * 0.6, "Compute-\nbound →", fontsize=10, + color="tomato", alpha=0.8) + + ax1.set_xlabel("m (FMA iterations per element, log scale)") + ax1.set_ylabel("Achieved Bandwidth [GB/s]", color="royalblue") + ax2.set_ylabel("Achieved Throughput [GFLOP/s]", color="tomato") + ax1.tick_params(axis="y", colors="royalblue") + ax2.tick_params(axis="y", colors="tomato") + + lines = [ln1, ln2] + labels = [ln.get_label() for ln in lines] + ax1.legend(lines, labels, loc="center left") + + device = info.get("device_name", "GPU") + ax1.set_title(f"Bandwidth & Compute Throughput vs Arithmetic Intensity — {device}") + ax1.grid(True, which="both", ls="--", lw=0.6, alpha=0.4) + + fig.savefig(RESULTS / "bw_vs_m.png", bbox_inches="tight") + plt.close(fig) + print("Saved: bw_vs_m.png") + + +# ───────────────────────────────────────────────────────────────────────────── +# Plot 3 — Bandwidth saturation vs vector size +# ───────────────────────────────────────────────────────────────────────────── + +def plot_bw_saturation(df: pd.DataFrame, info: dict): + """ + For m=1 (pure memory-bound kernel), shows that small N can't saturate + HBM bandwidth — you need enough concurrent memory requests. + """ + peak_bw = info["theoretical_bw_gbs"] + emp_bw = info.get("empirical_bw_gbs", peak_bw) + + fig, ax = plt.subplots(figsize=(10, 5), layout="constrained") + + ax.semilogx(df["N"], df["achieved_bw_gbs"], "o-", color="royalblue", + lw=2, ms=5, label="Achieved bandwidth") + ax.axhline(emp_bw, color="royalblue", ls="--", lw=1.5, + label=f"Empirical peak BW ({emp_bw:.0f} GB/s)") + ax.axhline(peak_bw, color="gray", ls=":", lw=1.2, + label=f"Theoretical peak BW ({peak_bw:.0f} GB/s)") + + ax.set_xlabel("Vector length N (log scale)") + ax.set_ylabel("Achieved Bandwidth [GB/s]") + device = info.get("device_name", "GPU") + ax.set_title(f"Memory Bandwidth Saturation (m=1) — {device}\n" + f"Bandwidth rises as N grows, saturating once enough warps fill the SM pipeline") + ax.legend() + ax.grid(True, which="both", ls="--", lw=0.6, alpha=0.5) + ax.xaxis.set_major_formatter(ticker.FuncFormatter( + lambda x, _: f"{int(x):,}")) + + fig.savefig(RESULTS / "bw_saturation.png", bbox_inches="tight") + plt.close(fig) + print("Saved: bw_saturation.png") + + +# ───────────────────────────────────────────────────────────────────────────── +# Plot 4 — Kernel time vs m with regime labels +# ───────────────────────────────────────────────────────────────────────────── + +def plot_time_vs_m(df: pd.DataFrame, info: dict): + """ + Log-log plot of kernel time vs m. + • Memory-bound regime: time is roughly constant (bottleneck is BW) + • Compute-bound regime: time ∝ m (slope ≈ 1 on log-log) + """ + fig, ax = plt.subplots(figsize=(10, 5), layout="constrained") + + m = df["m"].values + t = df["kernel_time_ms"].values + + ax.loglog(m, t, "o-", color="darkorchid", lw=2, ms=5) + + # Ideal memory-bound reference (flat line at t[0]) + ax.axhline(t[0], color="royalblue", ls="--", lw=1.5, alpha=0.7, + label=f"Memory-bound limit ({t[0]:.2f} ms)") + + # Ideal compute-bound reference: t ∝ m, anchored at last measured point + m_ref = m[-1] + t_ref = t[-1] + m_line = np.array([m[len(m)//3], m[-1]], dtype=float) + t_line = t_ref * (m_line / m_ref) + ax.loglog(m_line, t_line, color="tomato", ls="--", lw=1.5, alpha=0.7, + label="Compute-bound (slope 1)") + + ax.set_xlabel("m (FMA iterations per element)") + ax.set_ylabel("Kernel time [ms]") + device = info.get("device_name", "GPU") + ax.set_title(f"Kernel Time vs m — {device}\n" + "Flat = memory-bound; rising linearly (slope 1 on log-log) = compute-bound") + ax.legend() + ax.grid(True, which="both", ls="--", lw=0.6, alpha=0.5) + + fig.savefig(RESULTS / "time_vs_m.png", bbox_inches="tight") + plt.close(fig) + print("Saved: time_vs_m.png") + + +# ───────────────────────────────────────────────────────────────────────────── +# Plot 5 — Efficiency heatmap (bandwidth util + compute util vs m) +# ───────────────────────────────────────────────────────────────────────────── + +def plot_efficiency(df: pd.DataFrame, info: dict): + """ + Normalised utilisation: + bw_util = achieved_bw / empirical_peak_bw ∈ [0, 1] + flop_util = achieved_gflops / peak_fp32 ∈ [0, 1] + """ + emp_bw = info.get("empirical_bw_gbs", info["theoretical_bw_gbs"]) + peak_flops = info["theoretical_fp32_gflops"] + + bw_util = df["achieved_bw_gbs"] / emp_bw + flop_util = df["achieved_gflops"] / peak_flops + m = df["m"].values + + fig, ax = plt.subplots(figsize=(10, 5), layout="constrained") + + ax.semilogx(m, bw_util * 100, "o-", color="royalblue", lw=2, ms=5, + label="Memory BW utilisation (%)") + ax.semilogx(m, flop_util * 100, "s-", color="tomato", lw=2, ms=5, + label="FP32 compute utilisation (%)") + + ax.axhline(100, color="gray", ls="--", lw=1, alpha=0.6) + ax.set_ylim(0, 115) + ax.set_xlabel("m (FMA iterations per element)") + ax.set_ylabel("Hardware Utilisation [%]") + device = info.get("device_name", "GPU") + ax.set_title(f"Hardware Utilisation vs Arithmetic Intensity — {device}") + ax.legend() + ax.grid(True, which="both", ls="--", lw=0.6, alpha=0.4) + + fig.savefig(RESULTS / "efficiency.png", bbox_inches="tight") + plt.close(fig) + print("Saved: efficiency.png") + + +# ───────────────────────────────────────────────────────────────────────────── +# Main +# ───────────────────────────────────────────────────────────────────────────── + +def main(): + info = load_device_info(RESULTS / "device_info.csv") + df_roof = load_csv("roofline.csv") + df_bw_sat = load_csv("bw_saturation.csv") + + print(f"Device: {info.get('device_name', '?')}") + print(f"Theoretical BW: {info['theoretical_bw_gbs']:.0f} GB/s") + print(f"Empirical BW: {info.get('empirical_bw_gbs', float('nan')):.0f} GB/s") + print(f"Peak FP32: {info['theoretical_fp32_gflops']:.0f} GFLOP/s") + print(f"Ridge point: {info['ridge_point_flop_per_byte']:.1f} FLOP/byte\n") + + plot_roofline(df_roof, info) + plot_bw_vs_m(df_roof, info) + plot_time_vs_m(df_roof, info) + plot_efficiency(df_roof, info) + plot_bw_saturation(df_bw_sat, info) + + +if __name__ == "__main__": + main() diff --git a/cuda-memory-compute-bound/scripts/compile.sh b/cuda-memory-compute-bound/scripts/compile.sh new file mode 100755 index 0000000..a52b83f --- /dev/null +++ b/cuda-memory-compute-bound/scripts/compile.sh @@ -0,0 +1,41 @@ +#!/bin/bash +# Compile the memory-bound vs compute-bound benchmark +# Usage: ./compile.sh [sm_XX] (default: sm_80 for A100) +# ./compile.sh sm_86 (RTX 3090 / A10) +# ./compile.sh sm_90 (H100) + +set -e + +THIS_DIR=$(dirname "$(realpath "$0")") +ROOT_DIR=$(dirname "${THIS_DIR}") +ARCH="${1:-sm_80}" + +echo "Building for arch=${ARCH}" + +# ── Recreate build dirs ──────────────────────────────────────────────────── +rm -rf "${ROOT_DIR}/build" +mkdir -p "${ROOT_DIR}/build/obj" "${ROOT_DIR}/build/bin" + +INCLUDES="-I${ROOT_DIR}/include" +CUDA_INCLUDES="-I/usr/include" +CUDA_LIB_DIRS="-L/usr/lib/x86_64-linux-gnu/" +CUDA_LIB="-lcudart" +FLAGS="-O3 -arch=${ARCH} --use_fast_math" + +# ── Compile objects ──────────────────────────────────────────────────────── +nvcc ${FLAGS} ${INCLUDES} ${CUDA_INCLUDES} \ + -c "${ROOT_DIR}/src/kernels.cu" \ + -o "${ROOT_DIR}/build/obj/kernels.o" + +nvcc ${FLAGS} ${INCLUDES} ${CUDA_INCLUDES} \ + -c "${ROOT_DIR}/src/main.cu" \ + -o "${ROOT_DIR}/build/obj/main.o" + +# ── Link ─────────────────────────────────────────────────────────────────── +nvcc ${FLAGS} \ + "${ROOT_DIR}/build/obj/kernels.o" \ + "${ROOT_DIR}/build/obj/main.o" \ + ${CUDA_LIB_DIRS} ${CUDA_LIB} \ + -o "${ROOT_DIR}/build/bin/main" + +echo "Binary: ${ROOT_DIR}/build/bin/main" diff --git a/cuda-memory-compute-bound/scripts/execute.sh b/cuda-memory-compute-bound/scripts/execute.sh new file mode 100755 index 0000000..4689b4a --- /dev/null +++ b/cuda-memory-compute-bound/scripts/execute.sh @@ -0,0 +1,50 @@ +#!/bin/bash +# Run the benchmark and produce CSV results. +# Optional: pass --ncu to also profile with Nsight Compute. +# +# Usage: +# ./execute.sh # just benchmark +# ./execute.sh --ncu # benchmark + ncu profile of the m=1 and m=256 kernels + +set -e + +THIS_DIR=$(dirname "$(realpath "$0")") +ROOT_DIR=$(dirname "${THIS_DIR}") +BINARY="${ROOT_DIR}/build/bin/main" +RESULTS="${ROOT_DIR}/results" + +if [ ! -f "${BINARY}" ]; then + echo "Binary not found. Run scripts/compile.sh first." + exit 1 +fi + +mkdir -p "${RESULTS}" + +echo "=== Running benchmark ===" +"${BINARY}" \ + --results_dir "${RESULTS}" \ + --threads_per_block 256 \ + --N_roofline 33554432 \ + --bw_repeats 10 + +echo "" +echo "=== Results saved to ${RESULTS}/ ===" + +# ── Optional: Nsight Compute profiling ──────────────────────────────────── +if [[ "$1" == "--ncu" ]]; then + echo "" + echo "=== Nsight Compute profiling (m=1 and m=512) ===" + + # Profile the memory-bound case (m=1) + ncu --set full \ + --export "${RESULTS}/ncu_m1" \ + --force-overwrite \ + "${BINARY}" \ + --results_dir "${RESULTS}" \ + --threads_per_block 256 \ + --N_roofline 33554432 \ + --bw_repeats 1 2>/dev/null || \ + echo "ncu not available — skipping (install CUDA Toolkit for profiling)" + + echo "NCU reports: ${RESULTS}/ncu_m1.ncu-rep" +fi diff --git a/cuda-memory-compute-bound/src/kernels.cu b/cuda-memory-compute-bound/src/kernels.cu new file mode 100644 index 0000000..408b78a --- /dev/null +++ b/cuda-memory-compute-bound/src/kernels.cu @@ -0,0 +1,102 @@ +#include "kernels.h" + +// ───────────────────────────────────────────────────────────────────────────── +// Vector-add + FMA benchmark kernel +// +// For small m (low AI) → memory-bandwidth-bound +// For large m (high AI) → compute-throughput-bound +// +// Key: the GPU hides FMA latency (~4-cycle pipeline) by scheduling other warps +// while a given warp waits. With enough active warps (large N), the throughput +// asymptotically approaches 2 FMAs/cycle even with a single accumulator chain +// per thread. +// ───────────────────────────────────────────────────────────────────────────── +__global__ void vector_add_fma_kernel( + const float* __restrict__ a, + const float* __restrict__ b, + float* __restrict__ c, + int N, int m) +{ + int idx = blockIdx.x * blockDim.x + threadIdx.x; + if (idx >= N) return; + + // One memory-bound load+add (1 FLOP, 8 bytes read) + float val = a[idx] + b[idx]; + + // m compute-intensive FMAs: val = val * alpha + beta (2 FLOPs each) + // alpha/beta chosen so val stays bounded for any starting value in [0,2]: + // val → converges to beta/(1-alpha) = 0.0001/(0.00001) = 10 + // #pragma unroll 1 prevents the compiler from unrolling (keeps m general) + #pragma unroll 1 + for (int j = 0; j < m; ++j) { + val = fmaf(val, 1.0f - 1e-5f, 1e-4f); + } + + // One write (4 bytes written) + c[idx] = val; +} + +// ───────────────────────────────────────────────────────────────────────────── +// Pure-copy kernel — measures empirical peak memory bandwidth +// ───────────────────────────────────────────────────────────────────────────── +__global__ void bandwidth_test_kernel( + const float* __restrict__ a, + float* __restrict__ b, + int N) +{ + int idx = blockIdx.x * blockDim.x + threadIdx.x; + if (idx < N) b[idx] = a[idx]; +} + +// ───────────────────────────────────────────────────────────────────────────── +// Tiny kernel just to force driver/context init before benchmarking +// ───────────────────────────────────────────────────────────────────────────── +__global__ void warmup_kernel(float* x, int N) +{ + int idx = blockIdx.x * blockDim.x + threadIdx.x; + if (idx < N) x[idx] = static_cast(idx); +} + +// ───────────────────────────────────────────────────────────────────────────── +// Host wrappers +// ───────────────────────────────────────────────────────────────────────────── + +void vector_add_fma( + const float* d_a, const float* d_b, float* d_c, + int N, int m, int threads_per_block, + cudaEvent_t start, cudaEvent_t stop, float& elapsed_ms) +{ + int blocks = (N + threads_per_block - 1) / threads_per_block; + + CHECK_CUDA_ERROR(cudaEventRecord(start)); + vector_add_fma_kernel<<>>(d_a, d_b, d_c, N, m); + CHECK_CUDA_ERROR(cudaEventRecord(stop)); + CHECK_LAST_CUDA_ERROR(); + CHECK_CUDA_ERROR(cudaEventSynchronize(stop)); + CHECK_CUDA_ERROR(cudaEventElapsedTime(&elapsed_ms, start, stop)); +} + +void bandwidth_test( + const float* d_a, float* d_b, + int N, int threads_per_block, + cudaEvent_t start, cudaEvent_t stop, float& elapsed_ms) +{ + int blocks = (N + threads_per_block - 1) / threads_per_block; + + CHECK_CUDA_ERROR(cudaEventRecord(start)); + bandwidth_test_kernel<<>>(d_a, d_b, N); + CHECK_CUDA_ERROR(cudaEventRecord(stop)); + CHECK_LAST_CUDA_ERROR(); + CHECK_CUDA_ERROR(cudaEventSynchronize(stop)); + CHECK_CUDA_ERROR(cudaEventElapsedTime(&elapsed_ms, start, stop)); +} + +void warmup() +{ + float* d; + CHECK_CUDA_ERROR(cudaMalloc(&d, 4096 * sizeof(float))); + warmup_kernel<<<16, 256>>>(d, 4096); + CHECK_LAST_CUDA_ERROR(); + CHECK_CUDA_ERROR(cudaDeviceSynchronize()); + CHECK_CUDA_ERROR(cudaFree(d)); +} diff --git a/cuda-memory-compute-bound/src/main.cu b/cuda-memory-compute-bound/src/main.cu new file mode 100644 index 0000000..11a9525 --- /dev/null +++ b/cuda-memory-compute-bound/src/main.cu @@ -0,0 +1,356 @@ +// +// main.cu — Memory-Bound vs Compute-Bound: A GPU Deep Dive with Vector Addition +// +// Benchmarks two related quantities as the arithmetic intensity (AI) rises: +// • Achieved memory bandwidth (GB/s) +// • Achieved FP32 throughput (GFLOP/s) +// +// The kernel does: +// val = a[i] + b[i] ← 1 FLOP, 8 bytes read +// for j in 0..m: val = fma(val,α,β) ← 2m FLOPs, no extra memory +// c[i] = val ← 0 FLOPs, 4 bytes written +// +// AI = (1 + 2m) / 12 [FLOPs / byte] +// +// Sweeps: +// 1. ROOFLINE — fixed large N, m ∈ {1,2,4,…,16384} +// Produces the roofline-model scatter plot. +// 2. BW_SAT — fixed m=1, N ∈ {2^10 … 2^27} +// Shows how bandwidth saturates as the working set grows. +// +// Output: two CSV files in / +// roofline.csv +// bw_saturation.csv +// +// Usage: +// ./main --results_dir ../results [--threads_per_block 256] +// [--N_roofline 33554432] [--bw_repeats 10] + +#include "kernels.h" +#include "utils.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// ───────────────────────────────────────────────────────────────────────────── +// Theoretical peak values from device properties +// ───────────────────────────────────────────────────────────────────────────── + +// CUDA cores per SM keyed on (major, minor) compute capability +static int cores_per_sm(int major, int minor) +{ + struct Entry { int major, minor, cores; }; + static constexpr Entry table[] = { + {3, 0, 192}, {3, 2, 192}, {3, 5, 192}, {3, 7, 192}, + {5, 0, 128}, {5, 2, 128}, + {6, 0, 64}, {6, 1, 128}, {6, 2, 128}, + {7, 0, 64}, {7, 2, 64}, {7, 5, 64}, + {8, 0, 64}, {8, 6, 128}, {8, 7, 128}, {8, 9, 128}, + {9, 0, 128}, + }; + for (auto& e : table) + if (e.major == major && e.minor == minor) return e.cores; + return 128; // safe default +} + +// Theoretical peak HBM / GDDR bandwidth in GB/s +static double theoretical_bw_gbs(const cudaDeviceProp& p) +{ + // memoryClockRate is in kHz; memoryBusWidth is in bits + // Bandwidth = clock_Hz × (bus_width_bits / 8) × 2 (DDR) / 1e9 + return static_cast(p.memoryClockRate) * 1e3 // → Hz + * (p.memoryBusWidth / 8.0) // → bytes/cycle + * 2.0 // DDR + / 1e9; +} + +// Theoretical peak FP32 in GFLOP/s +static double theoretical_fp32_gflops(const cudaDeviceProp& p) +{ + // Each CUDA core can issue 1 FMA (= 2 FLOPs) per clock cycle + return 2.0 + * p.multiProcessorCount + * cores_per_sm(p.major, p.minor) + * static_cast(p.clockRate) * 1e3 // → Hz + / 1e9; +} + +// ───────────────────────────────────────────────────────────────────────────── +// Helpers +// ───────────────────────────────────────────────────────────────────────────── + +static float time_kernel_ms( + const float* d_a, const float* d_b, float* d_c, + int N, int m, int tpb, int repeats, + cudaEvent_t ev0, cudaEvent_t ev1) +{ + float best = 1e30f; + for (int r = 0; r < repeats; ++r) { + float t; + vector_add_fma(d_a, d_b, d_c, N, m, tpb, ev0, ev1, t); + best = std::min(best, t); + } + return best; +} + +static float time_bw_kernel_ms( + const float* d_a, float* d_c, + int N, int tpb, int repeats, + cudaEvent_t ev0, cudaEvent_t ev1) +{ + float best = 1e30f; + for (int r = 0; r < repeats; ++r) { + float t; + bandwidth_test(d_a, d_c, N, tpb, ev0, ev1, t); + best = std::min(best, t); + } + return best; +} + +// ───────────────────────────────────────────────────────────────────────────── +// Arg parsing +// ───────────────────────────────────────────────────────────────────────────── + +struct Config { + std::string results_dir = "../results"; + int threads_per_block = 256; + long N_roofline = 1L << 25; // 33 554 432 floats = 128 MB / array + int bw_repeats = 10; +}; + +static Config parse_args(int argc, char** argv) +{ + Config cfg; + for (int i = 1; i < argc; ++i) { + std::string a = argv[i]; + if (a == "--results_dir" && i + 1 < argc) cfg.results_dir = argv[++i]; + if (a == "--threads_per_block"&& i + 1 < argc) cfg.threads_per_block = std::stoi(argv[++i]); + if (a == "--N_roofline" && i + 1 < argc) cfg.N_roofline = std::stol(argv[++i]); + if (a == "--bw_repeats" && i + 1 < argc) cfg.bw_repeats = std::stoi(argv[++i]); + if (a == "--help") { + std::cout << "Usage: ./main [--results_dir DIR] [--threads_per_block INT]\n" + << " [--N_roofline INT] [--bw_repeats INT]\n"; + std::exit(0); + } + } + return cfg; +} + +// ───────────────────────────────────────────────────────────────────────────── +// main +// ───────────────────────────────────────────────────────────────────────────── + +int main(int argc, char** argv) +{ + Config cfg = parse_args(argc, argv); + + // ── Device info ───────────────────────────────────────────────────────── + cudaDeviceProp prop; + CHECK_CUDA_ERROR(cudaGetDeviceProperties(&prop, 0)); + + double th_bw = theoretical_bw_gbs(prop); + double th_flops = theoretical_fp32_gflops(prop); + double ridge_pt = th_flops / th_bw; // FLOPs/byte at the ridge + + std::cout << "Device : " << prop.name << "\n" + << "Compute cap. : " << prop.major << "." << prop.minor << "\n" + << "SMs : " << prop.multiProcessorCount << "\n" + << "CUDA cores/SM : " << cores_per_sm(prop.major, prop.minor) << "\n" + << "GPU clock : " << prop.clockRate / 1e6 << " GHz\n" + << "Mem clock : " << prop.memoryClockRate / 1e6 << " GHz\n" + << "Bus width : " << prop.memoryBusWidth << " bits\n" + << "Theoretical BW : " << th_bw << " GB/s\n" + << "Theoretical FLOPS: " << th_flops << " GFLOP/s\n" + << "Ridge point : " << ridge_pt << " FLOP/byte\n\n"; + + // ── Make results directory ─────────────────────────────────────────────── + std::filesystem::create_directories(cfg.results_dir); + + // ── Shared CUDA events ─────────────────────────────────────────────────── + cudaEvent_t ev0, ev1; + CHECK_CUDA_ERROR(cudaEventCreate(&ev0)); + CHECK_CUDA_ERROR(cudaEventCreate(&ev1)); + + // ── Warm up ────────────────────────────────────────────────────────────── + warmup(); + warmup(); + + // ── Measure empirical peak bandwidth (pure-copy kernel) ────────────────── + { + const long N_bw = cfg.N_roofline; + float *d_a, *d_b; + CHECK_CUDA_ERROR(cudaMalloc(&d_a, N_bw * sizeof(float))); + CHECK_CUDA_ERROR(cudaMalloc(&d_b, N_bw * sizeof(float))); + CHECK_CUDA_ERROR(cudaMemset(d_a, 1, N_bw * sizeof(float))); + + float t = time_bw_kernel_ms(d_a, d_b, static_cast(N_bw), + cfg.threads_per_block, cfg.bw_repeats, ev0, ev1); + + // copy reads 1 array + writes 1 array = 2 × N × sizeof(float) bytes + double bytes = 2.0 * N_bw * sizeof(float); + double emp_bw = bytes / (t * 1e-3) / 1e9; + + std::cout << "Empirical peak BW: " << emp_bw << " GB/s (vs " + << th_bw << " GB/s theoretical)\n\n"; + + // Save device info + measured peak to a small metadata file + std::ofstream meta(cfg.results_dir + "/device_info.csv"); + meta << "key,value\n" + << "device_name," << prop.name << "\n" + << "compute_capability," << prop.major << "." << prop.minor << "\n" + << "num_sms," << prop.multiProcessorCount << "\n" + << "cores_per_sm," << cores_per_sm(prop.major, prop.minor) << "\n" + << "gpu_clock_ghz," << prop.clockRate / 1e6 << "\n" + << "mem_clock_ghz," << prop.memoryClockRate / 1e6 << "\n" + << "mem_bus_width_bits," << prop.memoryBusWidth << "\n" + << "theoretical_bw_gbs," << th_bw << "\n" + << "theoretical_fp32_gflops," << th_flops << "\n" + << "ridge_point_flop_per_byte," << ridge_pt << "\n" + << "empirical_bw_gbs," << emp_bw << "\n"; + + CHECK_CUDA_ERROR(cudaFree(d_a)); + CHECK_CUDA_ERROR(cudaFree(d_b)); + } + + // ═════════════════════════════════════════════════════════════════════════ + // SWEEP 1: ROOFLINE — fixed N, sweep m + // + // For each m we record: + // • kernel_time_ms + // • achieved_bandwidth_gbs = 12·N / (time_s · 1e9) + // • achieved_gflops = (1+2m)·N / (time_s · 1e9) + // • arithmetic_intensity = (1+2m) / 12 + // ═════════════════════════════════════════════════════════════════════════ + { + const long N = cfg.N_roofline; + const double bytes_per_elem = 3.0 * sizeof(float); // 2 reads + 1 write + + // m values span several orders of magnitude to cross the ridge point + std::vector m_values; + for (int m = 1; m <= 16384; m *= 2) m_values.push_back(m); + // add fine-grained points near the expected ridge + for (int m : {3, 6, 12, 24, 48, 96, 192, 384, 768, 1536, 3072, 6144, 12288}) + m_values.push_back(m); + std::sort(m_values.begin(), m_values.end()); + m_values.erase(std::unique(m_values.begin(), m_values.end()), m_values.end()); + + float *d_a, *d_b, *d_c; + CHECK_CUDA_ERROR(cudaMalloc(&d_a, N * sizeof(float))); + CHECK_CUDA_ERROR(cudaMalloc(&d_b, N * sizeof(float))); + CHECK_CUDA_ERROR(cudaMalloc(&d_c, N * sizeof(float))); + + // Fill with random data (on host, then copy) + std::vector h(N); + std::mt19937 rng(42); + std::uniform_real_distribution dist(0.0f, 1.0f); + std::generate(h.begin(), h.end(), [&]{ return dist(rng); }); + CHECK_CUDA_ERROR(cudaMemcpy(d_a, h.data(), N * sizeof(float), cudaMemcpyHostToDevice)); + CHECK_CUDA_ERROR(cudaMemcpy(d_b, h.data(), N * sizeof(float), cudaMemcpyHostToDevice)); + + std::ofstream csv(cfg.results_dir + "/roofline.csv"); + csv << "m,N,kernel_time_ms,arithmetic_intensity,achieved_bw_gbs,achieved_gflops," + "theoretical_bw_gbs,theoretical_fp32_gflops,ridge_point\n"; + + std::cout << "=== ROOFLINE SWEEP (N=" << N << ") ===\n"; + std::cout << " m AI(FLOPs/B) BW(GB/s) GFLOP/s\n"; + + for (int m : m_values) { + float t_ms = time_kernel_ms(d_a, d_b, d_c, static_cast(N), + m, cfg.threads_per_block, + cfg.bw_repeats, ev0, ev1); + + double t_s = t_ms * 1e-3; + double flops = static_cast(N) * (1.0 + 2.0 * m); + double bytes = static_cast(N) * bytes_per_elem; + double ai = flops / bytes; + double bw_gbs = bytes / t_s / 1e9; + double gflops = flops / t_s / 1e9; + + csv << m << "," << N << "," << t_ms << "," + << ai << "," << bw_gbs << "," << gflops << "," + << th_bw << "," << th_flops << "," << ridge_pt << "\n"; + + std::cout << " m=" << m + << " AI=" << ai + << " BW=" << bw_gbs << " GB/s" + << " GFLOP/s=" << gflops << "\n"; + } + + csv.close(); + CHECK_CUDA_ERROR(cudaFree(d_a)); + CHECK_CUDA_ERROR(cudaFree(d_b)); + CHECK_CUDA_ERROR(cudaFree(d_c)); + std::cout << "\n"; + } + + // ═════════════════════════════════════════════════════════════════════════ + // SWEEP 2: BANDWIDTH SATURATION — m=1, sweep N + // + // Shows that small vectors can't saturate HBM bandwidth (the GPU doesn't + // have enough warps in flight to fill the memory controllers), while large + // vectors approach the empirical peak. + // ═════════════════════════════════════════════════════════════════════════ + { + const int m = 1; // minimal compute — purely memory-bound + const double bytes_per_elem = 3.0 * sizeof(float); + + // N from 1024 to N_roofline in powers of 2 + std::vector n_values; + for (long n = 1024; n <= cfg.N_roofline; n *= 2) + n_values.push_back(n); + + std::ofstream csv(cfg.results_dir + "/bw_saturation.csv"); + csv << "m,N,sizeMB,kernel_time_ms,arithmetic_intensity," + "achieved_bw_gbs,achieved_gflops," + "theoretical_bw_gbs,theoretical_fp32_gflops\n"; + + std::cout << "=== BANDWIDTH SATURATION SWEEP (m=" << m << ") ===\n"; + std::cout << " N sizeMB BW(GB/s)\n"; + + for (long N : n_values) { + float *d_a, *d_b, *d_c; + CHECK_CUDA_ERROR(cudaMalloc(&d_a, N * sizeof(float))); + CHECK_CUDA_ERROR(cudaMalloc(&d_b, N * sizeof(float))); + CHECK_CUDA_ERROR(cudaMalloc(&d_c, N * sizeof(float))); + CHECK_CUDA_ERROR(cudaMemset(d_a, 1, N * sizeof(float))); + CHECK_CUDA_ERROR(cudaMemset(d_b, 2, N * sizeof(float))); + + float t_ms = time_kernel_ms(d_a, d_b, d_c, static_cast(N), + m, cfg.threads_per_block, + cfg.bw_repeats, ev0, ev1); + + double t_s = t_ms * 1e-3; + double flops = static_cast(N) * (1.0 + 2.0 * m); + double bytes = static_cast(N) * bytes_per_elem; + double ai = flops / bytes; + double bw_gbs = bytes / t_s / 1e9; + double gflops = flops / t_s / 1e9; + double sizeMB = static_cast(N) * sizeof(float) / (1024.0 * 1024.0); + + csv << m << "," << N << "," << sizeMB << "," << t_ms << "," + << ai << "," << bw_gbs << "," << gflops << "," + << th_bw << "," << th_flops << "\n"; + + std::cout << " N=" << N << " " << sizeMB << " MB BW=" << bw_gbs << " GB/s\n"; + + CHECK_CUDA_ERROR(cudaFree(d_a)); + CHECK_CUDA_ERROR(cudaFree(d_b)); + CHECK_CUDA_ERROR(cudaFree(d_c)); + } + + csv.close(); + } + + CHECK_CUDA_ERROR(cudaEventDestroy(ev0)); + CHECK_CUDA_ERROR(cudaEventDestroy(ev1)); + + std::cout << "\nResults written to: " << cfg.results_dir << "/\n"; + return 0; +} From 89940480ddac64df0bcf01f017b114d77615fced Mon Sep 17 00:00:00 2001 From: Sebas <12949040+SebastiaAgramunt@users.noreply.github.com> Date: Sat, 18 Apr 2026 03:16:58 +0000 Subject: [PATCH 02/21] kernel performance --- cuda-memory-compute-bound/include/kernels.h | 43 --- cuda-memory-compute-bound/include/utils.h | 29 -- cuda-memory-compute-bound/scripts/analyze.py | 326 ----------------- cuda-memory-compute-bound/scripts/compile.sh | 41 --- cuda-memory-compute-bound/scripts/execute.sh | 50 --- cuda-memory-compute-bound/src/kernels.cu | 102 ------ cuda-memory-compute-bound/src/main.cu | 356 ------------------- cuda-mma/Makefile | 45 +++ cuda-mma/include/cuda_check.cuh | 40 +++ cuda-mma/include/kernels.cuh | 38 ++ cuda-mma/include/timer.h | 29 ++ cuda-mma/include/utils.h | 17 + cuda-mma/main.cu | 242 +++++++++++++ cuda-mma/src/kernels.cu | 108 ++++++ cuda-mma/src/utils.cu | 38 ++ 15 files changed, 557 insertions(+), 947 deletions(-) delete mode 100644 cuda-memory-compute-bound/include/kernels.h delete mode 100644 cuda-memory-compute-bound/include/utils.h delete mode 100755 cuda-memory-compute-bound/scripts/analyze.py delete mode 100755 cuda-memory-compute-bound/scripts/compile.sh delete mode 100755 cuda-memory-compute-bound/scripts/execute.sh delete mode 100644 cuda-memory-compute-bound/src/kernels.cu delete mode 100644 cuda-memory-compute-bound/src/main.cu create mode 100644 cuda-mma/Makefile create mode 100644 cuda-mma/include/cuda_check.cuh create mode 100644 cuda-mma/include/kernels.cuh create mode 100644 cuda-mma/include/timer.h create mode 100644 cuda-mma/include/utils.h create mode 100644 cuda-mma/main.cu create mode 100644 cuda-mma/src/kernels.cu create mode 100644 cuda-mma/src/utils.cu diff --git a/cuda-memory-compute-bound/include/kernels.h b/cuda-memory-compute-bound/include/kernels.h deleted file mode 100644 index a1b7f71..0000000 --- a/cuda-memory-compute-bound/include/kernels.h +++ /dev/null @@ -1,43 +0,0 @@ -#ifndef KERNELS_H -#define KERNELS_H - -#include "utils.h" -#include - -// ─── Benchmark kernel ──────────────────────────────────────────────────────── -// -// Each thread: -// 1. Loads a[i] and b[i] → 2 × sizeof(float) = 8 bytes read -// 2. Computes val = a[i] + b[i] (1 FLOP) -// 3. Runs m FMA iterations: val = val * alpha + beta (2 FLOPs each) -// 4. Stores c[i] → 1 × sizeof(float) = 4 bytes written -// -// Totals per element -// Bytes: 3 × 4 = 12 -// FLOPs: 1 + 2·m -// Arithmetic intensity: (1 + 2·m) / 12 [FLOPs / byte] -// -// The GPU hides FMA latency (~4 cycles) through inter-warp parallelism, so for -// large N the kernel will approach peak FP32 throughput once AI exceeds the -// ridge point. -void vector_add_fma( - const float* d_a, const float* d_b, float* d_c, - int N, int m, int threads_per_block, - cudaEvent_t start, cudaEvent_t stop, float& elapsed_ms); - -// ─── Peak-bandwidth kernel ──────────────────────────────────────────────────── -// -// Pure memory copy: c[i] = a[i] -// Bytes: 2 × 4 = 8 per element -// FLOPs: 0 → arithmetic intensity ≈ 0 -// -// Used to measure empirical peak HBM / GDDR bandwidth. -void bandwidth_test( - const float* d_a, float* d_b, - int N, int threads_per_block, - cudaEvent_t start, cudaEvent_t stop, float& elapsed_ms); - -// GPU warm-up (first kernel launch incurs driver-init overhead) -void warmup(); - -#endif diff --git a/cuda-memory-compute-bound/include/utils.h b/cuda-memory-compute-bound/include/utils.h deleted file mode 100644 index 3e70ca5..0000000 --- a/cuda-memory-compute-bound/include/utils.h +++ /dev/null @@ -1,29 +0,0 @@ -#ifndef UTILS_H -#define UTILS_H - -#include -#include - -inline void check(cudaError_t err, const char* const func, const char* const file, const int line) -{ - if (err != cudaSuccess) { - std::cerr << "CUDA Runtime Error at: " << file << ":" << line << "\n"; - std::cerr << cudaGetErrorString(err) << " " << func << "\n"; - std::exit(EXIT_FAILURE); - } -} - -inline void checkLast(const char* const file, const int line) -{ - cudaError_t const err{cudaGetLastError()}; - if (err != cudaSuccess) { - std::cerr << "CUDA Runtime Error at: " << file << ":" << line << "\n"; - std::cerr << cudaGetErrorString(err) << "\n"; - std::exit(EXIT_FAILURE); - } -} - -#define CHECK_CUDA_ERROR(expr) (check((expr), #expr, __FILE__, __LINE__)) -#define CHECK_LAST_CUDA_ERROR() (checkLast(__FILE__, __LINE__)) - -#endif diff --git a/cuda-memory-compute-bound/scripts/analyze.py b/cuda-memory-compute-bound/scripts/analyze.py deleted file mode 100755 index f0a01ee..0000000 --- a/cuda-memory-compute-bound/scripts/analyze.py +++ /dev/null @@ -1,326 +0,0 @@ -#!/usr/bin/env python3 -""" -analyze.py — Visualise memory-bound vs compute-bound GPU behaviour. - -Reads: - results/roofline.csv - results/bw_saturation.csv - results/device_info.csv - -Produces: - results/roofline.png — classic roofline model - results/bw_vs_m.png — bandwidth & GFLOP/s vs m - results/bw_saturation.png — bandwidth vs vector size (N) - results/time_breakdown.png — compute time vs m with regime labels -""" - -import pathlib -import numpy as np -import pandas as pd -import matplotlib.pyplot as plt -import matplotlib.ticker as ticker - -THIS_DIR = pathlib.Path(__file__).parent.resolve() -ROOT_DIR = THIS_DIR.parent.resolve() -RESULTS = ROOT_DIR / "results" - -plt.rcParams.update({ - "font.size": 13, - "axes.titlesize": 15, - "axes.labelsize": 13, - "legend.fontsize": 11, - "figure.dpi": 120, -}) - -# ───────────────────────────────────────────────────────────────────────────── -# Load data -# ───────────────────────────────────────────────────────────────────────────── - -def load_device_info(path: pathlib.Path) -> dict: - df = pd.read_csv(path, index_col="key") - return {k: float(v) if v.replace(".", "", 1).isdigit() else v - for k, v in df["value"].items()} - - -def load_csv(name: str) -> pd.DataFrame: - path = RESULTS / name - if not path.exists(): - raise FileNotFoundError(f"Missing {path} — run execute.sh first") - return pd.read_csv(path) - - -# ───────────────────────────────────────────────────────────────────────────── -# Plot 1 — Roofline model -# ───────────────────────────────────────────────────────────────────────────── - -def plot_roofline(df: pd.DataFrame, info: dict): - """ - Classic roofline: GFLOP/s on y-axis, arithmetic intensity on x-axis. - - Two hardware ceilings: - • Memory-bandwidth ceiling: GFLOP/s = BW_GB/s × AI - • Compute ceiling: GFLOP/s = peak_FP32 - - Points below the ridge are memory-bound; above are compute-bound. - """ - peak_bw = info["theoretical_bw_gbs"] # GB/s - peak_flops = info["theoretical_fp32_gflops"] # GFLOP/s - emp_bw = info.get("empirical_bw_gbs", peak_bw) - ridge_pt = info["ridge_point_flop_per_byte"] - - ai_range = np.logspace(-2, 4, 500) - - # Roofline ceiling - roof_mem = peak_bw * ai_range # memory-bandwidth ceiling (GFLOP/s) - roof_compute = np.full_like(ai_range, peak_flops) - roofline = np.minimum(roof_mem, roof_compute) - - fig, ax = plt.subplots(figsize=(10, 6), layout="constrained") - - ax.loglog(ai_range, roofline, "k-", lw=2.5, label="Roofline (theoretical)", zorder=3) - ax.loglog(ai_range, emp_bw * ai_range, "k--", lw=1.5, - label=f"Empirical BW ({emp_bw:.0f} GB/s)", zorder=3) - ax.axhline(peak_flops, color="steelblue", lw=1.5, ls=":", - label=f"Peak FP32 ({peak_flops:.0f} GFLOP/s)") - - # Mark ridge point - ax.axvline(ridge_pt, color="gray", lw=1.2, ls="--") - ax.text(ridge_pt * 1.08, peak_flops * 0.55, f"Ridge\n{ridge_pt:.1f} FLOP/byte", - fontsize=10, color="gray", va="center") - - # Shade regimes - ax.axvspan(ai_range[0], ridge_pt, alpha=0.06, color="royalblue", - label="Memory-bound region") - ax.axvspan(ridge_pt, ai_range[-1], alpha=0.06, color="tomato", - label="Compute-bound region") - - # Measured points, coloured by m - m_vals = df["m"].values - sc = ax.scatter(df["arithmetic_intensity"], df["achieved_gflops"], - c=np.log2(m_vals.astype(float)), cmap="plasma", - s=70, zorder=5, edgecolors="k", linewidths=0.4) - - # Annotate a few interesting m values - for m_label in [1, 32, 128, 1024, 8192]: - row = df[df["m"] == m_label] - if row.empty: - continue - ax.annotate(f"m={m_label}", - xy=(row["arithmetic_intensity"].values[0], - row["achieved_gflops"].values[0]), - xytext=(6, 4), textcoords="offset points", - fontsize=8, color="black") - - cbar = fig.colorbar(sc, ax=ax, pad=0.02) - cbar.set_label("log₂(m) — FMA iterations per element") - - device = info.get("device_name", "GPU") - ax.set_xlabel("Arithmetic Intensity [FLOP / byte]") - ax.set_ylabel("Achieved Throughput [GFLOP/s]") - ax.set_title(f"Roofline Model — {device}\nvector_add + m FMAs per element") - ax.legend(loc="lower right") - ax.grid(True, which="both", ls="--", lw=0.6, alpha=0.5) - ax.set_xlim(0.05, 5000) - ax.set_ylim(1, peak_flops * 3) - - fig.savefig(RESULTS / "roofline.png", bbox_inches="tight") - plt.close(fig) - print("Saved: roofline.png") - - -# ───────────────────────────────────────────────────────────────────────────── -# Plot 2 — Bandwidth and GFLOP/s vs m -# ───────────────────────────────────────────────────────────────────────────── - -def plot_bw_vs_m(df: pd.DataFrame, info: dict): - """ - Two y-axes: left = achieved GB/s, right = achieved GFLOP/s. - Dashed horizontal lines show hardware ceilings. - The cross-over from BW-saturation to FLOP-saturation is the ridge point. - """ - peak_bw = info["theoretical_bw_gbs"] - peak_flops = info["theoretical_fp32_gflops"] - emp_bw = info.get("empirical_bw_gbs", peak_bw) - - fig, ax1 = plt.subplots(figsize=(11, 6), layout="constrained") - ax2 = ax1.twinx() - - m = df["m"].values - - ln1, = ax1.semilogx(m, df["achieved_bw_gbs"], "o-", color="royalblue", - lw=2, ms=5, label="Achieved bandwidth") - ln2, = ax2.semilogx(m, df["achieved_gflops"], "s-", color="tomato", - lw=2, ms=5, label="Achieved GFLOP/s") - - ax1.axhline(emp_bw, color="royalblue", ls="--", lw=1.5, - label=f"Empirical peak BW ({emp_bw:.0f} GB/s)") - ax2.axhline(peak_flops, color="tomato", ls="--", lw=1.5, - label=f"Peak FP32 ({peak_flops:.0f} GFLOP/s)") - - # Regime annotations - ax1.axvspan(m[0], m[-1] // 4, alpha=0.05, color="royalblue") - ax1.axvspan(m[-1] // 4, m[-1], alpha=0.05, color="tomato") - ax1.text(1.5, emp_bw * 0.6, "← Memory-\nbound", fontsize=10, - color="royalblue", alpha=0.8) - ax1.text(m[-1] * 0.3, emp_bw * 0.6, "Compute-\nbound →", fontsize=10, - color="tomato", alpha=0.8) - - ax1.set_xlabel("m (FMA iterations per element, log scale)") - ax1.set_ylabel("Achieved Bandwidth [GB/s]", color="royalblue") - ax2.set_ylabel("Achieved Throughput [GFLOP/s]", color="tomato") - ax1.tick_params(axis="y", colors="royalblue") - ax2.tick_params(axis="y", colors="tomato") - - lines = [ln1, ln2] - labels = [ln.get_label() for ln in lines] - ax1.legend(lines, labels, loc="center left") - - device = info.get("device_name", "GPU") - ax1.set_title(f"Bandwidth & Compute Throughput vs Arithmetic Intensity — {device}") - ax1.grid(True, which="both", ls="--", lw=0.6, alpha=0.4) - - fig.savefig(RESULTS / "bw_vs_m.png", bbox_inches="tight") - plt.close(fig) - print("Saved: bw_vs_m.png") - - -# ───────────────────────────────────────────────────────────────────────────── -# Plot 3 — Bandwidth saturation vs vector size -# ───────────────────────────────────────────────────────────────────────────── - -def plot_bw_saturation(df: pd.DataFrame, info: dict): - """ - For m=1 (pure memory-bound kernel), shows that small N can't saturate - HBM bandwidth — you need enough concurrent memory requests. - """ - peak_bw = info["theoretical_bw_gbs"] - emp_bw = info.get("empirical_bw_gbs", peak_bw) - - fig, ax = plt.subplots(figsize=(10, 5), layout="constrained") - - ax.semilogx(df["N"], df["achieved_bw_gbs"], "o-", color="royalblue", - lw=2, ms=5, label="Achieved bandwidth") - ax.axhline(emp_bw, color="royalblue", ls="--", lw=1.5, - label=f"Empirical peak BW ({emp_bw:.0f} GB/s)") - ax.axhline(peak_bw, color="gray", ls=":", lw=1.2, - label=f"Theoretical peak BW ({peak_bw:.0f} GB/s)") - - ax.set_xlabel("Vector length N (log scale)") - ax.set_ylabel("Achieved Bandwidth [GB/s]") - device = info.get("device_name", "GPU") - ax.set_title(f"Memory Bandwidth Saturation (m=1) — {device}\n" - f"Bandwidth rises as N grows, saturating once enough warps fill the SM pipeline") - ax.legend() - ax.grid(True, which="both", ls="--", lw=0.6, alpha=0.5) - ax.xaxis.set_major_formatter(ticker.FuncFormatter( - lambda x, _: f"{int(x):,}")) - - fig.savefig(RESULTS / "bw_saturation.png", bbox_inches="tight") - plt.close(fig) - print("Saved: bw_saturation.png") - - -# ───────────────────────────────────────────────────────────────────────────── -# Plot 4 — Kernel time vs m with regime labels -# ───────────────────────────────────────────────────────────────────────────── - -def plot_time_vs_m(df: pd.DataFrame, info: dict): - """ - Log-log plot of kernel time vs m. - • Memory-bound regime: time is roughly constant (bottleneck is BW) - • Compute-bound regime: time ∝ m (slope ≈ 1 on log-log) - """ - fig, ax = plt.subplots(figsize=(10, 5), layout="constrained") - - m = df["m"].values - t = df["kernel_time_ms"].values - - ax.loglog(m, t, "o-", color="darkorchid", lw=2, ms=5) - - # Ideal memory-bound reference (flat line at t[0]) - ax.axhline(t[0], color="royalblue", ls="--", lw=1.5, alpha=0.7, - label=f"Memory-bound limit ({t[0]:.2f} ms)") - - # Ideal compute-bound reference: t ∝ m, anchored at last measured point - m_ref = m[-1] - t_ref = t[-1] - m_line = np.array([m[len(m)//3], m[-1]], dtype=float) - t_line = t_ref * (m_line / m_ref) - ax.loglog(m_line, t_line, color="tomato", ls="--", lw=1.5, alpha=0.7, - label="Compute-bound (slope 1)") - - ax.set_xlabel("m (FMA iterations per element)") - ax.set_ylabel("Kernel time [ms]") - device = info.get("device_name", "GPU") - ax.set_title(f"Kernel Time vs m — {device}\n" - "Flat = memory-bound; rising linearly (slope 1 on log-log) = compute-bound") - ax.legend() - ax.grid(True, which="both", ls="--", lw=0.6, alpha=0.5) - - fig.savefig(RESULTS / "time_vs_m.png", bbox_inches="tight") - plt.close(fig) - print("Saved: time_vs_m.png") - - -# ───────────────────────────────────────────────────────────────────────────── -# Plot 5 — Efficiency heatmap (bandwidth util + compute util vs m) -# ───────────────────────────────────────────────────────────────────────────── - -def plot_efficiency(df: pd.DataFrame, info: dict): - """ - Normalised utilisation: - bw_util = achieved_bw / empirical_peak_bw ∈ [0, 1] - flop_util = achieved_gflops / peak_fp32 ∈ [0, 1] - """ - emp_bw = info.get("empirical_bw_gbs", info["theoretical_bw_gbs"]) - peak_flops = info["theoretical_fp32_gflops"] - - bw_util = df["achieved_bw_gbs"] / emp_bw - flop_util = df["achieved_gflops"] / peak_flops - m = df["m"].values - - fig, ax = plt.subplots(figsize=(10, 5), layout="constrained") - - ax.semilogx(m, bw_util * 100, "o-", color="royalblue", lw=2, ms=5, - label="Memory BW utilisation (%)") - ax.semilogx(m, flop_util * 100, "s-", color="tomato", lw=2, ms=5, - label="FP32 compute utilisation (%)") - - ax.axhline(100, color="gray", ls="--", lw=1, alpha=0.6) - ax.set_ylim(0, 115) - ax.set_xlabel("m (FMA iterations per element)") - ax.set_ylabel("Hardware Utilisation [%]") - device = info.get("device_name", "GPU") - ax.set_title(f"Hardware Utilisation vs Arithmetic Intensity — {device}") - ax.legend() - ax.grid(True, which="both", ls="--", lw=0.6, alpha=0.4) - - fig.savefig(RESULTS / "efficiency.png", bbox_inches="tight") - plt.close(fig) - print("Saved: efficiency.png") - - -# ───────────────────────────────────────────────────────────────────────────── -# Main -# ───────────────────────────────────────────────────────────────────────────── - -def main(): - info = load_device_info(RESULTS / "device_info.csv") - df_roof = load_csv("roofline.csv") - df_bw_sat = load_csv("bw_saturation.csv") - - print(f"Device: {info.get('device_name', '?')}") - print(f"Theoretical BW: {info['theoretical_bw_gbs']:.0f} GB/s") - print(f"Empirical BW: {info.get('empirical_bw_gbs', float('nan')):.0f} GB/s") - print(f"Peak FP32: {info['theoretical_fp32_gflops']:.0f} GFLOP/s") - print(f"Ridge point: {info['ridge_point_flop_per_byte']:.1f} FLOP/byte\n") - - plot_roofline(df_roof, info) - plot_bw_vs_m(df_roof, info) - plot_time_vs_m(df_roof, info) - plot_efficiency(df_roof, info) - plot_bw_saturation(df_bw_sat, info) - - -if __name__ == "__main__": - main() diff --git a/cuda-memory-compute-bound/scripts/compile.sh b/cuda-memory-compute-bound/scripts/compile.sh deleted file mode 100755 index a52b83f..0000000 --- a/cuda-memory-compute-bound/scripts/compile.sh +++ /dev/null @@ -1,41 +0,0 @@ -#!/bin/bash -# Compile the memory-bound vs compute-bound benchmark -# Usage: ./compile.sh [sm_XX] (default: sm_80 for A100) -# ./compile.sh sm_86 (RTX 3090 / A10) -# ./compile.sh sm_90 (H100) - -set -e - -THIS_DIR=$(dirname "$(realpath "$0")") -ROOT_DIR=$(dirname "${THIS_DIR}") -ARCH="${1:-sm_80}" - -echo "Building for arch=${ARCH}" - -# ── Recreate build dirs ──────────────────────────────────────────────────── -rm -rf "${ROOT_DIR}/build" -mkdir -p "${ROOT_DIR}/build/obj" "${ROOT_DIR}/build/bin" - -INCLUDES="-I${ROOT_DIR}/include" -CUDA_INCLUDES="-I/usr/include" -CUDA_LIB_DIRS="-L/usr/lib/x86_64-linux-gnu/" -CUDA_LIB="-lcudart" -FLAGS="-O3 -arch=${ARCH} --use_fast_math" - -# ── Compile objects ──────────────────────────────────────────────────────── -nvcc ${FLAGS} ${INCLUDES} ${CUDA_INCLUDES} \ - -c "${ROOT_DIR}/src/kernels.cu" \ - -o "${ROOT_DIR}/build/obj/kernels.o" - -nvcc ${FLAGS} ${INCLUDES} ${CUDA_INCLUDES} \ - -c "${ROOT_DIR}/src/main.cu" \ - -o "${ROOT_DIR}/build/obj/main.o" - -# ── Link ─────────────────────────────────────────────────────────────────── -nvcc ${FLAGS} \ - "${ROOT_DIR}/build/obj/kernels.o" \ - "${ROOT_DIR}/build/obj/main.o" \ - ${CUDA_LIB_DIRS} ${CUDA_LIB} \ - -o "${ROOT_DIR}/build/bin/main" - -echo "Binary: ${ROOT_DIR}/build/bin/main" diff --git a/cuda-memory-compute-bound/scripts/execute.sh b/cuda-memory-compute-bound/scripts/execute.sh deleted file mode 100755 index 4689b4a..0000000 --- a/cuda-memory-compute-bound/scripts/execute.sh +++ /dev/null @@ -1,50 +0,0 @@ -#!/bin/bash -# Run the benchmark and produce CSV results. -# Optional: pass --ncu to also profile with Nsight Compute. -# -# Usage: -# ./execute.sh # just benchmark -# ./execute.sh --ncu # benchmark + ncu profile of the m=1 and m=256 kernels - -set -e - -THIS_DIR=$(dirname "$(realpath "$0")") -ROOT_DIR=$(dirname "${THIS_DIR}") -BINARY="${ROOT_DIR}/build/bin/main" -RESULTS="${ROOT_DIR}/results" - -if [ ! -f "${BINARY}" ]; then - echo "Binary not found. Run scripts/compile.sh first." - exit 1 -fi - -mkdir -p "${RESULTS}" - -echo "=== Running benchmark ===" -"${BINARY}" \ - --results_dir "${RESULTS}" \ - --threads_per_block 256 \ - --N_roofline 33554432 \ - --bw_repeats 10 - -echo "" -echo "=== Results saved to ${RESULTS}/ ===" - -# ── Optional: Nsight Compute profiling ──────────────────────────────────── -if [[ "$1" == "--ncu" ]]; then - echo "" - echo "=== Nsight Compute profiling (m=1 and m=512) ===" - - # Profile the memory-bound case (m=1) - ncu --set full \ - --export "${RESULTS}/ncu_m1" \ - --force-overwrite \ - "${BINARY}" \ - --results_dir "${RESULTS}" \ - --threads_per_block 256 \ - --N_roofline 33554432 \ - --bw_repeats 1 2>/dev/null || \ - echo "ncu not available — skipping (install CUDA Toolkit for profiling)" - - echo "NCU reports: ${RESULTS}/ncu_m1.ncu-rep" -fi diff --git a/cuda-memory-compute-bound/src/kernels.cu b/cuda-memory-compute-bound/src/kernels.cu deleted file mode 100644 index 408b78a..0000000 --- a/cuda-memory-compute-bound/src/kernels.cu +++ /dev/null @@ -1,102 +0,0 @@ -#include "kernels.h" - -// ───────────────────────────────────────────────────────────────────────────── -// Vector-add + FMA benchmark kernel -// -// For small m (low AI) → memory-bandwidth-bound -// For large m (high AI) → compute-throughput-bound -// -// Key: the GPU hides FMA latency (~4-cycle pipeline) by scheduling other warps -// while a given warp waits. With enough active warps (large N), the throughput -// asymptotically approaches 2 FMAs/cycle even with a single accumulator chain -// per thread. -// ───────────────────────────────────────────────────────────────────────────── -__global__ void vector_add_fma_kernel( - const float* __restrict__ a, - const float* __restrict__ b, - float* __restrict__ c, - int N, int m) -{ - int idx = blockIdx.x * blockDim.x + threadIdx.x; - if (idx >= N) return; - - // One memory-bound load+add (1 FLOP, 8 bytes read) - float val = a[idx] + b[idx]; - - // m compute-intensive FMAs: val = val * alpha + beta (2 FLOPs each) - // alpha/beta chosen so val stays bounded for any starting value in [0,2]: - // val → converges to beta/(1-alpha) = 0.0001/(0.00001) = 10 - // #pragma unroll 1 prevents the compiler from unrolling (keeps m general) - #pragma unroll 1 - for (int j = 0; j < m; ++j) { - val = fmaf(val, 1.0f - 1e-5f, 1e-4f); - } - - // One write (4 bytes written) - c[idx] = val; -} - -// ───────────────────────────────────────────────────────────────────────────── -// Pure-copy kernel — measures empirical peak memory bandwidth -// ───────────────────────────────────────────────────────────────────────────── -__global__ void bandwidth_test_kernel( - const float* __restrict__ a, - float* __restrict__ b, - int N) -{ - int idx = blockIdx.x * blockDim.x + threadIdx.x; - if (idx < N) b[idx] = a[idx]; -} - -// ───────────────────────────────────────────────────────────────────────────── -// Tiny kernel just to force driver/context init before benchmarking -// ───────────────────────────────────────────────────────────────────────────── -__global__ void warmup_kernel(float* x, int N) -{ - int idx = blockIdx.x * blockDim.x + threadIdx.x; - if (idx < N) x[idx] = static_cast(idx); -} - -// ───────────────────────────────────────────────────────────────────────────── -// Host wrappers -// ───────────────────────────────────────────────────────────────────────────── - -void vector_add_fma( - const float* d_a, const float* d_b, float* d_c, - int N, int m, int threads_per_block, - cudaEvent_t start, cudaEvent_t stop, float& elapsed_ms) -{ - int blocks = (N + threads_per_block - 1) / threads_per_block; - - CHECK_CUDA_ERROR(cudaEventRecord(start)); - vector_add_fma_kernel<<>>(d_a, d_b, d_c, N, m); - CHECK_CUDA_ERROR(cudaEventRecord(stop)); - CHECK_LAST_CUDA_ERROR(); - CHECK_CUDA_ERROR(cudaEventSynchronize(stop)); - CHECK_CUDA_ERROR(cudaEventElapsedTime(&elapsed_ms, start, stop)); -} - -void bandwidth_test( - const float* d_a, float* d_b, - int N, int threads_per_block, - cudaEvent_t start, cudaEvent_t stop, float& elapsed_ms) -{ - int blocks = (N + threads_per_block - 1) / threads_per_block; - - CHECK_CUDA_ERROR(cudaEventRecord(start)); - bandwidth_test_kernel<<>>(d_a, d_b, N); - CHECK_CUDA_ERROR(cudaEventRecord(stop)); - CHECK_LAST_CUDA_ERROR(); - CHECK_CUDA_ERROR(cudaEventSynchronize(stop)); - CHECK_CUDA_ERROR(cudaEventElapsedTime(&elapsed_ms, start, stop)); -} - -void warmup() -{ - float* d; - CHECK_CUDA_ERROR(cudaMalloc(&d, 4096 * sizeof(float))); - warmup_kernel<<<16, 256>>>(d, 4096); - CHECK_LAST_CUDA_ERROR(); - CHECK_CUDA_ERROR(cudaDeviceSynchronize()); - CHECK_CUDA_ERROR(cudaFree(d)); -} diff --git a/cuda-memory-compute-bound/src/main.cu b/cuda-memory-compute-bound/src/main.cu deleted file mode 100644 index 11a9525..0000000 --- a/cuda-memory-compute-bound/src/main.cu +++ /dev/null @@ -1,356 +0,0 @@ -// -// main.cu — Memory-Bound vs Compute-Bound: A GPU Deep Dive with Vector Addition -// -// Benchmarks two related quantities as the arithmetic intensity (AI) rises: -// • Achieved memory bandwidth (GB/s) -// • Achieved FP32 throughput (GFLOP/s) -// -// The kernel does: -// val = a[i] + b[i] ← 1 FLOP, 8 bytes read -// for j in 0..m: val = fma(val,α,β) ← 2m FLOPs, no extra memory -// c[i] = val ← 0 FLOPs, 4 bytes written -// -// AI = (1 + 2m) / 12 [FLOPs / byte] -// -// Sweeps: -// 1. ROOFLINE — fixed large N, m ∈ {1,2,4,…,16384} -// Produces the roofline-model scatter plot. -// 2. BW_SAT — fixed m=1, N ∈ {2^10 … 2^27} -// Shows how bandwidth saturates as the working set grows. -// -// Output: two CSV files in / -// roofline.csv -// bw_saturation.csv -// -// Usage: -// ./main --results_dir ../results [--threads_per_block 256] -// [--N_roofline 33554432] [--bw_repeats 10] - -#include "kernels.h" -#include "utils.h" - -#include -#include -#include -#include -#include -#include -#include -#include -#include - -// ───────────────────────────────────────────────────────────────────────────── -// Theoretical peak values from device properties -// ───────────────────────────────────────────────────────────────────────────── - -// CUDA cores per SM keyed on (major, minor) compute capability -static int cores_per_sm(int major, int minor) -{ - struct Entry { int major, minor, cores; }; - static constexpr Entry table[] = { - {3, 0, 192}, {3, 2, 192}, {3, 5, 192}, {3, 7, 192}, - {5, 0, 128}, {5, 2, 128}, - {6, 0, 64}, {6, 1, 128}, {6, 2, 128}, - {7, 0, 64}, {7, 2, 64}, {7, 5, 64}, - {8, 0, 64}, {8, 6, 128}, {8, 7, 128}, {8, 9, 128}, - {9, 0, 128}, - }; - for (auto& e : table) - if (e.major == major && e.minor == minor) return e.cores; - return 128; // safe default -} - -// Theoretical peak HBM / GDDR bandwidth in GB/s -static double theoretical_bw_gbs(const cudaDeviceProp& p) -{ - // memoryClockRate is in kHz; memoryBusWidth is in bits - // Bandwidth = clock_Hz × (bus_width_bits / 8) × 2 (DDR) / 1e9 - return static_cast(p.memoryClockRate) * 1e3 // → Hz - * (p.memoryBusWidth / 8.0) // → bytes/cycle - * 2.0 // DDR - / 1e9; -} - -// Theoretical peak FP32 in GFLOP/s -static double theoretical_fp32_gflops(const cudaDeviceProp& p) -{ - // Each CUDA core can issue 1 FMA (= 2 FLOPs) per clock cycle - return 2.0 - * p.multiProcessorCount - * cores_per_sm(p.major, p.minor) - * static_cast(p.clockRate) * 1e3 // → Hz - / 1e9; -} - -// ───────────────────────────────────────────────────────────────────────────── -// Helpers -// ───────────────────────────────────────────────────────────────────────────── - -static float time_kernel_ms( - const float* d_a, const float* d_b, float* d_c, - int N, int m, int tpb, int repeats, - cudaEvent_t ev0, cudaEvent_t ev1) -{ - float best = 1e30f; - for (int r = 0; r < repeats; ++r) { - float t; - vector_add_fma(d_a, d_b, d_c, N, m, tpb, ev0, ev1, t); - best = std::min(best, t); - } - return best; -} - -static float time_bw_kernel_ms( - const float* d_a, float* d_c, - int N, int tpb, int repeats, - cudaEvent_t ev0, cudaEvent_t ev1) -{ - float best = 1e30f; - for (int r = 0; r < repeats; ++r) { - float t; - bandwidth_test(d_a, d_c, N, tpb, ev0, ev1, t); - best = std::min(best, t); - } - return best; -} - -// ───────────────────────────────────────────────────────────────────────────── -// Arg parsing -// ───────────────────────────────────────────────────────────────────────────── - -struct Config { - std::string results_dir = "../results"; - int threads_per_block = 256; - long N_roofline = 1L << 25; // 33 554 432 floats = 128 MB / array - int bw_repeats = 10; -}; - -static Config parse_args(int argc, char** argv) -{ - Config cfg; - for (int i = 1; i < argc; ++i) { - std::string a = argv[i]; - if (a == "--results_dir" && i + 1 < argc) cfg.results_dir = argv[++i]; - if (a == "--threads_per_block"&& i + 1 < argc) cfg.threads_per_block = std::stoi(argv[++i]); - if (a == "--N_roofline" && i + 1 < argc) cfg.N_roofline = std::stol(argv[++i]); - if (a == "--bw_repeats" && i + 1 < argc) cfg.bw_repeats = std::stoi(argv[++i]); - if (a == "--help") { - std::cout << "Usage: ./main [--results_dir DIR] [--threads_per_block INT]\n" - << " [--N_roofline INT] [--bw_repeats INT]\n"; - std::exit(0); - } - } - return cfg; -} - -// ───────────────────────────────────────────────────────────────────────────── -// main -// ───────────────────────────────────────────────────────────────────────────── - -int main(int argc, char** argv) -{ - Config cfg = parse_args(argc, argv); - - // ── Device info ───────────────────────────────────────────────────────── - cudaDeviceProp prop; - CHECK_CUDA_ERROR(cudaGetDeviceProperties(&prop, 0)); - - double th_bw = theoretical_bw_gbs(prop); - double th_flops = theoretical_fp32_gflops(prop); - double ridge_pt = th_flops / th_bw; // FLOPs/byte at the ridge - - std::cout << "Device : " << prop.name << "\n" - << "Compute cap. : " << prop.major << "." << prop.minor << "\n" - << "SMs : " << prop.multiProcessorCount << "\n" - << "CUDA cores/SM : " << cores_per_sm(prop.major, prop.minor) << "\n" - << "GPU clock : " << prop.clockRate / 1e6 << " GHz\n" - << "Mem clock : " << prop.memoryClockRate / 1e6 << " GHz\n" - << "Bus width : " << prop.memoryBusWidth << " bits\n" - << "Theoretical BW : " << th_bw << " GB/s\n" - << "Theoretical FLOPS: " << th_flops << " GFLOP/s\n" - << "Ridge point : " << ridge_pt << " FLOP/byte\n\n"; - - // ── Make results directory ─────────────────────────────────────────────── - std::filesystem::create_directories(cfg.results_dir); - - // ── Shared CUDA events ─────────────────────────────────────────────────── - cudaEvent_t ev0, ev1; - CHECK_CUDA_ERROR(cudaEventCreate(&ev0)); - CHECK_CUDA_ERROR(cudaEventCreate(&ev1)); - - // ── Warm up ────────────────────────────────────────────────────────────── - warmup(); - warmup(); - - // ── Measure empirical peak bandwidth (pure-copy kernel) ────────────────── - { - const long N_bw = cfg.N_roofline; - float *d_a, *d_b; - CHECK_CUDA_ERROR(cudaMalloc(&d_a, N_bw * sizeof(float))); - CHECK_CUDA_ERROR(cudaMalloc(&d_b, N_bw * sizeof(float))); - CHECK_CUDA_ERROR(cudaMemset(d_a, 1, N_bw * sizeof(float))); - - float t = time_bw_kernel_ms(d_a, d_b, static_cast(N_bw), - cfg.threads_per_block, cfg.bw_repeats, ev0, ev1); - - // copy reads 1 array + writes 1 array = 2 × N × sizeof(float) bytes - double bytes = 2.0 * N_bw * sizeof(float); - double emp_bw = bytes / (t * 1e-3) / 1e9; - - std::cout << "Empirical peak BW: " << emp_bw << " GB/s (vs " - << th_bw << " GB/s theoretical)\n\n"; - - // Save device info + measured peak to a small metadata file - std::ofstream meta(cfg.results_dir + "/device_info.csv"); - meta << "key,value\n" - << "device_name," << prop.name << "\n" - << "compute_capability," << prop.major << "." << prop.minor << "\n" - << "num_sms," << prop.multiProcessorCount << "\n" - << "cores_per_sm," << cores_per_sm(prop.major, prop.minor) << "\n" - << "gpu_clock_ghz," << prop.clockRate / 1e6 << "\n" - << "mem_clock_ghz," << prop.memoryClockRate / 1e6 << "\n" - << "mem_bus_width_bits," << prop.memoryBusWidth << "\n" - << "theoretical_bw_gbs," << th_bw << "\n" - << "theoretical_fp32_gflops," << th_flops << "\n" - << "ridge_point_flop_per_byte," << ridge_pt << "\n" - << "empirical_bw_gbs," << emp_bw << "\n"; - - CHECK_CUDA_ERROR(cudaFree(d_a)); - CHECK_CUDA_ERROR(cudaFree(d_b)); - } - - // ═════════════════════════════════════════════════════════════════════════ - // SWEEP 1: ROOFLINE — fixed N, sweep m - // - // For each m we record: - // • kernel_time_ms - // • achieved_bandwidth_gbs = 12·N / (time_s · 1e9) - // • achieved_gflops = (1+2m)·N / (time_s · 1e9) - // • arithmetic_intensity = (1+2m) / 12 - // ═════════════════════════════════════════════════════════════════════════ - { - const long N = cfg.N_roofline; - const double bytes_per_elem = 3.0 * sizeof(float); // 2 reads + 1 write - - // m values span several orders of magnitude to cross the ridge point - std::vector m_values; - for (int m = 1; m <= 16384; m *= 2) m_values.push_back(m); - // add fine-grained points near the expected ridge - for (int m : {3, 6, 12, 24, 48, 96, 192, 384, 768, 1536, 3072, 6144, 12288}) - m_values.push_back(m); - std::sort(m_values.begin(), m_values.end()); - m_values.erase(std::unique(m_values.begin(), m_values.end()), m_values.end()); - - float *d_a, *d_b, *d_c; - CHECK_CUDA_ERROR(cudaMalloc(&d_a, N * sizeof(float))); - CHECK_CUDA_ERROR(cudaMalloc(&d_b, N * sizeof(float))); - CHECK_CUDA_ERROR(cudaMalloc(&d_c, N * sizeof(float))); - - // Fill with random data (on host, then copy) - std::vector h(N); - std::mt19937 rng(42); - std::uniform_real_distribution dist(0.0f, 1.0f); - std::generate(h.begin(), h.end(), [&]{ return dist(rng); }); - CHECK_CUDA_ERROR(cudaMemcpy(d_a, h.data(), N * sizeof(float), cudaMemcpyHostToDevice)); - CHECK_CUDA_ERROR(cudaMemcpy(d_b, h.data(), N * sizeof(float), cudaMemcpyHostToDevice)); - - std::ofstream csv(cfg.results_dir + "/roofline.csv"); - csv << "m,N,kernel_time_ms,arithmetic_intensity,achieved_bw_gbs,achieved_gflops," - "theoretical_bw_gbs,theoretical_fp32_gflops,ridge_point\n"; - - std::cout << "=== ROOFLINE SWEEP (N=" << N << ") ===\n"; - std::cout << " m AI(FLOPs/B) BW(GB/s) GFLOP/s\n"; - - for (int m : m_values) { - float t_ms = time_kernel_ms(d_a, d_b, d_c, static_cast(N), - m, cfg.threads_per_block, - cfg.bw_repeats, ev0, ev1); - - double t_s = t_ms * 1e-3; - double flops = static_cast(N) * (1.0 + 2.0 * m); - double bytes = static_cast(N) * bytes_per_elem; - double ai = flops / bytes; - double bw_gbs = bytes / t_s / 1e9; - double gflops = flops / t_s / 1e9; - - csv << m << "," << N << "," << t_ms << "," - << ai << "," << bw_gbs << "," << gflops << "," - << th_bw << "," << th_flops << "," << ridge_pt << "\n"; - - std::cout << " m=" << m - << " AI=" << ai - << " BW=" << bw_gbs << " GB/s" - << " GFLOP/s=" << gflops << "\n"; - } - - csv.close(); - CHECK_CUDA_ERROR(cudaFree(d_a)); - CHECK_CUDA_ERROR(cudaFree(d_b)); - CHECK_CUDA_ERROR(cudaFree(d_c)); - std::cout << "\n"; - } - - // ═════════════════════════════════════════════════════════════════════════ - // SWEEP 2: BANDWIDTH SATURATION — m=1, sweep N - // - // Shows that small vectors can't saturate HBM bandwidth (the GPU doesn't - // have enough warps in flight to fill the memory controllers), while large - // vectors approach the empirical peak. - // ═════════════════════════════════════════════════════════════════════════ - { - const int m = 1; // minimal compute — purely memory-bound - const double bytes_per_elem = 3.0 * sizeof(float); - - // N from 1024 to N_roofline in powers of 2 - std::vector n_values; - for (long n = 1024; n <= cfg.N_roofline; n *= 2) - n_values.push_back(n); - - std::ofstream csv(cfg.results_dir + "/bw_saturation.csv"); - csv << "m,N,sizeMB,kernel_time_ms,arithmetic_intensity," - "achieved_bw_gbs,achieved_gflops," - "theoretical_bw_gbs,theoretical_fp32_gflops\n"; - - std::cout << "=== BANDWIDTH SATURATION SWEEP (m=" << m << ") ===\n"; - std::cout << " N sizeMB BW(GB/s)\n"; - - for (long N : n_values) { - float *d_a, *d_b, *d_c; - CHECK_CUDA_ERROR(cudaMalloc(&d_a, N * sizeof(float))); - CHECK_CUDA_ERROR(cudaMalloc(&d_b, N * sizeof(float))); - CHECK_CUDA_ERROR(cudaMalloc(&d_c, N * sizeof(float))); - CHECK_CUDA_ERROR(cudaMemset(d_a, 1, N * sizeof(float))); - CHECK_CUDA_ERROR(cudaMemset(d_b, 2, N * sizeof(float))); - - float t_ms = time_kernel_ms(d_a, d_b, d_c, static_cast(N), - m, cfg.threads_per_block, - cfg.bw_repeats, ev0, ev1); - - double t_s = t_ms * 1e-3; - double flops = static_cast(N) * (1.0 + 2.0 * m); - double bytes = static_cast(N) * bytes_per_elem; - double ai = flops / bytes; - double bw_gbs = bytes / t_s / 1e9; - double gflops = flops / t_s / 1e9; - double sizeMB = static_cast(N) * sizeof(float) / (1024.0 * 1024.0); - - csv << m << "," << N << "," << sizeMB << "," << t_ms << "," - << ai << "," << bw_gbs << "," << gflops << "," - << th_bw << "," << th_flops << "\n"; - - std::cout << " N=" << N << " " << sizeMB << " MB BW=" << bw_gbs << " GB/s\n"; - - CHECK_CUDA_ERROR(cudaFree(d_a)); - CHECK_CUDA_ERROR(cudaFree(d_b)); - CHECK_CUDA_ERROR(cudaFree(d_c)); - } - - csv.close(); - } - - CHECK_CUDA_ERROR(cudaEventDestroy(ev0)); - CHECK_CUDA_ERROR(cudaEventDestroy(ev1)); - - std::cout << "\nResults written to: " << cfg.results_dir << "/\n"; - return 0; -} diff --git a/cuda-mma/Makefile b/cuda-mma/Makefile new file mode 100644 index 0000000..f94d09d --- /dev/null +++ b/cuda-mma/Makefile @@ -0,0 +1,45 @@ +# --------------------------------------------------------------------------- +# 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 + +NVCC_FLAGS := $(CXX_FLAGS) -arch=$(ARCH) \ + -rdc=true \ + --generate-line-info \ + -Xcompiler -Wall \ + -I. -Iinclude + +TARGET := cuda_mma +SRCS := main.cu src/kernels.cu src/utils.cu +HDRS := include/cuda_check.cuh include/kernels.cuh \ + include/kernels.cuh include/timer.h include/utils.h + +.PHONY: all clean run profile + +all: $(TARGET) + +$(TARGET): $(SRCS) $(HDRS) + $(NVCC) $(NVCC_FLAGS) -o $@ $(SRCS) + +run: $(TARGET) + ./$(TARGET) + +# 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 -f $(TARGET) diff --git a/cuda-mma/include/cuda_check.cuh b/cuda-mma/include/cuda_check.cuh new file mode 100644 index 0000000..80f1368 --- /dev/null +++ b/cuda-mma/include/cuda_check.cuh @@ -0,0 +1,40 @@ +// cuda_check.cuh +#pragma once + +#include +#include +#include + +// ── 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) diff --git a/cuda-mma/include/kernels.cuh b/cuda-mma/include/kernels.cuh new file mode 100644 index 0000000..2fdf4a6 --- /dev/null +++ b/cuda-mma/include/kernels.cuh @@ -0,0 +1,38 @@ +#pragma once + +#include + + +__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 +__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 +__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) \ No newline at end of file diff --git a/cuda-mma/include/timer.h b/cuda-mma/include/timer.h new file mode 100644 index 0000000..89e3a62 --- /dev/null +++ b/cuda-mma/include/timer.h @@ -0,0 +1,29 @@ +#pragma once + +#include + +// --------------------------------------------------------------------------- +// 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; + } +}; diff --git a/cuda-mma/include/utils.h b/cuda-mma/include/utils.h new file mode 100644 index 0000000..fc4d6a4 --- /dev/null +++ b/cuda-mma/include/utils.h @@ -0,0 +1,17 @@ +#pragma once + +#include + +// 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); diff --git a/cuda-mma/main.cu b/cuda-mma/main.cu new file mode 100644 index 0000000..b681e53 --- /dev/null +++ b/cuda-mma/main.cu @@ -0,0 +1,242 @@ +#include +#include +#include +#include + +#include "cuda_check.cuh" +#include "kernels.cuh" +#include "utils.h" +#include "timer.h" + +// --------------------------------------------------------------------------- +// Benchmark config +// --------------------------------------------------------------------------- +#define WARMUP 3 +#define ITERS 5 + +// --------------------------------------------------------------------------- +// Benchmark helper +// Allocates device memory for a square S×S problem, runs the naive kernel +// WARMUP+ITERS times, and returns the average elapsed milliseconds. +// Host data is filled once at startup; each timed run resets d_C to zero so +// that C += A*B doesn't accumulate across iterations. +// --------------------------------------------------------------------------- +static float benchmark_naive(int S, + const float* h_A, + const float* h_B) +{ + size_t szA = (size_t)S * S * sizeof(float); + size_t szB = (size_t)S * S * sizeof(float); + size_t szC = (size_t)S * S * sizeof(float); + + float *d_A, *d_B, *d_C; + CUDA_CHECK(cudaMalloc(&d_A, szA)); + CUDA_CHECK(cudaMalloc(&d_B, szB)); + CUDA_CHECK(cudaMalloc(&d_C, szC)); + + CUDA_CHECK(cudaMemcpy(d_A, h_A, szA, cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemcpy(d_B, h_B, szB, cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemset(d_C, 0, szC)); + + dim3 gridDim(CEIL_DIV(S, 32), CEIL_DIV(S, 32)); + dim3 blockDim(32, 32, 1); + + // Warmup — not measured. + for (int i = 0; i < WARMUP; ++i) + sgemm_naive<<>>(S, S, S, 1.0f, d_A, d_B, 0.0f, d_C); + CUDA_CHECK(cudaDeviceSynchronize()); + + // Timed runs. + GpuTimer timer; + float total_ms = 0.0f; + for (int i = 0; i < ITERS; ++i) { + CUDA_CHECK(cudaMemset(d_C, 0, szC)); + timer.start(); + sgemm_naive<<>>(S, S, S, 1.0f, d_A, d_B, 0.0f, d_C); + total_ms += timer.stop(); + } + + CUDA_CHECK(cudaFree(d_A)); + CUDA_CHECK(cudaFree(d_B)); + CUDA_CHECK(cudaFree(d_C)); + + return total_ms / ITERS; +} + +static float benchmark_tiled(int S, + const float* h_A, + const float* h_B) +{ + size_t szA = (size_t)S * S * sizeof(float); + size_t szB = (size_t)S * S * sizeof(float); + size_t szC = (size_t)S * S * sizeof(float); + + float *d_A, *d_B, *d_C; + CUDA_CHECK(cudaMalloc(&d_A, szA)); + CUDA_CHECK(cudaMalloc(&d_B, szB)); + CUDA_CHECK(cudaMalloc(&d_C, szC)); + + CUDA_CHECK(cudaMemcpy(d_A, h_A, szA, cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemcpy(d_B, h_B, szB, cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemset(d_C, 0, szC)); + + dim3 blockDim(32, 32, 1); + dim3 gridDim(CEIL_DIV(S, blockDim.x / 32), CEIL_DIV(S, blockDim.y / 32)); + + // Warmup — not measured. + for (int i = 0; i < WARMUP; ++i) + sgemm_tiled<16><<>>(S, S, S, 1.0f, d_A, d_B, 0.0f, d_C); + CUDA_CHECK(cudaDeviceSynchronize()); + + // Timed runs. + GpuTimer timer; + float total_ms = 0.0f; + for (int i = 0; i < ITERS; ++i) { + CUDA_CHECK(cudaMemset(d_C, 0, szC)); + timer.start(); + sgemm_tiled<16><<>>(S, S, S, 1.0f, d_A, d_B, 0.0f, d_C); + total_ms += timer.stop(); + } + + CUDA_CHECK(cudaFree(d_A)); + CUDA_CHECK(cudaFree(d_B)); + CUDA_CHECK(cudaFree(d_C)); + + return total_ms / ITERS; +} + + +void benchmark_coalesced(int S, + const float* h_A, + const float* h_B) +{ + size_t szA = (size_t)S * S * sizeof(float); + size_t szB = (size_t)S * S * sizeof(float); + size_t szC = (size_t)S * S * sizeof(float); + + float *d_A, *d_B, *d_C; + CUDA_CHECK(cudaMalloc(&d_A, szA)); + CUDA_CHECK(cudaMalloc(&d_B, szB)); + CUDA_CHECK(cudaMalloc(&d_C, szC)); + + CUDA_CHECK(cudaMemcpy(d_A, h_A, szA, cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemcpy(d_B, h_B, szB, cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemset(d_C, 0, szC)); + + dim3 block(32 * 32); + dim3 grid((S + block.x / 16 - 1) / (block.x / 16), + (S + block.x / 16 - 1) / (block.x / 16)); + + // Warmup — not measured. + for (int i = 0; i < WARMUP; ++i) + sgemm_coalesced<16><<>>(S, S, S, 1.0f, d_A, d_B, 0.0f, d_C); + CUDA_CHECK(cudaDeviceSynchronize()); + + // Timed runs. + GpuTimer timer; + float total_ms = 0.0f; + for (int i = 0; i < ITERS; ++i) { + CUDA_CHECK(cudaMemset(d_C, 0, szC)); + timer.start(); + total_ms += timer.stop(); + } + + CUDA_CHECK(cudaFree(d_A)); + CUDA_CHECK(cudaFree(d_B)); + CUDA_CHECK(cudaFree(d_C)); + + printf("Coalesced kernel: %9.3f ms\n", total_ms / ITERS); +} + +// --------------------------------------------------------------------------- +// Main +// --------------------------------------------------------------------------- +int main() { + srand(42); + + // // ── Correctness check (small fixed size) ──────────────────────────────── + // // The CPU triple-loop is O(S^3), so we only run it for a small S. + // const int S_verify = 256; + // { + // float* h_A = new float[S_verify * S_verify]; + // float* h_B = new float[S_verify * S_verify]; + // float* h_C = new float[S_verify * S_verify](); // zero-init + // float* h_C_ref = new float[S_verify * S_verify](); + + // fill_random(h_A, S_verify * S_verify); + // fill_random(h_B, S_verify * S_verify); + // matmul_cpu(h_A, h_B, h_C_ref, S_verify, S_verify, S_verify); + + // float *d_A, *d_B, *d_C; + // size_t sz = (size_t)S_verify * S_verify * sizeof(float); + // CUDA_CHECK(cudaMalloc(&d_A, sz)); + // CUDA_CHECK(cudaMalloc(&d_B, sz)); + // CUDA_CHECK(cudaMalloc(&d_C, sz)); + // CUDA_CHECK(cudaMemcpy(d_A, h_A, sz, cudaMemcpyHostToDevice)); + // CUDA_CHECK(cudaMemcpy(d_B, h_B, sz, cudaMemcpyHostToDevice)); + // CUDA_CHECK(cudaMemset(d_C, 0, sz)); + + // constexpr int BS = 32; + // dim3 block(BS * BS); // BLOCKSIZE² threads: threadIdx.x / BS = row, threadIdx.x % BS = col + // dim3 grid((S_verify + BS - 1) / BS, (S_verify + BS - 1) / BS); + // sgemm_coalesced<<>>(S_verify, S_verify, S_verify, 1.0f, d_A, d_B, 0.0f, d_C); + // CHECK_LAST_ERROR(); + // CHECK_SYNC(); + // CUDA_CHECK(cudaMemcpy(h_C, d_C, sz, cudaMemcpyDeviceToHost)); + + // printf("Correctness check (S=%d): ", S_verify); + // if (verify(h_C_ref, h_C, S_verify * S_verify)) + // printf("PASSED\n\n"); + // else + // printf("FAILED\n\n"); + + // CUDA_CHECK(cudaFree(d_A)); + // CUDA_CHECK(cudaFree(d_B)); + // CUDA_CHECK(cudaFree(d_C)); + // delete[] h_A; delete[] h_B; delete[] h_C; delete[] h_C_ref; + // } + + + const int sizes[] = {128, 256, 512, 1024, 2048, 4096, 8192, 16384, 32768}; + const int N_SIZES = sizeof(sizes) / sizeof(sizes[0]); + + // Allocate host matrices for the largest size once. + 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); + + printf("Naive kernel roofline sweep (%d iters, %d warmup)\n", ITERS, WARMUP); + printf("%-6s %9s %9s %11s %9s\n", + "Size", "Time(ms)", "GFLOP/s", "BW(GB/s)", "AI(F/B)"); + printf("------ --------- --------- ----------- ---------\n"); + + for (int i = 0; i < N_SIZES; ++i) { + int S = sizes[i]; + + float avg_ms = benchmark_naive(S, h_A, h_B); + float avg_ms_tiled = benchmark_tiled(S, h_A, h_B); + + // Each term is promoted to double before multiplying to prevent int32 + // overflow at large S (see note above). + double flops = 2.0 * S * S * S; + double bytes = ((double)S * S // read A + + (double)S * S // read B + + 2.0 * S * S) // read + write C + * sizeof(float); + 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); + + // printf(" Tiled kernel: %9.3f ms\n", avg_ms_tiled); + } + + delete[] h_A; + delete[] h_B; + + return 0; +} diff --git a/cuda-mma/src/kernels.cu b/cuda-mma/src/kernels.cu new file mode 100644 index 0000000..bb530c1 --- /dev/null +++ b/cuda-mma/src/kernels.cu @@ -0,0 +1,108 @@ +#include + +#include "kernels.cuh" + +// Single Precision Matrix Multiplication Kernels SGEMM: C = alpha * A * B + beta * C + + +// Naive implementation: 1 thread per output element, no shared memory, non-coalesced accesses. +__global__ void sgemm_naive( + 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) + { + + int row = blockIdx.y * blockDim.y + threadIdx.y; // 0 .. M-1 + int col = blockIdx.x * blockDim.x + threadIdx.x; // 0 .. N-1 + + if (row >= M || col >= N) return; + + float acc = 0.0f; + for (size_t i = 0; i < K; ++i) + acc += A[row * K + i] * B[i * N + col]; + + C[row * N + col] = alpha * acc + beta * C[row * N + col]; +} + + +// Coalesced access version: 1 thread per output element, no shared memory, but coalesced accesses to A and B. +template +__global__ void sgemm_coalesced( + size_t M, + size_t N, + size_t K, + float alpha, + const float *A, + const float *B, + float beta, + float *C) { + const int cRow = blockIdx.x * BLOCKSIZE + (threadIdx.x / BLOCKSIZE); + const int cCol = blockIdx.y * BLOCKSIZE + (threadIdx.x % BLOCKSIZE); + + if (cRow < M && cCol < N) { + float tmp = 0.0; + for (size_t i = 0; i < K; ++i) { + tmp += A[cRow * K + i] * B[i * N + cCol]; + } + C[cRow * N + cCol] = alpha * tmp + beta * C[cRow * N + cCol]; + } +} + +// Tiled version: 1 thread per output element, shared memory for tiles of A and B. +template +__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) +{ + __shared__ float sA[TILE_SIZE][TILE_SIZE]; + __shared__ float sB[TILE_SIZE][TILE_SIZE]; + + int row = blockIdx.y * TILE_SIZE + threadIdx.y; + int col = blockIdx.x * TILE_SIZE + threadIdx.x; + + float acc = 0.0f; + + // Iterate over tiles along the K dimension. + for (size_t t = 0; t < (K + TILE_SIZE - 1) / TILE_SIZE; ++t) { + size_t aCol = t * TILE_SIZE + threadIdx.x; // column of A this thread loads + size_t bRow = t * TILE_SIZE + threadIdx.y; // row of B this thread loads + + // Boundary-safe loads: pad with 0 for out-of-bounds tiles. + sA[threadIdx.y][threadIdx.x] = (row < M && aCol < K) ? A[row * K + aCol] : 0.0f; + sB[threadIdx.y][threadIdx.x] = (bRow < K && col < N) ? B[bRow * N + col] : 0.0f; + __syncthreads(); + + // Accumulate partial dot product from shared memory. + #pragma unroll + for (int i = 0; i < TILE_SIZE; ++i) + acc += sA[threadIdx.y][i] * sB[i][threadIdx.x]; + __syncthreads(); + } + + if (row < M && col < N) + C[row * N + col] = alpha * acc + beta * C[row * N + col]; +} + +template __global__ void sgemm_tiled<16>(size_t M, size_t N, size_t K, float alpha, + const float *A, const float *B, + float beta, float *C); +template __global__ void sgemm_tiled<32>(size_t M, size_t N, size_t K, float alpha, + const float *A, const float *B, + float beta, float *C); +template __global__ void sgemm_coalesced<16>(size_t M, size_t N, size_t K, float alpha, + const float *A, const float *B, + float beta, float *C); +template __global__ void sgemm_coalesced<32>(size_t M, size_t N, size_t K, float alpha, + const float *A, const float *B, + float beta, float *C); diff --git a/cuda-mma/src/utils.cu b/cuda-mma/src/utils.cu new file mode 100644 index 0000000..ede93a7 --- /dev/null +++ b/cuda-mma/src/utils.cu @@ -0,0 +1,38 @@ +#include +#include +#include + +#include "utils.h" + +void fill_random(float* p, int n) { + for (int i = 0; i < n; ++i) + p[i] = (float)rand() / RAND_MAX * 2.0f - 1.0f; +} + +void matmul_cpu(const float* A, const float* B, float* C, int m, int n, int k) { + for (int row = 0; row < m; ++row) + for (int col = 0; col < n; ++col) { + float acc = 0.0f; + for (int i = 0; i < k; ++i) + acc += A[row * k + i] * B[i * n + col]; + C[row * n + col] += acc; + } +} + +bool verify(const float* ref, const float* gpu, int total_elements, + float atol, float rtol) { + int mismatches = 0; + for (int i = 0; i < total_elements; ++i) { + float diff = fabsf(ref[i] - gpu[i]); + float tol = atol + rtol * fabsf(ref[i]); + if (diff > tol) { + if (mismatches < 5) + fprintf(stderr, " Mismatch at [%d]: ref=%.6f gpu=%.6f diff=%.2e\n", + i, ref[i], gpu[i], diff); + ++mismatches; + } + } + if (mismatches > 0) + fprintf(stderr, " Total mismatches: %d / %d\n", mismatches, total_elements); + return mismatches == 0; +} From 35ef9e6a3e0c2529a93c3aace90d5a72058671f6 Mon Sep 17 00:00:00 2001 From: Sebas <12949040+SebastiaAgramunt@users.noreply.github.com> Date: Sat, 18 Apr 2026 07:17:11 +0000 Subject: [PATCH 03/21] working --- cuda-mma/.gitignore | 1 + cuda-mma/Makefile | 28 +++- cuda-mma/include/benchmarks.h | 5 + cuda-mma/main.cu | 233 ++++-------------------------- cuda-mma/src/benchmarks.cu | 125 ++++++++++++++++ cuda-mma/src/kernels.cu | 113 +++++++++++++-- cuda-mma/test/test_correctness.cu | 50 +++++++ 7 files changed, 324 insertions(+), 231 deletions(-) create mode 100644 cuda-mma/.gitignore create mode 100644 cuda-mma/include/benchmarks.h create mode 100644 cuda-mma/src/benchmarks.cu create mode 100644 cuda-mma/test/test_correctness.cu diff --git a/cuda-mma/.gitignore b/cuda-mma/.gitignore new file mode 100644 index 0000000..567609b --- /dev/null +++ b/cuda-mma/.gitignore @@ -0,0 +1 @@ +build/ diff --git a/cuda-mma/Makefile b/cuda-mma/Makefile index f94d09d..70e8cb7 100644 --- a/cuda-mma/Makefile +++ b/cuda-mma/Makefile @@ -17,21 +17,33 @@ NVCC_FLAGS := $(CXX_FLAGS) -arch=$(ARCH) \ -Xcompiler -Wall \ -I. -Iinclude -TARGET := cuda_mma -SRCS := main.cu src/kernels.cu src/utils.cu -HDRS := include/cuda_check.cuh include/kernels.cuh \ - include/kernels.cuh include/timer.h include/utils.h +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 +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 profile +.PHONY: all clean run test profile -all: $(TARGET) +all: $(TARGET) $(TEST) -$(TARGET): $(SRCS) $(HDRS) +$(BIN_DIR): + mkdir -p $(BIN_DIR) + +$(TARGET): $(SRCS) $(HDRS) | $(BIN_DIR) $(NVCC) $(NVCC_FLAGS) -o $@ $(SRCS) +$(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 \ @@ -42,4 +54,4 @@ profile: $(TARGET) ./$(TARGET) clean: - rm -f $(TARGET) + rm -rf build diff --git a/cuda-mma/include/benchmarks.h b/cuda-mma/include/benchmarks.h new file mode 100644 index 0000000..6aa05b8 --- /dev/null +++ b/cuda-mma/include/benchmarks.h @@ -0,0 +1,5 @@ +#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); diff --git a/cuda-mma/main.cu b/cuda-mma/main.cu index b681e53..ba9e25d 100644 --- a/cuda-mma/main.cu +++ b/cuda-mma/main.cu @@ -1,239 +1,58 @@ #include #include -#include -#include -#include "cuda_check.cuh" -#include "kernels.cuh" #include "utils.h" -#include "timer.h" +#include "benchmarks.h" -// --------------------------------------------------------------------------- -// Benchmark config -// --------------------------------------------------------------------------- #define WARMUP 3 #define ITERS 5 -// --------------------------------------------------------------------------- -// Benchmark helper -// Allocates device memory for a square S×S problem, runs the naive kernel -// WARMUP+ITERS times, and returns the average elapsed milliseconds. -// Host data is filled once at startup; each timed run resets d_C to zero so -// that C += A*B doesn't accumulate across iterations. -// --------------------------------------------------------------------------- -static float benchmark_naive(int S, - const float* h_A, - const float* h_B) -{ - size_t szA = (size_t)S * S * sizeof(float); - size_t szB = (size_t)S * S * sizeof(float); - size_t szC = (size_t)S * S * sizeof(float); - - float *d_A, *d_B, *d_C; - CUDA_CHECK(cudaMalloc(&d_A, szA)); - CUDA_CHECK(cudaMalloc(&d_B, szB)); - CUDA_CHECK(cudaMalloc(&d_C, szC)); - - CUDA_CHECK(cudaMemcpy(d_A, h_A, szA, cudaMemcpyHostToDevice)); - CUDA_CHECK(cudaMemcpy(d_B, h_B, szB, cudaMemcpyHostToDevice)); - CUDA_CHECK(cudaMemset(d_C, 0, szC)); - - dim3 gridDim(CEIL_DIV(S, 32), CEIL_DIV(S, 32)); - dim3 blockDim(32, 32, 1); - - // Warmup — not measured. - for (int i = 0; i < WARMUP; ++i) - sgemm_naive<<>>(S, S, S, 1.0f, d_A, d_B, 0.0f, d_C); - CUDA_CHECK(cudaDeviceSynchronize()); - - // Timed runs. - GpuTimer timer; - float total_ms = 0.0f; - for (int i = 0; i < ITERS; ++i) { - CUDA_CHECK(cudaMemset(d_C, 0, szC)); - timer.start(); - sgemm_naive<<>>(S, S, S, 1.0f, d_A, d_B, 0.0f, d_C); - total_ms += timer.stop(); - } - - CUDA_CHECK(cudaFree(d_A)); - CUDA_CHECK(cudaFree(d_B)); - CUDA_CHECK(cudaFree(d_C)); - - return total_ms / ITERS; -} - -static float benchmark_tiled(int S, - const float* h_A, - const float* h_B) -{ - size_t szA = (size_t)S * S * sizeof(float); - size_t szB = (size_t)S * S * sizeof(float); - size_t szC = (size_t)S * S * sizeof(float); - - float *d_A, *d_B, *d_C; - CUDA_CHECK(cudaMalloc(&d_A, szA)); - CUDA_CHECK(cudaMalloc(&d_B, szB)); - CUDA_CHECK(cudaMalloc(&d_C, szC)); - - CUDA_CHECK(cudaMemcpy(d_A, h_A, szA, cudaMemcpyHostToDevice)); - CUDA_CHECK(cudaMemcpy(d_B, h_B, szB, cudaMemcpyHostToDevice)); - CUDA_CHECK(cudaMemset(d_C, 0, szC)); - - dim3 blockDim(32, 32, 1); - dim3 gridDim(CEIL_DIV(S, blockDim.x / 32), CEIL_DIV(S, blockDim.y / 32)); +using BenchmarkFn = float (*)(int, const float*, const float*); - // Warmup — not measured. - for (int i = 0; i < WARMUP; ++i) - sgemm_tiled<16><<>>(S, S, S, 1.0f, d_A, d_B, 0.0f, d_C); - CUDA_CHECK(cudaDeviceSynchronize()); - - // Timed runs. - GpuTimer timer; - float total_ms = 0.0f; - for (int i = 0; i < ITERS; ++i) { - CUDA_CHECK(cudaMemset(d_C, 0, szC)); - timer.start(); - sgemm_tiled<16><<>>(S, S, S, 1.0f, d_A, d_B, 0.0f, d_C); - total_ms += timer.stop(); - } - - CUDA_CHECK(cudaFree(d_A)); - CUDA_CHECK(cudaFree(d_B)); - CUDA_CHECK(cudaFree(d_C)); - - return total_ms / ITERS; -} - - -void benchmark_coalesced(int S, - const float* h_A, - const float* h_B) +static void run_sweep(const char* label, BenchmarkFn fn, + const int* sizes, int n_sizes, + const float* h_A, const float* h_B) { - size_t szA = (size_t)S * S * sizeof(float); - size_t szB = (size_t)S * S * sizeof(float); - size_t szC = (size_t)S * S * sizeof(float); - - float *d_A, *d_B, *d_C; - CUDA_CHECK(cudaMalloc(&d_A, szA)); - CUDA_CHECK(cudaMalloc(&d_B, szB)); - CUDA_CHECK(cudaMalloc(&d_C, szC)); + 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"); - CUDA_CHECK(cudaMemcpy(d_A, h_A, szA, cudaMemcpyHostToDevice)); - CUDA_CHECK(cudaMemcpy(d_B, h_B, szB, cudaMemcpyHostToDevice)); - CUDA_CHECK(cudaMemset(d_C, 0, szC)); + for (int i = 0; i < n_sizes; ++i) { + int S = sizes[i]; - dim3 block(32 * 32); - dim3 grid((S + block.x / 16 - 1) / (block.x / 16), - (S + block.x / 16 - 1) / (block.x / 16)); + float avg_ms = fn(S, h_A, h_B); - // Warmup — not measured. - for (int i = 0; i < WARMUP; ++i) - sgemm_coalesced<16><<>>(S, S, S, 1.0f, d_A, d_B, 0.0f, d_C); - CUDA_CHECK(cudaDeviceSynchronize()); + double flops = 2.0 * S * S * S; + double bytes = ((double)S * S + + (double)S * S + + 2.0 * S * S) + * sizeof(float); + double ai = flops / bytes; + double gflops = flops / (avg_ms * 1e-3) / 1e9; + double bandwidth = bytes / (avg_ms * 1e-3) / 1e9; - // Timed runs. - GpuTimer timer; - float total_ms = 0.0f; - for (int i = 0; i < ITERS; ++i) { - CUDA_CHECK(cudaMemset(d_C, 0, szC)); - timer.start(); - total_ms += timer.stop(); + printf("%-6d %9.3f %9.1f %11.1f %9.2f\n", + S, avg_ms, gflops, bandwidth, ai); } - - CUDA_CHECK(cudaFree(d_A)); - CUDA_CHECK(cudaFree(d_B)); - CUDA_CHECK(cudaFree(d_C)); - - printf("Coalesced kernel: %9.3f ms\n", total_ms / ITERS); + printf("\n"); } -// --------------------------------------------------------------------------- -// Main -// --------------------------------------------------------------------------- int main() { srand(42); - // // ── Correctness check (small fixed size) ──────────────────────────────── - // // The CPU triple-loop is O(S^3), so we only run it for a small S. - // const int S_verify = 256; - // { - // float* h_A = new float[S_verify * S_verify]; - // float* h_B = new float[S_verify * S_verify]; - // float* h_C = new float[S_verify * S_verify](); // zero-init - // float* h_C_ref = new float[S_verify * S_verify](); - - // fill_random(h_A, S_verify * S_verify); - // fill_random(h_B, S_verify * S_verify); - // matmul_cpu(h_A, h_B, h_C_ref, S_verify, S_verify, S_verify); - - // float *d_A, *d_B, *d_C; - // size_t sz = (size_t)S_verify * S_verify * sizeof(float); - // CUDA_CHECK(cudaMalloc(&d_A, sz)); - // CUDA_CHECK(cudaMalloc(&d_B, sz)); - // CUDA_CHECK(cudaMalloc(&d_C, sz)); - // CUDA_CHECK(cudaMemcpy(d_A, h_A, sz, cudaMemcpyHostToDevice)); - // CUDA_CHECK(cudaMemcpy(d_B, h_B, sz, cudaMemcpyHostToDevice)); - // CUDA_CHECK(cudaMemset(d_C, 0, sz)); - - // constexpr int BS = 32; - // dim3 block(BS * BS); // BLOCKSIZE² threads: threadIdx.x / BS = row, threadIdx.x % BS = col - // dim3 grid((S_verify + BS - 1) / BS, (S_verify + BS - 1) / BS); - // sgemm_coalesced<<>>(S_verify, S_verify, S_verify, 1.0f, d_A, d_B, 0.0f, d_C); - // CHECK_LAST_ERROR(); - // CHECK_SYNC(); - // CUDA_CHECK(cudaMemcpy(h_C, d_C, sz, cudaMemcpyDeviceToHost)); - - // printf("Correctness check (S=%d): ", S_verify); - // if (verify(h_C_ref, h_C, S_verify * S_verify)) - // printf("PASSED\n\n"); - // else - // printf("FAILED\n\n"); - - // CUDA_CHECK(cudaFree(d_A)); - // CUDA_CHECK(cudaFree(d_B)); - // CUDA_CHECK(cudaFree(d_C)); - // delete[] h_A; delete[] h_B; delete[] h_C; delete[] h_C_ref; - // } - - const int sizes[] = {128, 256, 512, 1024, 2048, 4096, 8192, 16384, 32768}; const int N_SIZES = sizeof(sizes) / sizeof(sizes[0]); - // Allocate host matrices for the largest size once. 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); - printf("Naive kernel roofline sweep (%d iters, %d warmup)\n", ITERS, WARMUP); - printf("%-6s %9s %9s %11s %9s\n", - "Size", "Time(ms)", "GFLOP/s", "BW(GB/s)", "AI(F/B)"); - printf("------ --------- --------- ----------- ---------\n"); - - for (int i = 0; i < N_SIZES; ++i) { - int S = sizes[i]; - - float avg_ms = benchmark_naive(S, h_A, h_B); - float avg_ms_tiled = benchmark_tiled(S, h_A, h_B); - - // Each term is promoted to double before multiplying to prevent int32 - // overflow at large S (see note above). - double flops = 2.0 * S * S * S; - double bytes = ((double)S * S // read A - + (double)S * S // read B - + 2.0 * S * S) // read + write C - * sizeof(float); - 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); - - // printf(" Tiled kernel: %9.3f ms\n", avg_ms_tiled); - } + run_sweep("Naive", benchmark_naive, sizes, N_SIZES, h_A, h_B); + run_sweep("Tiled", benchmark_tiled, sizes, N_SIZES, h_A, h_B); + run_sweep("Coalesced", benchmark_coalesced, sizes, N_SIZES, h_A, h_B); delete[] h_A; delete[] h_B; diff --git a/cuda-mma/src/benchmarks.cu b/cuda-mma/src/benchmarks.cu new file mode 100644 index 0000000..4599d04 --- /dev/null +++ b/cuda-mma/src/benchmarks.cu @@ -0,0 +1,125 @@ +#include + +#include "cuda_check.cuh" +#include "kernels.cuh" +#include "timer.h" +#include "utils.h" +#include "benchmarks.h" + +#define WARMUP 3 +#define ITERS 5 + +float benchmark_naive(int S, const float* h_A, const float* h_B) +{ + size_t szA = (size_t)S * S * sizeof(float); + size_t szB = (size_t)S * S * sizeof(float); + size_t szC = (size_t)S * S * sizeof(float); + + float *d_A, *d_B, *d_C; + CUDA_CHECK(cudaMalloc(&d_A, szA)); + CUDA_CHECK(cudaMalloc(&d_B, szB)); + CUDA_CHECK(cudaMalloc(&d_C, szC)); + + CUDA_CHECK(cudaMemcpy(d_A, h_A, szA, cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemcpy(d_B, h_B, szB, cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemset(d_C, 0, szC)); + + dim3 gridDim(CEIL_DIV(S, 32), CEIL_DIV(S, 32)); + dim3 blockDim(32, 32, 1); + + for (int i = 0; i < WARMUP; ++i) + sgemm_naive<<>>(S, S, S, 1.0f, d_A, d_B, 0.0f, d_C); + CUDA_CHECK(cudaDeviceSynchronize()); + + GpuTimer timer; + float total_ms = 0.0f; + for (int i = 0; i < ITERS; ++i) { + CUDA_CHECK(cudaMemset(d_C, 0, szC)); + timer.start(); + sgemm_naive<<>>(S, S, S, 1.0f, d_A, d_B, 0.0f, d_C); + total_ms += timer.stop(); + } + + CUDA_CHECK(cudaFree(d_A)); + CUDA_CHECK(cudaFree(d_B)); + CUDA_CHECK(cudaFree(d_C)); + + return total_ms / ITERS; +} + +float benchmark_tiled(int S, const float* h_A, const float* h_B) +{ + size_t szA = (size_t)S * S * sizeof(float); + size_t szB = (size_t)S * S * sizeof(float); + size_t szC = (size_t)S * S * sizeof(float); + + float *d_A, *d_B, *d_C; + CUDA_CHECK(cudaMalloc(&d_A, szA)); + CUDA_CHECK(cudaMalloc(&d_B, szB)); + CUDA_CHECK(cudaMalloc(&d_C, szC)); + + CUDA_CHECK(cudaMemcpy(d_A, h_A, szA, cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemcpy(d_B, h_B, szB, cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemset(d_C, 0, szC)); + + dim3 blockDim(32, 32, 1); + dim3 gridDim(CEIL_DIV(S, blockDim.x), CEIL_DIV(S, blockDim.y)); + + for (int i = 0; i < WARMUP; ++i) + sgemm_tiled<32><<>>(S, S, S, 1.0f, d_A, d_B, 0.0f, d_C); + CUDA_CHECK(cudaDeviceSynchronize()); + + GpuTimer timer; + float total_ms = 0.0f; + for (int i = 0; i < ITERS; ++i) { + CUDA_CHECK(cudaMemset(d_C, 0, szC)); + timer.start(); + sgemm_tiled<32><<>>(S, S, S, 1.0f, d_A, d_B, 0.0f, d_C); + total_ms += timer.stop(); + } + + CUDA_CHECK(cudaFree(d_A)); + CUDA_CHECK(cudaFree(d_B)); + CUDA_CHECK(cudaFree(d_C)); + + return total_ms / ITERS; +} + +float benchmark_coalesced(int S, const float* h_A, const float* h_B) +{ + size_t szA = (size_t)S * S * sizeof(float); + size_t szB = (size_t)S * S * sizeof(float); + size_t szC = (size_t)S * S * sizeof(float); + + float *d_A, *d_B, *d_C; + CUDA_CHECK(cudaMalloc(&d_A, szA)); + CUDA_CHECK(cudaMalloc(&d_B, szB)); + CUDA_CHECK(cudaMalloc(&d_C, szC)); + + CUDA_CHECK(cudaMemcpy(d_A, h_A, szA, cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemcpy(d_B, h_B, szB, cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemset(d_C, 0, szC)); + + dim3 block(32 * 32); + dim3 grid((S + block.x / 16 - 1) / (block.x / 16), + (S + block.x / 16 - 1) / (block.x / 16)); + + for (int i = 0; i < WARMUP; ++i) + sgemm_coalesced<16><<>>(S, S, S, 1.0f, d_A, d_B, 0.0f, d_C); + CUDA_CHECK(cudaDeviceSynchronize()); + + GpuTimer timer; + float total_ms = 0.0f; + for (int i = 0; i < ITERS; ++i) { + CUDA_CHECK(cudaMemset(d_C, 0, szC)); + timer.start(); + sgemm_coalesced<16><<>>(S, S, S, 1.0f, d_A, d_B, 0.0f, d_C); + total_ms += timer.stop(); + } + + CUDA_CHECK(cudaFree(d_A)); + CUDA_CHECK(cudaFree(d_B)); + CUDA_CHECK(cudaFree(d_C)); + + return total_ms / ITERS; +} diff --git a/cuda-mma/src/kernels.cu b/cuda-mma/src/kernels.cu index bb530c1..b77070b 100644 --- a/cuda-mma/src/kernels.cu +++ b/cuda-mma/src/kernels.cu @@ -56,37 +56,41 @@ __global__ void sgemm_coalesced( // Tiled version: 1 thread per output element, shared memory for tiles of A and B. template __global__ void sgemm_tiled( - size_t M, // Number of rows in A and C + 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) + float *C) // [M x N] row-major (in-out: C = alpha*A*B + beta*C) { - __shared__ float sA[TILE_SIZE][TILE_SIZE]; - __shared__ float sB[TILE_SIZE][TILE_SIZE]; + // Declaring the tiles + __shared__ float A_tile[TILE_SIZE * TILE_SIZE]; + __shared__ float B_tile[TILE_SIZE * TILE_SIZE]; - int row = blockIdx.y * TILE_SIZE + threadIdx.y; - int col = blockIdx.x * TILE_SIZE + threadIdx.x; + // the element of C for this specific thread + size_t row = blockIdx.y * TILE_SIZE + threadIdx.y; + size_t col = blockIdx.x * TILE_SIZE + threadIdx.x; - float acc = 0.0f; + // we multiply along the K dimension, iterate over this + size_t n_tiles = (K + TILE_SIZE -1)/TILE_SIZE; + + float acc = .0f; + for(size_t tile=0; tile +// __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) +// { +// __shared__ float sA[TILE_SIZE][TILE_SIZE]; +// __shared__ float sB[TILE_SIZE][TILE_SIZE]; + +// int row = blockIdx.y * TILE_SIZE + threadIdx.y; +// int col = blockIdx.x * TILE_SIZE + threadIdx.x; + +// float acc = 0.0f; + +// // Iterate over tiles along the K dimension. +// for (size_t t = 0; t < (K + TILE_SIZE - 1) / TILE_SIZE; ++t) { +// size_t aCol = t * TILE_SIZE + threadIdx.x; // column of A this thread loads +// size_t bRow = t * TILE_SIZE + threadIdx.y; // row of B this thread loads + +// // Boundary-safe loads: pad with 0 for out-of-bounds tiles. +// sA[threadIdx.y][threadIdx.x] = (row < M && aCol < K) ? A[row * K + aCol] : 0.0f; +// sB[threadIdx.y][threadIdx.x] = (bRow < K && col < N) ? B[bRow * N + col] : 0.0f; +// __syncthreads(); + +// // Accumulate partial dot product from shared memory. +// #pragma unroll +// for (int i = 0; i < TILE_SIZE; ++i) +// acc += sA[threadIdx.y][i] * sB[i][threadIdx.x]; +// __syncthreads(); +// } + +// if (row < M && col < N) +// C[row * N + col] = alpha * acc + beta * C[row * N + col]; +// } + template __global__ void sgemm_tiled<16>(size_t M, size_t N, size_t K, float alpha, const float *A, const float *B, float beta, float *C); diff --git a/cuda-mma/test/test_correctness.cu b/cuda-mma/test/test_correctness.cu new file mode 100644 index 0000000..9ec89bf --- /dev/null +++ b/cuda-mma/test/test_correctness.cu @@ -0,0 +1,50 @@ +#include +#include + +#include "cuda_check.cuh" +#include "kernels.cuh" +#include "utils.h" + +int main() { + srand(42); + + const int S = 256; + float* h_A = new float[S * S]; + float* h_B = new float[S * S]; + float* h_C = new float[S * S](); + float* h_C_ref = new float[S * S](); + + fill_random(h_A, S * S); + fill_random(h_B, S * S); + matmul_cpu(h_A, h_B, h_C_ref, S, S, S); + + float *d_A, *d_B, *d_C; + size_t sz = (size_t)S * S * sizeof(float); + CUDA_CHECK(cudaMalloc(&d_A, sz)); + CUDA_CHECK(cudaMalloc(&d_B, sz)); + CUDA_CHECK(cudaMalloc(&d_C, sz)); + CUDA_CHECK(cudaMemcpy(d_A, h_A, sz, cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemcpy(d_B, h_B, sz, cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemset(d_C, 0, sz)); + + constexpr int BS = 32; + dim3 block(BS, BS); + dim3 grid((S + BS - 1) / BS, (S + BS - 1) / BS); + sgemm_tiled<<>>(S, S, S, 1.0f, d_A, d_B, 0.0f, d_C); + CHECK_LAST_ERROR(); + CHECK_SYNC(); + CUDA_CHECK(cudaMemcpy(h_C, d_C, sz, cudaMemcpyDeviceToHost)); + + printf("Correctness check (S=%d): ", S); + if (verify(h_C_ref, h_C, S * S)) + printf("PASSED\n"); + else + printf("FAILED\n"); + + CUDA_CHECK(cudaFree(d_A)); + CUDA_CHECK(cudaFree(d_B)); + CUDA_CHECK(cudaFree(d_C)); + delete[] h_A; delete[] h_B; delete[] h_C; delete[] h_C_ref; + + return 0; +} From 7b8b9b10c49332c88d8e247d72ba410db419dcbd Mon Sep 17 00:00:00 2001 From: Sebas <12949040+SebastiaAgramunt@users.noreply.github.com> Date: Sat, 18 Apr 2026 07:35:03 +0000 Subject: [PATCH 04/21] adding csvs --- cuda-mma/.gitignore | 1 + cuda-mma/main.cu | 18 ++++++++--- cuda-mma/scripts/plot_results.py | 54 ++++++++++++++++++++++++++++++++ 3 files changed, 69 insertions(+), 4 deletions(-) create mode 100644 cuda-mma/scripts/plot_results.py diff --git a/cuda-mma/.gitignore b/cuda-mma/.gitignore index 567609b..b97a986 100644 --- a/cuda-mma/.gitignore +++ b/cuda-mma/.gitignore @@ -1 +1,2 @@ build/ +output/ diff --git a/cuda-mma/main.cu b/cuda-mma/main.cu index ba9e25d..7fcb2d8 100644 --- a/cuda-mma/main.cu +++ b/cuda-mma/main.cu @@ -1,5 +1,6 @@ #include #include +#include #include "utils.h" #include "benchmarks.h" @@ -11,13 +12,17 @@ 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 float* h_A, const float* h_B, + const char* csv_path) { 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]; @@ -34,13 +39,18 @@ static void run_sweep(const char* label, BenchmarkFn fn, 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]); @@ -50,9 +60,9 @@ int main() { 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); - run_sweep("Tiled", benchmark_tiled, sizes, N_SIZES, h_A, h_B); - run_sweep("Coalesced", benchmark_coalesced, sizes, N_SIZES, h_A, h_B); + 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"); delete[] h_A; delete[] h_B; diff --git a/cuda-mma/scripts/plot_results.py b/cuda-mma/scripts/plot_results.py new file mode 100644 index 0000000..2386147 --- /dev/null +++ b/cuda-mma/scripts/plot_results.py @@ -0,0 +1,54 @@ +#!/usr/bin/env python3 +import os +import pandas as pd +import matplotlib.pyplot as plt +import matplotlib.ticker as ticker + +OUTPUT_DIR = os.path.join(os.path.dirname(__file__), "..", "output") + +kernels = { + "Naive": "naive.csv", + "Tiled": "tiled.csv", + "Coalesced": "coalesced.csv", +} + +def load(filename): + path = os.path.join(OUTPUT_DIR, filename) + if not os.path.exists(path): + return None + return pd.read_csv(path) + +fig, axes = plt.subplots(1, 3, figsize=(15, 5)) +ax_gflops, ax_bw, ax_ai = axes + +for label, filename in kernels.items(): + df = load(filename) + if df is None: + print(f"Warning: {filename} not found, skipping.") + continue + ax_gflops.plot(df["size"], df["gflops"], marker="o", label=label) + ax_bw.plot (df["size"], df["bandwidth_gbs"], marker="o", label=label) + ax_ai.plot (df["size"], df["arithmetic_intensity"], marker="o", label=label) + +for ax in axes: + ax.set_xscale("log", base=2) + ax.set_yscale("log") + ax.xaxis.set_major_formatter(ticker.FuncFormatter(lambda x, _: f"{int(x)}")) + ax.set_xlabel("Matrix size (S×S)") + ax.legend() + ax.grid(True, which="both", linestyle="--", linewidth=0.5) + +ax_gflops.set_title("Throughput") +ax_gflops.set_ylabel("GFLOP/s") + +ax_bw.set_title("Memory Bandwidth") +ax_bw.set_ylabel("GB/s") + +ax_ai.set_title("Arithmetic Intensity") +ax_ai.set_ylabel("FLOP/Byte") + +fig.tight_layout() +out_path = os.path.join(OUTPUT_DIR, "roofline.png") +plt.savefig(out_path, dpi=150) +print(f"Saved {out_path}") +plt.show() From 4d382ca008ebfcec70c2cd2e3d72488d7cad18f5 Mon Sep 17 00:00:00 2001 From: Sebas <12949040+SebastiaAgramunt@users.noreply.github.com> Date: Sat, 18 Apr 2026 18:02:53 +0000 Subject: [PATCH 05/21] benchmarks --- cuda-mma/scripts/plot_results.py | 11 ++++++----- cuda-mma/src/benchmarks.cu | 9 ++++----- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/cuda-mma/scripts/plot_results.py b/cuda-mma/scripts/plot_results.py index 2386147..a2c1082 100644 --- a/cuda-mma/scripts/plot_results.py +++ b/cuda-mma/scripts/plot_results.py @@ -1,10 +1,10 @@ #!/usr/bin/env python3 -import os +from path import Path import pandas as pd import matplotlib.pyplot as plt import matplotlib.ticker as ticker -OUTPUT_DIR = os.path.join(os.path.dirname(__file__), "..", "output") +OUTPUT_DIR = Path(__file__).parent.parent / "output" kernels = { "Naive": "naive.csv", @@ -13,12 +13,13 @@ } def load(filename): - path = os.path.join(OUTPUT_DIR, filename) - if not os.path.exists(path): + path = OUTPUT_DIR / filename + if not path.exists(): return None return pd.read_csv(path) fig, axes = plt.subplots(1, 3, figsize=(15, 5)) +fig.suptitle("CUDA SGEMM Roofline Analysis - Nvidia GPU A10") ax_gflops, ax_bw, ax_ai = axes for label, filename in kernels.items(): @@ -48,7 +49,7 @@ def load(filename): ax_ai.set_ylabel("FLOP/Byte") fig.tight_layout() -out_path = os.path.join(OUTPUT_DIR, "roofline.png") +out_path = OUTPUT_DIR / "roofline.png" plt.savefig(out_path, dpi=150) print(f"Saved {out_path}") plt.show() diff --git a/cuda-mma/src/benchmarks.cu b/cuda-mma/src/benchmarks.cu index 4599d04..a67905b 100644 --- a/cuda-mma/src/benchmarks.cu +++ b/cuda-mma/src/benchmarks.cu @@ -100,12 +100,11 @@ float benchmark_coalesced(int S, const float* h_A, const float* h_B) CUDA_CHECK(cudaMemcpy(d_B, h_B, szB, cudaMemcpyHostToDevice)); CUDA_CHECK(cudaMemset(d_C, 0, szC)); - dim3 block(32 * 32); - dim3 grid((S + block.x / 16 - 1) / (block.x / 16), - (S + block.x / 16 - 1) / (block.x / 16)); + dim3 gridDim(CEIL_DIV(S, 32), CEIL_DIV(S, 32)); + dim3 blockDim(32 * 32); for (int i = 0; i < WARMUP; ++i) - sgemm_coalesced<16><<>>(S, S, S, 1.0f, d_A, d_B, 0.0f, d_C); + sgemm_coalesced<32><<>>(S, S, S, 1.0f, d_A, d_B, 0.0f, d_C); CUDA_CHECK(cudaDeviceSynchronize()); GpuTimer timer; @@ -113,7 +112,7 @@ float benchmark_coalesced(int S, const float* h_A, const float* h_B) for (int i = 0; i < ITERS; ++i) { CUDA_CHECK(cudaMemset(d_C, 0, szC)); timer.start(); - sgemm_coalesced<16><<>>(S, S, S, 1.0f, d_A, d_B, 0.0f, d_C); + sgemm_coalesced<32><<>>(S, S, S, 1.0f, d_A, d_B, 0.0f, d_C); total_ms += timer.stop(); } From ddf33fa2c006d649bd03abdabe6d3c5d2192f2fa Mon Sep 17 00:00:00 2001 From: Sebas <12949040+SebastiaAgramunt@users.noreply.github.com> Date: Tue, 21 Apr 2026 04:53:20 +0000 Subject: [PATCH 06/21] addig cuBLAS and cBLAS --- cuda-mma/Makefile | 2 +- cuda-mma/include/benchmarks.h | 2 + cuda-mma/main.cu | 12 +++- cuda-mma/scripts/install_env.sh | 8 +++ cuda-mma/scripts/plot_results.py | 28 +++++++--- cuda-mma/src/benchmarks.cu | 95 ++++++++++++++++++++++++++++++++ 6 files changed, 134 insertions(+), 13 deletions(-) create mode 100755 cuda-mma/scripts/install_env.sh diff --git a/cuda-mma/Makefile b/cuda-mma/Makefile index 70e8cb7..7187e4c 100644 --- a/cuda-mma/Makefile +++ b/cuda-mma/Makefile @@ -33,7 +33,7 @@ $(BIN_DIR): mkdir -p $(BIN_DIR) $(TARGET): $(SRCS) $(HDRS) | $(BIN_DIR) - $(NVCC) $(NVCC_FLAGS) -o $@ $(SRCS) + $(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) diff --git a/cuda-mma/include/benchmarks.h b/cuda-mma/include/benchmarks.h index 6aa05b8..fca2cba 100644 --- a/cuda-mma/include/benchmarks.h +++ b/cuda-mma/include/benchmarks.h @@ -3,3 +3,5 @@ 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); diff --git a/cuda-mma/main.cu b/cuda-mma/main.cu index 7fcb2d8..fa6b0c9 100644 --- a/cuda-mma/main.cu +++ b/cuda-mma/main.cu @@ -60,9 +60,15 @@ int main() { 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("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"); + + // CBLAS runs on CPU; cap at 4096 to keep runtime reasonable + const int cblas_sizes[] = {128, 256, 512, 1024, 2048, 4096}; + 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; diff --git a/cuda-mma/scripts/install_env.sh b/cuda-mma/scripts/install_env.sh new file mode 100755 index 0000000..3dc0b81 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/cuda-mma/scripts/plot_results.py b/cuda-mma/scripts/plot_results.py index a2c1082..e2d0a59 100644 --- a/cuda-mma/scripts/plot_results.py +++ b/cuda-mma/scripts/plot_results.py @@ -10,6 +10,8 @@ "Naive": "naive.csv", "Tiled": "tiled.csv", "Coalesced": "coalesced.csv", + "cuBLAS": "cublas.csv", + "CBLAS": "cblas.csv", } def load(filename): @@ -18,19 +20,25 @@ def load(filename): return None return pd.read_csv(path) -fig, axes = plt.subplots(1, 3, figsize=(15, 5)) -fig.suptitle("CUDA SGEMM Roofline Analysis - Nvidia GPU A10") -ax_gflops, ax_bw, ax_ai = axes +fig, axes = plt.subplots(1, 2, figsize=(12, 5)) +fig.suptitle("CUDA SGEMM Roofline Analysis - Nvidia GPU A6000") +ax_gflops, ax_bw = axes for label, filename in kernels.items(): df = load(filename) if df is None: print(f"Warning: {filename} not found, skipping.") continue - ax_gflops.plot(df["size"], df["gflops"], marker="o", label=label) - ax_bw.plot (df["size"], df["bandwidth_gbs"], marker="o", label=label) - ax_ai.plot (df["size"], df["arithmetic_intensity"], marker="o", label=label) + ax_gflops.plot(df["size"], df["gflops"], marker="o", label=label) + ax_bw.plot (df["size"], df["bandwidth_gbs"], marker="o", label=label) +def ai_fmt(x, _): + if x >= 1e9: return f"{x/1e9:.3g}G" + if x >= 1e6: return f"{x/1e6:.3g}M" + if x >= 1e3: return f"{x/1e3:.3g}k" + return f"{x:.3g}" + +# AI = S/8 FLOP/Byte = S/8 * 1e6 FLOP/MB (2S³ flops / 16S² bytes) for ax in axes: ax.set_xscale("log", base=2) ax.set_yscale("log") @@ -39,15 +47,17 @@ def load(filename): ax.legend() ax.grid(True, which="both", linestyle="--", linewidth=0.5) + sec = ax.secondary_xaxis("top", functions=(lambda x: x / 8 * 1e6, lambda x: x * 8 / 1e6)) + sec.set_xscale("log", base=2) + sec.xaxis.set_major_formatter(ticker.FuncFormatter(ai_fmt)) + sec.set_xlabel("Arithmetic Intensity (FLOP/MB)") + ax_gflops.set_title("Throughput") ax_gflops.set_ylabel("GFLOP/s") ax_bw.set_title("Memory Bandwidth") ax_bw.set_ylabel("GB/s") -ax_ai.set_title("Arithmetic Intensity") -ax_ai.set_ylabel("FLOP/Byte") - fig.tight_layout() out_path = OUTPUT_DIR / "roofline.png" plt.savefig(out_path, dpi=150) diff --git a/cuda-mma/src/benchmarks.cu b/cuda-mma/src/benchmarks.cu index a67905b..07e67b8 100644 --- a/cuda-mma/src/benchmarks.cu +++ b/cuda-mma/src/benchmarks.cu @@ -1,11 +1,24 @@ #include +#include +#include +#include #include "cuda_check.cuh" #include "kernels.cuh" #include "timer.h" #include "utils.h" #include "benchmarks.h" +#define CUBLAS_CHECK(call) \ + do { \ + cublasStatus_t _s = (call); \ + if (_s != CUBLAS_STATUS_SUCCESS) { \ + fprintf(stderr, "cuBLAS error %d at %s:%d\n", _s, \ + __FILE__, __LINE__); \ + exit(EXIT_FAILURE); \ + } \ + } while (0) + #define WARMUP 3 #define ITERS 5 @@ -122,3 +135,85 @@ float benchmark_coalesced(int S, const float* h_A, const float* h_B) return total_ms / ITERS; } + +float benchmark_cublas(int S, const float* h_A, const float* h_B) +{ + size_t szA = (size_t)S * S * sizeof(float); + size_t szB = (size_t)S * S * sizeof(float); + size_t szC = (size_t)S * S * sizeof(float); + + float *d_A, *d_B, *d_C; + CUDA_CHECK(cudaMalloc(&d_A, szA)); + CUDA_CHECK(cudaMalloc(&d_B, szB)); + CUDA_CHECK(cudaMalloc(&d_C, szC)); + + CUDA_CHECK(cudaMemcpy(d_A, h_A, szA, cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemcpy(d_B, h_B, szB, cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemset(d_C, 0, szC)); + + cublasHandle_t handle; + CUBLAS_CHECK(cublasCreate(&handle)); + + // cuBLAS is column-major. For row-major C=A*B, use the identity + // C^T = B^T * A^T, so pass B first with leading dimension S. + const float alpha = 1.0f, beta = 0.0f; + + for (int i = 0; i < WARMUP; ++i) + CUBLAS_CHECK(cublasSgemm(handle, + CUBLAS_OP_N, CUBLAS_OP_N, + S, S, S, + &alpha, + d_B, S, + d_A, S, + &beta, + d_C, S)); + CUDA_CHECK(cudaDeviceSynchronize()); + + GpuTimer timer; + float total_ms = 0.0f; + for (int i = 0; i < ITERS; ++i) { + CUDA_CHECK(cudaMemset(d_C, 0, szC)); + timer.start(); + CUBLAS_CHECK(cublasSgemm(handle, + CUBLAS_OP_N, CUBLAS_OP_N, + S, S, S, + &alpha, + d_B, S, + d_A, S, + &beta, + d_C, S)); + total_ms += timer.stop(); + } + + CUBLAS_CHECK(cublasDestroy(handle)); + CUDA_CHECK(cudaFree(d_A)); + CUDA_CHECK(cudaFree(d_B)); + CUDA_CHECK(cudaFree(d_C)); + + return total_ms / ITERS; +} + +float benchmark_cblas(int S, const float* h_A, const float* h_B) +{ + size_t sz = (size_t)S * S; + float* h_C = new float[sz](); + + using clock = std::chrono::high_resolution_clock; + + for (int i = 0; i < WARMUP; ++i) + cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, + S, S, S, 1.0f, h_A, S, h_B, S, 0.0f, h_C, S); + + float total_ms = 0.0f; + for (int i = 0; i < ITERS; ++i) { + std::fill(h_C, h_C + sz, 0.0f); + auto t0 = clock::now(); + cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, + S, S, S, 1.0f, h_A, S, h_B, S, 0.0f, h_C, S); + auto t1 = clock::now(); + total_ms += std::chrono::duration(t1 - t0).count(); + } + + delete[] h_C; + return total_ms / ITERS; +} From a17281d8271822611ae3a17582cc8b1365aaddbb Mon Sep 17 00:00:00 2001 From: Sebas <12949040+SebastiaAgramunt@users.noreply.github.com> Date: Tue, 21 Apr 2026 05:21:38 +0000 Subject: [PATCH 07/21] reverting changes --- cuda-mma/main.cu | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/cuda-mma/main.cu b/cuda-mma/main.cu index fa6b0c9..fa8db0e 100644 --- a/cuda-mma/main.cu +++ b/cuda-mma/main.cu @@ -60,15 +60,15 @@ int main() { 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("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"); // CBLAS runs on CPU; cap at 4096 to keep runtime reasonable - const int cblas_sizes[] = {128, 256, 512, 1024, 2048, 4096}; + 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"); + run_sweep("CBLAS", benchmark_cblas, cblas_sizes, N_CBLAS, h_A, h_B, "output/cblas.csv"); delete[] h_A; delete[] h_B; From 87371c5b267fcf8f5340dd642311b2f40a3439d5 Mon Sep 17 00:00:00 2001 From: Sebas <12949040+SebastiaAgramunt@users.noreply.github.com> Date: Wed, 22 Apr 2026 04:45:16 +0000 Subject: [PATCH 08/21] adding coarsened threads --- cuda-mma/include/benchmarks.h | 1 + cuda-mma/include/kernels.cuh | 14 +++++++ cuda-mma/main.cu | 1 + cuda-mma/scripts/plot_results.py | 1 + cuda-mma/src/benchmarks.cu | 40 +++++++++++++++++++ cuda-mma/src/kernels.cu | 68 ++++++++++++++++++++++++++++++++ 6 files changed, 125 insertions(+) diff --git a/cuda-mma/include/benchmarks.h b/cuda-mma/include/benchmarks.h index fca2cba..8d3e317 100644 --- a/cuda-mma/include/benchmarks.h +++ b/cuda-mma/include/benchmarks.h @@ -5,3 +5,4 @@ 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); diff --git a/cuda-mma/include/kernels.cuh b/cuda-mma/include/kernels.cuh index 2fdf4a6..797faa7 100644 --- a/cuda-mma/include/kernels.cuh +++ b/cuda-mma/include/kernels.cuh @@ -35,4 +35,18 @@ __global__ void sgemm_tiled( 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 +__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) \ No newline at end of file diff --git a/cuda-mma/main.cu b/cuda-mma/main.cu index fa8db0e..9c13db3 100644 --- a/cuda-mma/main.cu +++ b/cuda-mma/main.cu @@ -64,6 +64,7 @@ int main() { 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"); // CBLAS runs on CPU; cap at 4096 to keep runtime reasonable const int cblas_sizes[] = {128, 256, 512, 1024, 2048, 4096, 8192}; diff --git a/cuda-mma/scripts/plot_results.py b/cuda-mma/scripts/plot_results.py index e2d0a59..1ad6ae2 100644 --- a/cuda-mma/scripts/plot_results.py +++ b/cuda-mma/scripts/plot_results.py @@ -11,6 +11,7 @@ "Tiled": "tiled.csv", "Coalesced": "coalesced.csv", "cuBLAS": "cublas.csv", + "Coarsened": "coarsened.csv", "CBLAS": "cblas.csv", } diff --git a/cuda-mma/src/benchmarks.cu b/cuda-mma/src/benchmarks.cu index 07e67b8..688fac0 100644 --- a/cuda-mma/src/benchmarks.cu +++ b/cuda-mma/src/benchmarks.cu @@ -217,3 +217,43 @@ float benchmark_cblas(int S, const float* h_A, const float* h_B) delete[] h_C; return total_ms / ITERS; } + +float benchmark_coarsened(int S, const float* h_A, const float* h_B) +{ + constexpr int BM = 64, BN = 64, BK = 8, TM = 8; + + size_t szA = (size_t)S * S * sizeof(float); + size_t szB = (size_t)S * S * sizeof(float); + size_t szC = (size_t)S * S * sizeof(float); + + float *d_A, *d_B, *d_C; + CUDA_CHECK(cudaMalloc(&d_A, szA)); + CUDA_CHECK(cudaMalloc(&d_B, szB)); + CUDA_CHECK(cudaMalloc(&d_C, szC)); + + CUDA_CHECK(cudaMemcpy(d_A, h_A, szA, cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemcpy(d_B, h_B, szB, cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemset(d_C, 0, szC)); + + dim3 gridDim(CEIL_DIV(S, BN), CEIL_DIV(S, BM)); + dim3 blockDim(BM * BN / TM); // 512 threads + + for (int i = 0; i < WARMUP; ++i) + sgemm_coarsened<<>>(S, S, S, 1.0f, d_A, d_B, 0.0f, d_C); + CUDA_CHECK(cudaDeviceSynchronize()); + + GpuTimer timer; + float total_ms = 0.0f; + for (int i = 0; i < ITERS; ++i) { + CUDA_CHECK(cudaMemset(d_C, 0, szC)); + timer.start(); + sgemm_coarsened<<>>(S, S, S, 1.0f, d_A, d_B, 0.0f, d_C); + total_ms += timer.stop(); + } + + CUDA_CHECK(cudaFree(d_A)); + CUDA_CHECK(cudaFree(d_B)); + CUDA_CHECK(cudaFree(d_C)); + + return total_ms / ITERS; +} diff --git a/cuda-mma/src/kernels.cu b/cuda-mma/src/kernels.cu index b77070b..0d6709a 100644 --- a/cuda-mma/src/kernels.cu +++ b/cuda-mma/src/kernels.cu @@ -187,3 +187,71 @@ template __global__ void sgemm_coalesced<16>(size_t M, size_t N, size_t K, float template __global__ void sgemm_coalesced<32>(size_t M, size_t N, size_t K, float alpha, const float *A, const float *B, float beta, float *C); + + +// Thread-coarsened tiled SGEMM. +// Block tile BM×BN, K-step BK, each thread handles TM rows. +// Constraint: BM*BK == blockDim.x (512 with defaults BM=64,BN=64,BK=8,TM=8) +// i.e. BM*BN/TM threads, each loading one element of As and one of Bs. +template +__global__ void sgemm_coarsened( + size_t M, size_t N, size_t K, + float alpha, const float *A, const float *B, + float beta, float *C) +{ + const int blockRow = blockIdx.y; + const int blockCol = blockIdx.x; + + // Each thread's output position within the block tile + const int threadRow = threadIdx.x / BN; // row group [0, BM/TM) + const int threadCol = threadIdx.x % BN; // column [0, BN) + + // Loading indices: each thread loads one element into As and one into Bs + const int innerRowA = threadIdx.x / BK; // [0, BM) + const int innerColA = threadIdx.x % BK; // [0, BK) + const int innerRowB = threadIdx.x / BN; // [0, BK) + const int innerColB = threadIdx.x % BN; // [0, BN) + + __shared__ float As[BM * BK]; + __shared__ float Bs[BK * BN]; + + float threadResults[TM] = {}; + + for (int tileIdx = 0; tileIdx < (int)((K + BK - 1) / BK); ++tileIdx) { + // Load one element of A tile + size_t gRowA = (size_t)blockRow * BM + innerRowA; + size_t gColA = (size_t)tileIdx * BK + innerColA; + As[innerRowA * BK + innerColA] = (gRowA < M && gColA < K) ? A[gRowA * K + gColA] : 0.0f; + + // Load one element of B tile + size_t gRowB = (size_t)tileIdx * BK + innerRowB; + size_t gColB = (size_t)blockCol * BN + innerColB; + Bs[innerRowB * BN + innerColB] = (gRowB < K && gColB < N) ? B[gRowB * N + gColB] : 0.0f; + + __syncthreads(); + + // Accumulate: load Bs value once, reuse across all TM rows + #pragma unroll + for (int dotIdx = 0; dotIdx < BK; ++dotIdx) { + float bVal = Bs[dotIdx * BN + threadCol]; + #pragma unroll + for (int tm = 0; tm < TM; ++tm) + threadResults[tm] += As[(threadRow * TM + tm) * BK + dotIdx] * bVal; + } + + __syncthreads(); + } + + // Write TM results to global memory + #pragma unroll + for (int tm = 0; tm < TM; ++tm) { + size_t gRow = (size_t)blockRow * BM + threadRow * TM + tm; + size_t gCol = (size_t)blockCol * BN + threadCol; + if (gRow < M && gCol < N) + C[gRow * N + gCol] = alpha * threadResults[tm] + beta * C[gRow * N + gCol]; + } +} + +template __global__ void sgemm_coarsened<64, 64, 8, 8>(size_t M, size_t N, size_t K, + float alpha, const float *A, + const float *B, float beta, float *C); From c50a504d006199c66a5550a31737a97e44a4fb3f Mon Sep 17 00:00:00 2001 From: Sebas <12949040+SebastiaAgramunt@users.noreply.github.com> Date: Wed, 22 Apr 2026 06:38:52 +0000 Subject: [PATCH 09/21] plot result --- cuda-mma/scripts/plot_results.py | 57 ++++++++++++++++++++++++++++---- 1 file changed, 51 insertions(+), 6 deletions(-) diff --git a/cuda-mma/scripts/plot_results.py b/cuda-mma/scripts/plot_results.py index 1ad6ae2..e68c674 100644 --- a/cuda-mma/scripts/plot_results.py +++ b/cuda-mma/scripts/plot_results.py @@ -7,12 +7,21 @@ OUTPUT_DIR = Path(__file__).parent.parent / "output" kernels = { + # "cBLAS": "cblas.csv", "Naive": "naive.csv", - "Tiled": "tiled.csv", "Coalesced": "coalesced.csv", - "cuBLAS": "cublas.csv", + "Tiled": "tiled.csv", "Coarsened": "coarsened.csv", - "CBLAS": "cblas.csv", + "cuBLAS": "cublas.csv", +} + +colors = { + "Naive": "C0", + "Coalesced": "C1", + "Tiled": "C2", + "Coarsened": "C3", + "cuBLAS": "C4", + "cBLAS": "C5", # C0–C4 used above } def load(filename): @@ -22,7 +31,7 @@ def load(filename): return pd.read_csv(path) fig, axes = plt.subplots(1, 2, figsize=(12, 5)) -fig.suptitle("CUDA SGEMM Roofline Analysis - Nvidia GPU A6000") +fig.suptitle("CUDA SGEMM Roofline Analysis - Nvidia GPU A100") ax_gflops, ax_bw = axes for label, filename in kernels.items(): @@ -30,8 +39,8 @@ def load(filename): if df is None: print(f"Warning: {filename} not found, skipping.") continue - ax_gflops.plot(df["size"], df["gflops"], marker="o", label=label) - ax_bw.plot (df["size"], df["bandwidth_gbs"], marker="o", label=label) + ax_gflops.plot(df["size"], df["gflops"], marker="o", label=label, color=colors[label]) + ax_bw.plot (df["size"], df["bandwidth_gbs"], marker="o", label=label, color=colors[label]) def ai_fmt(x, _): if x >= 1e9: return f"{x/1e9:.3g}G" @@ -64,3 +73,39 @@ def ai_fmt(x, _): plt.savefig(out_path, dpi=150) print(f"Saved {out_path}") plt.show() + +# --- Naive vs cBLAS throughput comparison --- +naive_vs_cblas = { + "Naive": "naive.csv", + "cBLAS": "cblas.csv", +} + +fig2, ax2 = plt.subplots(figsize=(7, 5)) +fig2.suptitle("SGEMM Roofline Analysis - Nvidia GPU A100 vs AMD EPYC 7J13 64-Core") + +for label, filename in naive_vs_cblas.items(): + df = load(filename) + if df is None: + print(f"Warning: {filename} not found, skipping.") + continue + ax2.plot(df["size"], df["gflops"], marker="o", label=label, color=colors[label]) + +ax2.set_xscale("log", base=2) +ax2.set_yscale("log") +ax2.xaxis.set_major_formatter(ticker.FuncFormatter(lambda x, _: f"{int(x)}")) +ax2.set_xlabel("Matrix size (S×S)") +ax2.set_ylabel("GFLOP/s") +ax2.set_title("Throughput") +ax2.legend() +ax2.grid(True, which="both", linestyle="--", linewidth=0.5) + +sec2 = ax2.secondary_xaxis("top", functions=(lambda x: x / 8 * 1e6, lambda x: x * 8 / 1e6)) +sec2.set_xscale("log", base=2) +sec2.xaxis.set_major_formatter(ticker.FuncFormatter(ai_fmt)) +sec2.set_xlabel("Arithmetic Intensity (FLOP/MB)") + +fig2.tight_layout() +out_path2 = OUTPUT_DIR / "roofline_naive_cblas.png" +plt.savefig(out_path2, dpi=150) +print(f"Saved {out_path2}") +plt.show() From 9f68d17dbcd19be59a57a8593c232836cb37de84 Mon Sep 17 00:00:00 2001 From: Sebas <12949040+SebastiaAgramunt@users.noreply.github.com> Date: Wed, 29 Apr 2026 05:11:43 +0000 Subject: [PATCH 10/21] cutlass et al --- cuda-mma/.gitignore | 1 + cuda-mma/Makefile | 7 ++- cuda-mma/include/benchmarks.h | 2 + cuda-mma/main.cu | 1 + cuda-mma/scripts/download_vendor.sh | 28 +++++++++ cuda-mma/scripts/plot_results.py | 4 +- cuda-mma/src/benchmark_cutlass.cu | 88 +++++++++++++++++++++++++++++ cuda-mma/src/benchmarks.cu | 16 +++++- cuda-mma/src/kernels.cu | 77 ------------------------- 9 files changed, 141 insertions(+), 83 deletions(-) create mode 100755 cuda-mma/scripts/download_vendor.sh create mode 100644 cuda-mma/src/benchmark_cutlass.cu diff --git a/cuda-mma/.gitignore b/cuda-mma/.gitignore index b97a986..ee8b6d5 100644 --- a/cuda-mma/.gitignore +++ b/cuda-mma/.gitignore @@ -1,2 +1,3 @@ build/ output/ +vendor/ diff --git a/cuda-mma/Makefile b/cuda-mma/Makefile index 7187e4c..b80b142 100644 --- a/cuda-mma/Makefile +++ b/cuda-mma/Makefile @@ -11,16 +11,19 @@ CXX_FLAGS := -std=c++17 -O3 # 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 \ - -I. -Iinclude + -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 +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 diff --git a/cuda-mma/include/benchmarks.h b/cuda-mma/include/benchmarks.h index 8d3e317..329ab74 100644 --- a/cuda-mma/include/benchmarks.h +++ b/cuda-mma/include/benchmarks.h @@ -6,3 +6,5 @@ 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(int S, const float* h_A, const float* h_B); + diff --git a/cuda-mma/main.cu b/cuda-mma/main.cu index 9c13db3..185a997 100644 --- a/cuda-mma/main.cu +++ b/cuda-mma/main.cu @@ -65,6 +65,7 @@ int main() { 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", benchmark_cutlass, sizes, N_SIZES, h_A, h_B, "output/cutlass.csv"); // CBLAS runs on CPU; cap at 4096 to keep runtime reasonable const int cblas_sizes[] = {128, 256, 512, 1024, 2048, 4096, 8192}; diff --git a/cuda-mma/scripts/download_vendor.sh b/cuda-mma/scripts/download_vendor.sh new file mode 100755 index 0000000..773bd02 --- /dev/null +++ b/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." diff --git a/cuda-mma/scripts/plot_results.py b/cuda-mma/scripts/plot_results.py index e68c674..ccff442 100644 --- a/cuda-mma/scripts/plot_results.py +++ b/cuda-mma/scripts/plot_results.py @@ -13,6 +13,7 @@ "Tiled": "tiled.csv", "Coarsened": "coarsened.csv", "cuBLAS": "cublas.csv", + "CUTLASS": "cutlass.csv", } colors = { @@ -21,7 +22,8 @@ "Tiled": "C2", "Coarsened": "C3", "cuBLAS": "C4", - "cBLAS": "C5", # C0–C4 used above + "cBLAS": "C6", + "CUTLASS": "C5", } def load(filename): diff --git a/cuda-mma/src/benchmark_cutlass.cu b/cuda-mma/src/benchmark_cutlass.cu new file mode 100644 index 0000000..84ffb3f --- /dev/null +++ b/cuda-mma/src/benchmark_cutlass.cu @@ -0,0 +1,88 @@ +#include +#include + +#include + +#include "cuda_check.cuh" +#include "timer.h" +#include "utils.h" +#include "benchmarks.h" + +#define WARMUP 3 +#define ITERS 5 + +// SGEMM via CUTLASS using TF32 Tensor Cores on Ampere (sm_80). +// CUTLASS 4.x requires tfloat32_t element types to select TF32 MMA instructions; +// float inputs are bitcast to tf32 precision (rounds the mantissa to 10 bits). +using CutlassGemm = cutlass::gemm::device::Gemm< + cutlass::tfloat32_t, cutlass::layout::RowMajor, // A (tf32 ~ float with 10-bit mantissa) + cutlass::tfloat32_t, cutlass::layout::RowMajor, // B + float, cutlass::layout::RowMajor, // C / D (full float output) + float, // accumulator + cutlass::arch::OpClassTensorOp, // TF32 tensor cores + cutlass::arch::Sm80, // Ampere + cutlass::gemm::GemmShape<128, 128, 32>, // threadblock tile + cutlass::gemm::GemmShape<64, 64, 32>, // warp tile + cutlass::gemm::GemmShape<16, 8, 8> // TF32 MMA instruction shape +>; + +float benchmark_cutlass(int S, const float* h_A, const float* h_B) +{ + size_t sz = (size_t)S * S * sizeof(float); + + float *d_A, *d_B, *d_C; + CUDA_CHECK(cudaMalloc(&d_A, sz)); + CUDA_CHECK(cudaMalloc(&d_B, sz)); + CUDA_CHECK(cudaMalloc(&d_C, sz)); + + CUDA_CHECK(cudaMemcpy(d_A, h_A, sz, cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemcpy(d_B, h_B, sz, cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemset(d_C, 0, sz)); + + const float alpha = 1.0f, beta = 0.0f; + + // Reinterpret float* as tfloat32_t* — same bit width, bitcast is intentional. + auto* tf32_A = reinterpret_cast(d_A); + auto* tf32_B = reinterpret_cast(d_B); + + CutlassGemm gemm_op; + CutlassGemm::Arguments args( + {S, S, S}, // problem size M, N, K + {tf32_A, S}, // A ref: pointer + leading dim + {tf32_B, S}, // B ref + {d_C, S}, // C ref (source for beta*C) + {d_C, S}, // D ref (output destination) + {alpha, beta} // epilogue: alpha*A*B + beta*C + ); + + // Allocate workspace (zero for split_k_slices=1) + size_t workspace_bytes = gemm_op.get_workspace_size(args); + void* workspace = nullptr; + if (workspace_bytes) + CUDA_CHECK(cudaMalloc(&workspace, workspace_bytes)); + + for (int i = 0; i < WARMUP; ++i) { + cutlass::Status s = gemm_op(args, workspace); + if (s != cutlass::Status::kSuccess) { + fprintf(stderr, "CUTLASS error %d\n", (int)s); + exit(EXIT_FAILURE); + } + } + CUDA_CHECK(cudaDeviceSynchronize()); + + GpuTimer timer; + float total_ms = 0.0f; + for (int i = 0; i < ITERS; ++i) { + CUDA_CHECK(cudaMemset(d_C, 0, sz)); + timer.start(); + gemm_op(args, workspace); + total_ms += timer.stop(); + } + + if (workspace) CUDA_CHECK(cudaFree(workspace)); + CUDA_CHECK(cudaFree(d_A)); + CUDA_CHECK(cudaFree(d_B)); + CUDA_CHECK(cudaFree(d_C)); + + return total_ms / ITERS; +} diff --git a/cuda-mma/src/benchmarks.cu b/cuda-mma/src/benchmarks.cu index 688fac0..8c2b6fc 100644 --- a/cuda-mma/src/benchmarks.cu +++ b/cuda-mma/src/benchmarks.cu @@ -136,6 +136,18 @@ float benchmark_coalesced(int S, const float* h_A, const float* h_B) return total_ms / ITERS; } +// Persistent handle for cuBLAS +static cublasHandle_t cublas_handle() +{ + static cublasHandle_t h = []() { + cublasHandle_t handle; + CUBLAS_CHECK(cublasCreate(&handle)); + CUBLAS_CHECK(cublasSetMathMode(handle, CUBLAS_TF32_TENSOR_OP_MATH)); + return handle; + }(); + return h; +} + float benchmark_cublas(int S, const float* h_A, const float* h_B) { size_t szA = (size_t)S * S * sizeof(float); @@ -151,8 +163,7 @@ float benchmark_cublas(int S, const float* h_A, const float* h_B) CUDA_CHECK(cudaMemcpy(d_B, h_B, szB, cudaMemcpyHostToDevice)); CUDA_CHECK(cudaMemset(d_C, 0, szC)); - cublasHandle_t handle; - CUBLAS_CHECK(cublasCreate(&handle)); + cublasHandle_t handle = cublas_handle(); // cuBLAS is column-major. For row-major C=A*B, use the identity // C^T = B^T * A^T, so pass B first with leading dimension S. @@ -185,7 +196,6 @@ float benchmark_cublas(int S, const float* h_A, const float* h_B) total_ms += timer.stop(); } - CUBLAS_CHECK(cublasDestroy(handle)); CUDA_CHECK(cudaFree(d_A)); CUDA_CHECK(cudaFree(d_B)); CUDA_CHECK(cudaFree(d_C)); diff --git a/cuda-mma/src/kernels.cu b/cuda-mma/src/kernels.cu index 0d6709a..80eaa0a 100644 --- a/cuda-mma/src/kernels.cu +++ b/cuda-mma/src/kernels.cu @@ -98,83 +98,6 @@ __global__ void sgemm_tiled( C[row * N + col] = alpha * acc + beta * C[row * N + col]; } - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -// // Tiled version: 1 thread per output element, shared memory for tiles of A and B. -// template -// __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) -// { -// __shared__ float sA[TILE_SIZE][TILE_SIZE]; -// __shared__ float sB[TILE_SIZE][TILE_SIZE]; - -// int row = blockIdx.y * TILE_SIZE + threadIdx.y; -// int col = blockIdx.x * TILE_SIZE + threadIdx.x; - -// float acc = 0.0f; - -// // Iterate over tiles along the K dimension. -// for (size_t t = 0; t < (K + TILE_SIZE - 1) / TILE_SIZE; ++t) { -// size_t aCol = t * TILE_SIZE + threadIdx.x; // column of A this thread loads -// size_t bRow = t * TILE_SIZE + threadIdx.y; // row of B this thread loads - -// // Boundary-safe loads: pad with 0 for out-of-bounds tiles. -// sA[threadIdx.y][threadIdx.x] = (row < M && aCol < K) ? A[row * K + aCol] : 0.0f; -// sB[threadIdx.y][threadIdx.x] = (bRow < K && col < N) ? B[bRow * N + col] : 0.0f; -// __syncthreads(); - -// // Accumulate partial dot product from shared memory. -// #pragma unroll -// for (int i = 0; i < TILE_SIZE; ++i) -// acc += sA[threadIdx.y][i] * sB[i][threadIdx.x]; -// __syncthreads(); -// } - -// if (row < M && col < N) -// C[row * N + col] = alpha * acc + beta * C[row * N + col]; -// } - template __global__ void sgemm_tiled<16>(size_t M, size_t N, size_t K, float alpha, const float *A, const float *B, float beta, float *C); From ab453c0ad6ed32b82a08e1fd7edd0748fe74d141 Mon Sep 17 00:00:00 2001 From: Sebas <12949040+SebastiaAgramunt@users.noreply.github.com> Date: Sat, 20 Jun 2026 06:11:51 +0000 Subject: [PATCH 11/21] better plot --- cuda-mma/scripts/plot_results.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cuda-mma/scripts/plot_results.py b/cuda-mma/scripts/plot_results.py index ccff442..dd34e08 100644 --- a/cuda-mma/scripts/plot_results.py +++ b/cuda-mma/scripts/plot_results.py @@ -32,7 +32,7 @@ def load(filename): return None return pd.read_csv(path) -fig, axes = plt.subplots(1, 2, figsize=(12, 5)) +fig, axes = plt.subplots(2, 1, figsize=(8, 12)) fig.suptitle("CUDA SGEMM Roofline Analysis - Nvidia GPU A100") ax_gflops, ax_bw = axes From 8733821980b6406703c51ec9ab5d7330f2d31e31 Mon Sep 17 00:00:00 2001 From: Sebas <12949040+SebastiaAgramunt@users.noreply.github.com> Date: Sat, 20 Jun 2026 21:47:06 +0000 Subject: [PATCH 12/21] improvements --- cuda-mma/scripts/plot_results.py | 82 +++++++++++++++++++------------- cuda-mma/src/kernels.cu | 41 ++++++++++------ 2 files changed, 76 insertions(+), 47 deletions(-) diff --git a/cuda-mma/scripts/plot_results.py b/cuda-mma/scripts/plot_results.py index dd34e08..1adadcb 100644 --- a/cuda-mma/scripts/plot_results.py +++ b/cuda-mma/scripts/plot_results.py @@ -76,38 +76,54 @@ def ai_fmt(x, _): print(f"Saved {out_path}") plt.show() -# --- Naive vs cBLAS throughput comparison --- -naive_vs_cblas = { - "Naive": "naive.csv", - "cBLAS": "cblas.csv", -} +def plot_comparison(series, title, out_name): + fig, ax = plt.subplots(figsize=(7, 5)) + fig.suptitle(title) -fig2, ax2 = plt.subplots(figsize=(7, 5)) -fig2.suptitle("SGEMM Roofline Analysis - Nvidia GPU A100 vs AMD EPYC 7J13 64-Core") + for label, filename in series.items(): + df = load(filename) + if df is None: + print(f"Warning: {filename} not found, skipping.") + continue + ax.plot(df["size"], df["gflops"], marker="o", label=label, color=colors[label]) -for label, filename in naive_vs_cblas.items(): - df = load(filename) - if df is None: - print(f"Warning: {filename} not found, skipping.") - continue - ax2.plot(df["size"], df["gflops"], marker="o", label=label, color=colors[label]) - -ax2.set_xscale("log", base=2) -ax2.set_yscale("log") -ax2.xaxis.set_major_formatter(ticker.FuncFormatter(lambda x, _: f"{int(x)}")) -ax2.set_xlabel("Matrix size (S×S)") -ax2.set_ylabel("GFLOP/s") -ax2.set_title("Throughput") -ax2.legend() -ax2.grid(True, which="both", linestyle="--", linewidth=0.5) - -sec2 = ax2.secondary_xaxis("top", functions=(lambda x: x / 8 * 1e6, lambda x: x * 8 / 1e6)) -sec2.set_xscale("log", base=2) -sec2.xaxis.set_major_formatter(ticker.FuncFormatter(ai_fmt)) -sec2.set_xlabel("Arithmetic Intensity (FLOP/MB)") - -fig2.tight_layout() -out_path2 = OUTPUT_DIR / "roofline_naive_cblas.png" -plt.savefig(out_path2, dpi=150) -print(f"Saved {out_path2}") -plt.show() + ax.set_xscale("log", base=2) + ax.set_yscale("log") + ax.xaxis.set_major_formatter(ticker.FuncFormatter(lambda x, _: f"{int(x)}")) + ax.set_xlabel("Matrix size (S×S)") + ax.set_ylabel("GFLOP/s") + ax.set_title("Throughput") + ax.legend() + ax.grid(True, which="both", linestyle="--", linewidth=0.5) + + sec = ax.secondary_xaxis("top", functions=(lambda x: x / 8 * 1e6, lambda x: x * 8 / 1e6)) + sec.set_xscale("log", base=2) + sec.xaxis.set_major_formatter(ticker.FuncFormatter(ai_fmt)) + sec.set_xlabel("Arithmetic Intensity (FLOP/MB)") + + fig.tight_layout() + out_path = OUTPUT_DIR / out_name + plt.savefig(out_path, dpi=150) + print(f"Saved {out_path}") + plt.show() + +# --- Naive vs cBLAS throughput comparison --- +plot_comparison( + {"Naive": "naive.csv", "cBLAS": "cblas.csv"}, + "SGEMM Roofline Analysis - Nvidia GPU A100 vs AMD EPYC 7J13 64-Core", + "roofline_naive_cblas.png", +) + +# --- Coalesced vs Tiled throughput comparison --- +plot_comparison( + {"Coalesced": "coalesced.csv", "Tiled": "tiled.csv"}, + "SGEMM Roofline Analysis - Nvidia GPU A100", + "roofline_coalesced_tiled.png", +) + +# --- Naive vs Coalesced throughput comparison --- +plot_comparison( + {"Naive": "naive.csv", "Coalesced": "coalesced.csv"}, + "SGEMM Roofline Analysis - Nvidia GPU A100", + "roofline_naive_coalesced.png", +) diff --git a/cuda-mma/src/kernels.cu b/cuda-mma/src/kernels.cu index 80eaa0a..bfd5836 100644 --- a/cuda-mma/src/kernels.cu +++ b/cuda-mma/src/kernels.cu @@ -17,16 +17,22 @@ __global__ void sgemm_naive( float *C) // [M x N] row-major (in-out: C = alpha*A*B + beta*C) { - int row = blockIdx.y * blockDim.y + threadIdx.y; // 0 .. M-1 - int col = blockIdx.x * blockDim.x + threadIdx.x; // 0 .. N-1 + const uint x = blockIdx.x * blockDim.x + threadIdx.x; // 0 .. M-1 (Rows of C) + const uint y = blockIdx.y * blockDim.y + threadIdx.y; // 0 .. N-1 (Columns of C) - if (row >= M || col >= N) return; + // Warp access e.g. warp 0: + // x is consecutive for threads 0..31, y is the same for threads 0..31 + // each warp accesses a column of C, with non-coalesced accesses to A and C, but broadcast access to B + // Ref: https://siboehm.com/articles/22/CUDA-MMM + + if (x >= M || y >= N) return; float acc = 0.0f; - for (size_t i = 0; i < K; ++i) - acc += A[row * K + i] * B[i * N + col]; + for (size_t i = 0; i < K; ++i){ + acc += A[x * K + i] * B[i * N + y]; + } - C[row * N + col] = alpha * acc + beta * C[row * N + col]; + C[x * N + y] = alpha * acc + beta * C[x * N + y]; } @@ -41,16 +47,23 @@ __global__ void sgemm_coalesced( const float *B, float beta, float *C) { - const int cRow = blockIdx.x * BLOCKSIZE + (threadIdx.x / BLOCKSIZE); - const int cCol = blockIdx.y * BLOCKSIZE + (threadIdx.x % BLOCKSIZE); - if (cRow < M && cCol < N) { - float tmp = 0.0; - for (size_t i = 0; i < K; ++i) { - tmp += A[cRow * K + i] * B[i * N + cCol]; + const uint x = blockIdx.x * BLOCKSIZE + (threadIdx.x / BLOCKSIZE); // 0 .. M-1 (Rows of C) + const uint y = blockIdx.y * BLOCKSIZE + (threadIdx.x % BLOCKSIZE); // 0 .. N-1 (Columns of C) + + // Warp access e.g. warp 0: + // x is the same for threads 0..31, y is consecutive for threads 0..31 + // each warp accesses a BLOCKSIZE (32) tile of C, with coalesced accesses to A and C but broadcast acces to B + // Ref: https://siboehm.com/articles/22/CUDA-MMM + + if (x >= M || y >= N) return; + + float acc = 0.0f; + for (size_t i = 0; i < K; ++i){ + acc += A[x * K + i] * B[i * N + y]; } - C[cRow * N + cCol] = alpha * tmp + beta * C[cRow * N + cCol]; - } + + C[x * N + y] = alpha * acc + beta * C[x * N + y]; } // Tiled version: 1 thread per output element, shared memory for tiles of A and B. From 9fced8cf610dbdaa59516dfc939d28acc9685383 Mon Sep 17 00:00:00 2001 From: Sebas <12949040+SebastiaAgramunt@users.noreply.github.com> Date: Sun, 21 Jun 2026 07:18:11 +0000 Subject: [PATCH 13/21] fixes --- cuda-mma/scripts/plot_results.py | 14 ++++++++++++++ cuda-mma/src/benchmark_cutlass.cu | 2 +- cuda-mma/src/kernels.cu | 12 ++++++------ 3 files changed, 21 insertions(+), 7 deletions(-) diff --git a/cuda-mma/scripts/plot_results.py b/cuda-mma/scripts/plot_results.py index 1adadcb..97dcfc7 100644 --- a/cuda-mma/scripts/plot_results.py +++ b/cuda-mma/scripts/plot_results.py @@ -127,3 +127,17 @@ def plot_comparison(series, title, out_name): "SGEMM Roofline Analysis - Nvidia GPU A100", "roofline_naive_coalesced.png", ) + +# --- Tiled vs cuBLAS throughput comparison --- +plot_comparison( + {"Tiled": "tiled.csv", "cuBLAS": "cublas.csv"}, + "SGEMM Roofline Analysis - Nvidia GPU A100", + "roofline_tiled_cublas.png", +) + +# --- cuBLAS vs CUTLASS throughput comparison --- +plot_comparison( + {"cuBLAS": "cublas.csv", "CUTLASS": "cutlass.csv"}, + "SGEMM Roofline Analysis - Nvidia GPU A100", + "roofline_cublas_cutlass.png", +) diff --git a/cuda-mma/src/benchmark_cutlass.cu b/cuda-mma/src/benchmark_cutlass.cu index 84ffb3f..2f8909d 100644 --- a/cuda-mma/src/benchmark_cutlass.cu +++ b/cuda-mma/src/benchmark_cutlass.cu @@ -18,7 +18,7 @@ using CutlassGemm = cutlass::gemm::device::Gemm< cutlass::tfloat32_t, cutlass::layout::RowMajor, // A (tf32 ~ float with 10-bit mantissa) cutlass::tfloat32_t, cutlass::layout::RowMajor, // B float, cutlass::layout::RowMajor, // C / D (full float output) - float, // accumulator + float, // accumulator cutlass::arch::OpClassTensorOp, // TF32 tensor cores cutlass::arch::Sm80, // Ampere cutlass::gemm::GemmShape<128, 128, 32>, // threadblock tile diff --git a/cuda-mma/src/kernels.cu b/cuda-mma/src/kernels.cu index bfd5836..d05e262 100644 --- a/cuda-mma/src/kernels.cu +++ b/cuda-mma/src/kernels.cu @@ -83,8 +83,8 @@ __global__ void sgemm_tiled( __shared__ float B_tile[TILE_SIZE * TILE_SIZE]; // the element of C for this specific thread - size_t row = blockIdx.y * TILE_SIZE + threadIdx.y; - size_t col = blockIdx.x * TILE_SIZE + threadIdx.x; + size_t y = blockIdx.y * TILE_SIZE + threadIdx.y; // Rows of C + size_t x = blockIdx.x * TILE_SIZE + threadIdx.x; // Columns of C // we multiply along the K dimension, iterate over this size_t n_tiles = (K + TILE_SIZE -1)/TILE_SIZE; @@ -96,8 +96,8 @@ __global__ void sgemm_tiled( size_t aCol = tile * TILE_SIZE + threadIdx.x; size_t bRow = tile * TILE_SIZE + threadIdx.y; - A_tile[threadIdx.y * TILE_SIZE + threadIdx.x] = (row < M && aCol < K) ? A[row * K + aCol] : 0.0f; - B_tile[threadIdx.y * TILE_SIZE + threadIdx.x] = (bRow < K && col < N) ? B[bRow * N + col] : 0.0f; + A_tile[threadIdx.y * TILE_SIZE + threadIdx.x] = (y < M && aCol < K) ? A[y * K + aCol] : 0.0f; + B_tile[threadIdx.y * TILE_SIZE + threadIdx.x] = (bRow < K && x < N) ? B[bRow * N + x] : 0.0f; __syncthreads(); @@ -107,8 +107,8 @@ __global__ void sgemm_tiled( __syncthreads(); } - if (row < M && col < N) - C[row * N + col] = alpha * acc + beta * C[row * N + col]; + if (y < M && x < N) + C[y * N + x] = alpha * acc + beta * C[y * N + x]; } template __global__ void sgemm_tiled<16>(size_t M, size_t N, size_t K, float alpha, From f3e4b6281e47503d4a7b4db8b52d5d1c63af05bb Mon Sep 17 00:00:00 2001 From: Sebas <12949040+SebastiaAgramunt@users.noreply.github.com> Date: Mon, 22 Jun 2026 00:55:28 +0000 Subject: [PATCH 14/21] adding plots --- cuda-mma/include/benchmarks.h | 3 +- cuda-mma/main.cu | 15 ++-- cuda-mma/scripts/plot_results.py | 4 +- cuda-mma/src/benchmark_cutlass.cu | 131 +++++++++++++++++++++++------- 4 files changed, 113 insertions(+), 40 deletions(-) diff --git a/cuda-mma/include/benchmarks.h b/cuda-mma/include/benchmarks.h index 329ab74..a0cfc7a 100644 --- a/cuda-mma/include/benchmarks.h +++ b/cuda-mma/include/benchmarks.h @@ -6,5 +6,6 @@ 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(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); diff --git a/cuda-mma/main.cu b/cuda-mma/main.cu index 185a997..38ea3ad 100644 --- a/cuda-mma/main.cu +++ b/cuda-mma/main.cu @@ -13,7 +13,9 @@ 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) + 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", @@ -29,10 +31,8 @@ static void run_sweep(const char* label, BenchmarkFn fn, float avg_ms = fn(S, h_A, h_B); double flops = 2.0 * S * S * S; - double bytes = ((double)S * S - + (double)S * S - + 2.0 * S * S) - * sizeof(float); + 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; @@ -65,8 +65,11 @@ int main() { 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", benchmark_cutlass, sizes, N_SIZES, h_A, h_B, "output/cutlass.csv"); + run_sweep("CUTLASS_fp32", benchmark_cutlass_fp32, sizes, N_SIZES, h_A, h_B, "output/cutlass_fp32.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]); diff --git a/cuda-mma/scripts/plot_results.py b/cuda-mma/scripts/plot_results.py index 97dcfc7..653e4e8 100644 --- a/cuda-mma/scripts/plot_results.py +++ b/cuda-mma/scripts/plot_results.py @@ -13,7 +13,7 @@ "Tiled": "tiled.csv", "Coarsened": "coarsened.csv", "cuBLAS": "cublas.csv", - "CUTLASS": "cutlass.csv", + "CUTLASS": "cutlass_fp32.csv", } colors = { @@ -137,7 +137,7 @@ def plot_comparison(series, title, out_name): # --- cuBLAS vs CUTLASS throughput comparison --- plot_comparison( - {"cuBLAS": "cublas.csv", "CUTLASS": "cutlass.csv"}, + {"cuBLAS": "cublas.csv", "CUTLASS": "cutlass_fp32.csv"}, "SGEMM Roofline Analysis - Nvidia GPU A100", "roofline_cublas_cutlass.png", ) diff --git a/cuda-mma/src/benchmark_cutlass.cu b/cuda-mma/src/benchmark_cutlass.cu index 2f8909d..7a4f409 100644 --- a/cuda-mma/src/benchmark_cutlass.cu +++ b/cuda-mma/src/benchmark_cutlass.cu @@ -1,7 +1,9 @@ #include #include +#include #include +#include #include "cuda_check.cuh" #include "timer.h" @@ -11,45 +13,47 @@ #define WARMUP 3 #define ITERS 5 -// SGEMM via CUTLASS using TF32 Tensor Cores on Ampere (sm_80). -// CUTLASS 4.x requires tfloat32_t element types to select TF32 MMA instructions; -// float inputs are bitcast to tf32 precision (rounds the mantissa to 10 bits). -using CutlassGemm = cutlass::gemm::device::Gemm< - cutlass::tfloat32_t, cutlass::layout::RowMajor, // A (tf32 ~ float with 10-bit mantissa) - cutlass::tfloat32_t, cutlass::layout::RowMajor, // B - float, cutlass::layout::RowMajor, // C / D (full float output) - float, // accumulator - cutlass::arch::OpClassTensorOp, // TF32 tensor cores - cutlass::arch::Sm80, // Ampere - cutlass::gemm::GemmShape<128, 128, 32>, // threadblock tile - cutlass::gemm::GemmShape<64, 64, 32>, // warp tile - cutlass::gemm::GemmShape<16, 8, 8> // TF32 MMA instruction shape ->; - -float benchmark_cutlass(int S, const float* h_A, const float* h_B) + +float benchmark_cutlass_fp16(int S, const float* h_A, const float* h_B) { - size_t sz = (size_t)S * S * sizeof(float); + // SGEMM via CUTLASS using FP16 Tensor Cores on Ampere (sm_80). + using CutlassGemm = cutlass::gemm::device::Gemm< + cutlass::half_t, cutlass::layout::RowMajor, // A (fp16 half precision) and row major + cutlass::half_t, cutlass::layout::RowMajor, // B (fp16 half precision) and row major + float, cutlass::layout::RowMajor, // C / D (full float output) + float, // accumulator + cutlass::arch::OpClassTensorOp, // FP16 tensor cores + cutlass::arch::Sm80 // Ampere architecture for our A100 GPU + >; - float *d_A, *d_B, *d_C; - CUDA_CHECK(cudaMalloc(&d_A, sz)); - CUDA_CHECK(cudaMalloc(&d_B, sz)); - CUDA_CHECK(cudaMalloc(&d_C, sz)); + size_t n = (size_t)S * S; + size_t szAB = n * sizeof(cutlass::half_t); + size_t szC = n * sizeof(float); - CUDA_CHECK(cudaMemcpy(d_A, h_A, sz, cudaMemcpyHostToDevice)); - CUDA_CHECK(cudaMemcpy(d_B, h_B, sz, cudaMemcpyHostToDevice)); - CUDA_CHECK(cudaMemset(d_C, 0, sz)); + // Narrow the float operands to half on the host (outside the timed region). + std::vector h_A_half(n), h_B_half(n); + for (size_t i = 0; i < n; ++i) { + h_A_half[i] = cutlass::half_t(h_A[i]); + h_B_half[i] = cutlass::half_t(h_B[i]); + } - const float alpha = 1.0f, beta = 0.0f; + cutlass::half_t *d_A, *d_B; + float *d_C; + CUDA_CHECK(cudaMalloc(&d_A, szAB)); + CUDA_CHECK(cudaMalloc(&d_B, szAB)); + CUDA_CHECK(cudaMalloc(&d_C, szC)); - // Reinterpret float* as tfloat32_t* — same bit width, bitcast is intentional. - auto* tf32_A = reinterpret_cast(d_A); - auto* tf32_B = reinterpret_cast(d_B); + CUDA_CHECK(cudaMemcpy(d_A, h_A_half.data(), szAB, cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemcpy(d_B, h_B_half.data(), szAB, cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemset(d_C, 0, szC)); + + const float alpha = 1.0f, beta = 0.0f; CutlassGemm gemm_op; CutlassGemm::Arguments args( {S, S, S}, // problem size M, N, K - {tf32_A, S}, // A ref: pointer + leading dim - {tf32_B, S}, // B ref + {d_A, S}, // A ref: pointer + leading dim + {d_B, S}, // B ref {d_C, S}, // C ref (source for beta*C) {d_C, S}, // D ref (output destination) {alpha, beta} // epilogue: alpha*A*B + beta*C @@ -73,7 +77,7 @@ float benchmark_cutlass(int S, const float* h_A, const float* h_B) GpuTimer timer; float total_ms = 0.0f; for (int i = 0; i < ITERS; ++i) { - CUDA_CHECK(cudaMemset(d_C, 0, sz)); + CUDA_CHECK(cudaMemset(d_C, 0, szC)); timer.start(); gemm_op(args, workspace); total_ms += timer.stop(); @@ -86,3 +90,68 @@ float benchmark_cutlass(int S, const float* h_A, const float* h_B) return total_ms / ITERS; } + +float benchmark_cutlass_fp32(int S, const float* h_A, const float* h_B) +{ + using CutlassGemmFP32 = cutlass::gemm::device::Gemm< + float, cutlass::layout::RowMajor, + float, cutlass::layout::RowMajor, + float, cutlass::layout::RowMajor, + float, + cutlass::arch::OpClassSimt, // CUDA cores, not tensor cores + cutlass::arch::Sm80 + >; + size_t n = (size_t)S * S; + size_t sz = n * sizeof(float); + + // No precision narrowing needed — upload float directly. + float *d_A, *d_B, *d_C; + CUDA_CHECK(cudaMalloc(&d_A, sz)); + CUDA_CHECK(cudaMalloc(&d_B, sz)); + CUDA_CHECK(cudaMalloc(&d_C, sz)); + CUDA_CHECK(cudaMemcpy(d_A, h_A, sz, cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemcpy(d_B, h_B, sz, cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemset(d_C, 0, sz)); + + const float alpha = 1.0f, beta = 0.0f; + + CutlassGemmFP32 gemm_op; + CutlassGemmFP32::Arguments args( + {S, S, S}, + {d_A, S}, + {d_B, S}, + {d_C, S}, + {d_C, S}, + {alpha, beta} + ); + + size_t workspace_bytes = gemm_op.get_workspace_size(args); + void* workspace = nullptr; + if (workspace_bytes) + CUDA_CHECK(cudaMalloc(&workspace, workspace_bytes)); + + for (int i = 0; i < WARMUP; ++i) { + cutlass::Status s = gemm_op(args, workspace); + if (s != cutlass::Status::kSuccess) { + fprintf(stderr, "CUTLASS error %d\n", (int)s); + exit(EXIT_FAILURE); + } + } + CUDA_CHECK(cudaDeviceSynchronize()); + + GpuTimer timer; + float total_ms = 0.0f; + for (int i = 0; i < ITERS; ++i) { + CUDA_CHECK(cudaMemset(d_C, 0, sz)); + timer.start(); + gemm_op(args, workspace); + total_ms += timer.stop(); + } + + if (workspace) CUDA_CHECK(cudaFree(workspace)); + CUDA_CHECK(cudaFree(d_A)); + CUDA_CHECK(cudaFree(d_B)); + CUDA_CHECK(cudaFree(d_C)); + + return total_ms / ITERS; +} \ No newline at end of file From d1f7226d43afa7318747419e92026df1dba6f1d8 Mon Sep 17 00:00:00 2001 From: Sebas <12949040+SebastiaAgramunt@users.noreply.github.com> Date: Mon, 22 Jun 2026 06:58:48 +0000 Subject: [PATCH 15/21] finally --- cuda-mma/include/benchmarks.h | 1 + cuda-mma/main.cu | 3 + cuda-mma/scripts/plot_results.py | 32 +++++++++-- cuda-mma/src/benchmark_cutlass.cu | 95 ++++++++++++++++++++++++++++--- 4 files changed, 118 insertions(+), 13 deletions(-) diff --git a/cuda-mma/include/benchmarks.h b/cuda-mma/include/benchmarks.h index a0cfc7a..2b5d26f 100644 --- a/cuda-mma/include/benchmarks.h +++ b/cuda-mma/include/benchmarks.h @@ -8,4 +8,5 @@ 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); diff --git a/cuda-mma/main.cu b/cuda-mma/main.cu index 38ea3ad..0946738 100644 --- a/cuda-mma/main.cu +++ b/cuda-mma/main.cu @@ -67,6 +67,9 @@ int main() { 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)); diff --git a/cuda-mma/scripts/plot_results.py b/cuda-mma/scripts/plot_results.py index 653e4e8..f767c85 100644 --- a/cuda-mma/scripts/plot_results.py +++ b/cuda-mma/scripts/plot_results.py @@ -11,19 +11,19 @@ "Naive": "naive.csv", "Coalesced": "coalesced.csv", "Tiled": "tiled.csv", - "Coarsened": "coarsened.csv", "cuBLAS": "cublas.csv", - "CUTLASS": "cutlass_fp32.csv", + "CUTLASS (TF32 TensorOp)": "cutlass_tf32.csv", } colors = { "Naive": "C0", "Coalesced": "C1", "Tiled": "C2", - "Coarsened": "C3", "cuBLAS": "C4", "cBLAS": "C6", - "CUTLASS": "C5", + "CUTLASS (FP32 SIMT)": "C5", + "CUTLASS (TF32 TensorOp)": "C7", + "CUTLASS (FP16 TensorOp)": "C8", } def load(filename): @@ -76,7 +76,7 @@ def ai_fmt(x, _): print(f"Saved {out_path}") plt.show() -def plot_comparison(series, title, out_name): +def plot_comparison(series, title, out_name, note=None): fig, ax = plt.subplots(figsize=(7, 5)) fig.suptitle(title) @@ -102,6 +102,12 @@ def plot_comparison(series, title, out_name): sec.set_xlabel("Arithmetic Intensity (FLOP/MB)") fig.tight_layout() + + if note: + fig.subplots_adjust(bottom=0.18) + fig.text(0.5, 0.02, note, ha="center", va="bottom", fontsize=8, + style="italic", wrap=True) + out_path = OUTPUT_DIR / out_name plt.savefig(out_path, dpi=150) print(f"Saved {out_path}") @@ -137,7 +143,21 @@ def plot_comparison(series, title, out_name): # --- cuBLAS vs CUTLASS throughput comparison --- plot_comparison( - {"cuBLAS": "cublas.csv", "CUTLASS": "cutlass_fp32.csv"}, + {"cuBLAS": "cublas.csv", "CUTLASS (FP32 SIMT)": "cutlass_fp32.csv"}, "SGEMM Roofline Analysis - Nvidia GPU A100", "roofline_cublas_cutlass.png", ) + +# --- CUTLASS FP32 vs TF32 vs FP16 vs cuBLAS throughput comparison --- +plot_comparison( + { + "CUTLASS (FP32 SIMT)": "cutlass_fp32.csv", + "CUTLASS (TF32 TensorOp)": "cutlass_tf32.csv", + "CUTLASS (FP16 TensorOp)": "cutlass_fp16.csv", + "cuBLAS": "cublas.csv", + }, + "SGEMM Roofline Analysis - Nvidia GPU A100", + "roofline_cutlass_variants.png", + note="Note: top axis assumes 4-byte operands (AI = S/8); FP16 TensorOp uses " + "fp16 A/B and its true AI is ~1.33x higher than shown here.", +) diff --git a/cuda-mma/src/benchmark_cutlass.cu b/cuda-mma/src/benchmark_cutlass.cu index 7a4f409..34b7aa3 100644 --- a/cuda-mma/src/benchmark_cutlass.cu +++ b/cuda-mma/src/benchmark_cutlass.cu @@ -13,7 +13,7 @@ #define WARMUP 3 #define ITERS 5 - +// CUTLASS version using FP16 Tensor Cores; A/B are narrowed to half precision on host before upload. Fast, but lower precision and requires extra host-side copy to narrow operands. float benchmark_cutlass_fp16(int S, const float* h_A, const float* h_B) { // SGEMM via CUTLASS using FP16 Tensor Cores on Ampere (sm_80). @@ -91,15 +91,96 @@ float benchmark_cutlass_fp16(int S, const float* h_A, const float* h_B) return total_ms / ITERS; } +// CUTLASS version using TF32 Tensor Cores; A/B are bitcast (still 4 bytes), not narrowed. Fast +float benchmark_cutlass_tf32(int S, const float* h_A, const float* h_B) +{ + // SGEMM via CUTLASS using TF32 Tensor Cores on Ampere (sm_80). + // tfloat32_t is bit-identical width to float (4 bytes); float inputs are + // bitcast, not copied, which truncates the mantissa to 10 bits for the MMA. + // Explicit tile shapes are required: CUTLASS's generic Sm80/tf32 default + // config hits a broken epilogue dispatch (FragmentIteratorComplexTensorOp) + // in this CUTLASS version. + using CutlassGemmTF32 = cutlass::gemm::device::Gemm< + cutlass::tfloat32_t, cutlass::layout::RowMajor, // A (tf32 ~ float with 10-bit mantissa) + cutlass::tfloat32_t, cutlass::layout::RowMajor, // B + float, cutlass::layout::RowMajor, // C / D (full float output) + float, // accumulator + cutlass::arch::OpClassTensorOp, // TF32 tensor cores + cutlass::arch::Sm80, // Ampere + cutlass::gemm::GemmShape<128, 128, 32>, // threadblock tile + cutlass::gemm::GemmShape<64, 64, 32>, // warp tile + cutlass::gemm::GemmShape<16, 8, 8> // TF32 MMA instruction shape + >; + + size_t n = (size_t)S * S; + size_t sz = n * sizeof(float); + + float *d_A, *d_B, *d_C; + CUDA_CHECK(cudaMalloc(&d_A, sz)); + CUDA_CHECK(cudaMalloc(&d_B, sz)); + CUDA_CHECK(cudaMalloc(&d_C, sz)); + CUDA_CHECK(cudaMemcpy(d_A, h_A, sz, cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemcpy(d_B, h_B, sz, cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemset(d_C, 0, sz)); + + const float alpha = 1.0f, beta = 0.0f; + + // Reinterpret float* as tfloat32_t* — same bit width, bitcast is intentional. + auto* tf32_A = reinterpret_cast(d_A); + auto* tf32_B = reinterpret_cast(d_B); + + CutlassGemmTF32 gemm_op; + CutlassGemmTF32::Arguments args( + {S, S, S}, // problem size M, N, K + {tf32_A, S}, // A ref: pointer + leading dim + {tf32_B, S}, // B ref + {d_C, S}, // C ref (source for beta*C) + {d_C, S}, // D ref (output destination) + {alpha, beta} // epilogue: alpha*A*B + beta*C + ); + + // Allocate workspace (zero for split_k_slices=1) + size_t workspace_bytes = gemm_op.get_workspace_size(args); + void* workspace = nullptr; + if (workspace_bytes) + CUDA_CHECK(cudaMalloc(&workspace, workspace_bytes)); + + for (int i = 0; i < WARMUP; ++i) { + cutlass::Status s = gemm_op(args, workspace); + if (s != cutlass::Status::kSuccess) { + fprintf(stderr, "CUTLASS error %d\n", (int)s); + exit(EXIT_FAILURE); + } + } + CUDA_CHECK(cudaDeviceSynchronize()); + + GpuTimer timer; + float total_ms = 0.0f; + for (int i = 0; i < ITERS; ++i) { + CUDA_CHECK(cudaMemset(d_C, 0, sz)); + timer.start(); + gemm_op(args, workspace); + total_ms += timer.stop(); + } + + if (workspace) CUDA_CHECK(cudaFree(workspace)); + CUDA_CHECK(cudaFree(d_A)); + CUDA_CHECK(cudaFree(d_B)); + CUDA_CHECK(cudaFree(d_C)); + + return total_ms / ITERS; +} + +// CUTLASS version using regular fp32 CUDA cores (no precision narrowing, no tensor cores). Slow float benchmark_cutlass_fp32(int S, const float* h_A, const float* h_B) { using CutlassGemmFP32 = cutlass::gemm::device::Gemm< - float, cutlass::layout::RowMajor, - float, cutlass::layout::RowMajor, - float, cutlass::layout::RowMajor, - float, - cutlass::arch::OpClassSimt, // CUDA cores, not tensor cores - cutlass::arch::Sm80 + float, cutlass::layout::RowMajor, // A and row major + float, cutlass::layout::RowMajor, // B and row major + float, cutlass::layout::RowMajor, // C / D (full float output) + float, // accumulator + cutlass::arch::OpClassSimt, // CUDA cores, not tensor cores + cutlass::arch::Sm80 // Ampere architecture for our A100 GPU >; size_t n = (size_t)S * S; size_t sz = n * sizeof(float); From 54c0281474708bbe6cc2d7a7e6ad93e06abdf8fc Mon Sep 17 00:00:00 2001 From: Sebas <12949040+SebastiaAgramunt@users.noreply.github.com> Date: Mon, 22 Jun 2026 07:02:14 +0000 Subject: [PATCH 16/21] readme --- README.md | 1 + 1 file changed, 1 insertion(+) diff --git a/README.md b/README.md index c9b3ae2..d4210cb 100644 --- a/README.md +++ b/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) From 1c764882dab26319883608a34a65043cf6164386 Mon Sep 17 00:00:00 2001 From: Sebas <12949040+SebastiaAgramunt@users.noreply.github.com> Date: Thu, 9 Jul 2026 14:40:21 +0000 Subject: [PATCH 17/21] Add vLLM benchmark scripts Profiling install/benchmark helpers for vLLM llama3-8b; nsys-rep results and .deb installers stay untracked via .gitignore. Co-Authored-By: Claude Sonnet 5 --- .gitignore | 4 +++- vLLM-profiling/bench.py | 13 +++++++++++++ vLLM-profiling/benchmark.sh | 19 +++++++++++++++++++ vLLM-profiling/install.sh | 30 ++++++++++++++++++++++++++++++ 4 files changed, 65 insertions(+), 1 deletion(-) create mode 100644 vLLM-profiling/bench.py create mode 100755 vLLM-profiling/benchmark.sh create mode 100755 vLLM-profiling/install.sh diff --git a/.gitignore b/.gitignore index 21ccaa1..0d66ca8 100644 --- a/.gitignore +++ b/.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 \ No newline at end of file +*.deb diff --git a/vLLM-profiling/bench.py b/vLLM-profiling/bench.py new file mode 100644 index 0000000..26e369b --- /dev/null +++ b/vLLM-profiling/bench.py @@ -0,0 +1,13 @@ +# bench.py +from vllm import LLM, SamplingParams +import json + +llm = LLM(model="meta-llama/Meta-Llama-3-8B-Instruct", gpu_memory_utilization=0.85) + +prompts = [f"Explain the concept of {topic} in three sentences." + for topic in ["entropy", "recursion", "inflation", "photosynthesis"] * 8] # 32 prompts + +sampling = SamplingParams(temperature=0.7, max_tokens=128) + +outputs = llm.generate(prompts, sampling) +print(f"Generated {len(outputs)} completions") \ No newline at end of file diff --git a/vLLM-profiling/benchmark.sh b/vLLM-profiling/benchmark.sh new file mode 100755 index 0000000..8feb026 --- /dev/null +++ b/vLLM-profiling/benchmark.sh @@ -0,0 +1,19 @@ +#!/bin/bash + +THIS_DIR=$(dirname "$(realpath "$0")") +NSIGHT_SYS_VERSION=2026.3.1.157-3804839 +CAPPED_NSIGHT_SYS_VERSION=$(echo ${NSIGHT_SYS_VERSION} | cut -d. -f1-3) +NSIGHT_SYSTEMS_CLI_PATH=/opt/nvidia/nsight-systems-cli + +source ${THIS_DIR}/.env +source ${THIS_DIR}/.venv/bin/activate + +python ${THIS_DIR}/bench.py + +# now the real profiled run +${NSIGHT_SYSTEMS_CLI_PATH}/${CAPPED_NSIGHT_SYS_VERSION}/bin/nsys profile \ + --trace=cuda,nvtx,osrt \ + --cuda-graph-trace=node \ + --output=vllm_llama3_8b_bs32 \ + --force-overwrite=true \ + python ${THIS_DIR}/bench.py \ No newline at end of file diff --git a/vLLM-profiling/install.sh b/vLLM-profiling/install.sh new file mode 100755 index 0000000..2d0c72d --- /dev/null +++ b/vLLM-profiling/install.sh @@ -0,0 +1,30 @@ +#!/bin/bash + +THIS_DIR=$(dirname "$(realpath "$0")") +NSIGHT_SYS_VERSION=2026.3.1.157-3804839 +CAPPED_NSIGHT_SYS_VERSION=$(echo ${NSIGHT_SYS_VERSION} | cut -d. -f1-3) +NSIGHT_SYSTEMS_CLI_PATH=/opt/nvidia/nsight-systems-cli + +# create python virtual env +rm -rf ${THIS_DIR}/.venv +python3 -m venv ${THIS_DIR}/.venv + +# vLLM (pulls a matching torch/cuda build) +${THIS_DIR}/.venv/bin/python -m pip install --upgrade pip +${THIS_DIR}/.venv/bin/python -m pip install vllm huggingface_hub + +# Confirm the driver itself is healthy before touching packages +nvidia-smi || { echo "GPU/driver not responding — stop here, don't proceed"; exit 1; } + +wget https://developer.nvidia.com/downloads/assets/tools/secure/nsight-systems/2026_3/NsightSystems-linux-cli-public-${NSIGHT_SYS_VERSION}.deb -O ${THIS_DIR}/nsight-systems.deb +sudo apt install ${THIS_DIR}/nsight-systems.deb + +# # What did dpkg actually complain about? (rerun to see the real error, not swallowed by the script) +sudo dpkg -i ${THIS_DIR}/nsight-systems.deb + +# # If that reports missing deps, resolve them narrowly: +sudo apt-get install -f -y --no-install-recommends + +${NSIGHT_SYSTEMS_CLI_PATH}/${CAPPED_NSIGHT_SYS_VERSION}/bin/nsys --version + + From 996f53c4b5e7922522019ddf0250eb8598e9d1e2 Mon Sep 17 00:00:00 2001 From: Sebas <12949040+SebastiaAgramunt@users.noreply.github.com> Date: Thu, 9 Jul 2026 14:41:20 +0000 Subject: [PATCH 18/21] Add cuda-profiling: kernel and CUTLASS profiling harness vendor/, build/, and reports/ are excluded via cuda-profiling/.gitignore; run scripts/download_vendor.sh to fetch CUTLASS locally. Co-Authored-By: Claude Sonnet 5 --- cuda-profiling/.gitignore | 3 + cuda-profiling/Makefile | 102 ++++++++++ cuda-profiling/README.md | 70 +++++++ cuda-profiling/include/cuda_check.cuh | 40 ++++ cuda-profiling/include/kernels.cuh | 52 +++++ cuda-profiling/include/profiling.h | 31 +++ cuda-profiling/include/timer.h | 29 +++ cuda-profiling/include/utils.h | 6 + cuda-profiling/main.cu | 51 +++++ cuda-profiling/scripts/download_vendor.sh | 28 +++ cuda-profiling/src/kernels.cu | 193 +++++++++++++++++++ cuda-profiling/src/profile_cutlass.cu | 219 +++++++++++++++++++++ cuda-profiling/src/profile_kernels.cu | 224 ++++++++++++++++++++++ cuda-profiling/src/utils.cu | 6 + 14 files changed, 1054 insertions(+) create mode 100644 cuda-profiling/.gitignore create mode 100644 cuda-profiling/Makefile create mode 100644 cuda-profiling/README.md create mode 100644 cuda-profiling/include/cuda_check.cuh create mode 100644 cuda-profiling/include/kernels.cuh create mode 100644 cuda-profiling/include/profiling.h create mode 100644 cuda-profiling/include/timer.h create mode 100644 cuda-profiling/include/utils.h create mode 100644 cuda-profiling/main.cu create mode 100755 cuda-profiling/scripts/download_vendor.sh create mode 100644 cuda-profiling/src/kernels.cu create mode 100644 cuda-profiling/src/profile_cutlass.cu create mode 100644 cuda-profiling/src/profile_kernels.cu create mode 100644 cuda-profiling/src/utils.cu diff --git a/cuda-profiling/.gitignore b/cuda-profiling/.gitignore new file mode 100644 index 0000000..1c330bd --- /dev/null +++ b/cuda-profiling/.gitignore @@ -0,0 +1,3 @@ +build/ +reports/ +vendor/ diff --git a/cuda-profiling/Makefile b/cuda-profiling/Makefile new file mode 100644 index 0000000..d5b7dc2 --- /dev/null +++ b/cuda-profiling/Makefile @@ -0,0 +1,102 @@ +# --------------------------------------------------------------------------- +# Build + profiling harness for the cuda-mma kernels. +# --------------------------------------------------------------------------- + +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 + +# --generate-line-info correlates ncu's source/SASS views back to this code. +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_profiling +SRCS := main.cu src/kernels.cu src/utils.cu src/profile_kernels.cu src/profile_cutlass.cu +HDRS := include/cuda_check.cuh include/kernels.cuh include/timer.h \ + include/utils.h include/profiling.h + +REPORT_DIR := reports + +# Square matrix size used for profiling runs. +SIZE ?= 4096 +# One name from PROFILE_TABLE (see include/profiling.h): naive, coalesced, +# tiled, coarsened, cublas, cutlass_fp32, cutlass_tf32, cutlass_fp16. +# Required for `make ncu`; if empty, `make run`/`make nsys` profile all of them. +KERNEL ?= + +.PHONY: all clean run vendor ncu ncu-metrics nsys nsys-stats + +all: $(TARGET) + +$(BIN_DIR) $(REPORT_DIR): + mkdir -p $@ + +$(TARGET): $(SRCS) $(HDRS) | $(BIN_DIR) + $(NVCC) $(NVCC_FLAGS) -o $@ $(SRCS) -lcublas -ldl \ + -Xlinker -rpath,/usr/local/cuda/targets/x86_64-linux/lib + +vendor: + ./scripts/download_vendor.sh + +run: $(TARGET) + ./$(TARGET) $(SIZE) $(KERNEL) + +# --- Nsight Compute -------------------------------------------------------- +# Deep-dive a single kernel (occupancy, roofline, memory workload, etc). +# The NVTX range and the binary's own KERNEL filter both narrow execution to +# just that kernel, so `--set full`'s replay passes stay fast. +# make ncu KERNEL=tiled SIZE=2048 +ncu: $(TARGET) | $(REPORT_DIR) + @if [ -z "$(KERNEL)" ]; then \ + echo "Usage: make ncu KERNEL= [SIZE=4096]"; \ + echo "Available: naive coalesced tiled coarsened cublas cutlass_fp32 cutlass_tf32 cutlass_fp16"; \ + exit 1; \ + fi + ncu --set full \ + --nvtx --nvtx-include "$(KERNEL)/" \ + --launch-count 1 \ + -f -o $(REPORT_DIR)/ncu_$(KERNEL)_$(SIZE) \ + ./$(TARGET) $(SIZE) $(KERNEL) + +# Fast, fixed metric set across every kernel in one pass (CSV on stdout). +# Same metrics as cuda-mma's `make profile`, plus an SM throughput summary. +ncu-metrics: $(TARGET) | $(REPORT_DIR) + 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,\ + sm__throughput.avg.pct_of_peak_sustained_elapsed,\ + smsp__sass_thread_inst_executed_op_ffma_pred_on.sum \ + --nvtx --csv \ + ./$(TARGET) $(SIZE) $(KERNEL) > $(REPORT_DIR)/ncu_metrics_$(SIZE).csv + +# --- Nsight Systems --------------------------------------------------------- +# Whole-application timeline: every kernel launch, memcpy and NVTX range. +# Open reports/nsys_.nsys-rep in the Nsight Systems UI to inspect it. +nsys: $(TARGET) | $(REPORT_DIR) + nsys profile \ + --trace=cuda,nvtx,osrt \ + --force-overwrite=true \ + -o $(REPORT_DIR)/nsys_$(SIZE) \ + ./$(TARGET) $(SIZE) $(KERNEL) + +# Per-kernel duration summary extracted from the most recent `make nsys` run. +nsys-stats: | $(REPORT_DIR) + nsys stats --report cuda_gpu_kern_sum --format csv \ + --output $(REPORT_DIR)/nsys_$(SIZE)_kernels \ + $(REPORT_DIR)/nsys_$(SIZE).nsys-rep + +clean: + rm -rf build reports diff --git a/cuda-profiling/README.md b/cuda-profiling/README.md new file mode 100644 index 0000000..5ac70a5 --- /dev/null +++ b/cuda-profiling/README.md @@ -0,0 +1,70 @@ +# Profiling the cuda-mma SGEMM kernels with Nsight Compute and Nsight Systems + +Same SGEMM kernels as [cuda-mma](../cuda-mma) (naive, coalesced, tiled, +coarsened, cuBLAS, CUTLASS fp32/tf32/fp16), restructured for profiling +instead of roofline benchmarking: one binary, one NVTX range per kernel, +no CSV sweep. + +## Build + +```bash +./scripts/download_vendor.sh # fetches CUTLASS headers into vendor/ +make # builds build/bin/cuda_profiling +``` + +## Run directly + +```bash +./build/bin/cuda_profiling [SIZE] [KERNEL] +``` + +`SIZE` is the square matrix dimension (default 4096). `KERNEL` restricts +the run to one implementation — `naive`, `coalesced`, `tiled`, +`coarsened`, `cublas`, `cutlass_fp32`, `cutlass_tf32`, or `cutlass_fp16` +(see `include/profiling.h`). Omit it to run all of them in sequence. + +Every kernel's profiled launches are wrapped in an `NvtxRange` named after +it (`include/profiling.h`). That serves two purposes: Nsight Systems shows +it as a labeled block on the timeline, and Nsight Compute can be told to +collect metrics only inside that block with `--nvtx-include`. + +## Nsight Systems — whole-application timeline + +```bash +make nsys SIZE=4096 # -> reports/nsys_4096.nsys-rep +make nsys-stats SIZE=4096 # -> reports/nsys_4096_kernels.csv +``` + +Open the `.nsys-rep` in the Nsight Systems UI to see every kernel launch, +memcpy, and NVTX range on one timeline — useful for spotting gaps between +launches, H2D/D2H transfer overlap, and the relative wall-clock cost of +each implementation back-to-back. + +## Nsight Compute — per-kernel deep dive + +```bash +make ncu KERNEL=tiled SIZE=2048 # -> reports/ncu_tiled_2048.ncu-rep +``` + +This runs `--set full` scoped to the `tiled` NVTX range only +(`--nvtx-include "tiled/"`) and limited to the first launch +(`--launch-count 1`), so the replay passes a full section set requires +stay fast even though the binary itself launches every kernel after it +in the warmup loop. `KERNEL` is required for this target — `make ncu` +on its own prints the list of valid names. + +For a fast metrics-only sweep across every kernel in one pass (occupancy, +memory throughput) instead of one deep dive: + +```bash +make ncu-metrics SIZE=4096 # -> reports/ncu_metrics_4096.csv +``` + +## Notes + +- `--generate-line-info` is on by default in the `Makefile` so Nsight + Compute's source/SASS view can map hotspots back to this code. +- `KERNEL`/`SIZE` are plain Makefile variables — override on the command + line as shown above, or `export` them before calling `make run`. +- `make clean` removes `build/` and `reports/`; `vendor/` (CUTLASS) is + left alone since it's only fetched once. diff --git a/cuda-profiling/include/cuda_check.cuh b/cuda-profiling/include/cuda_check.cuh new file mode 100644 index 0000000..80f1368 --- /dev/null +++ b/cuda-profiling/include/cuda_check.cuh @@ -0,0 +1,40 @@ +// cuda_check.cuh +#pragma once + +#include +#include +#include + +// ── 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) diff --git a/cuda-profiling/include/kernels.cuh b/cuda-profiling/include/kernels.cuh new file mode 100644 index 0000000..75eeecc --- /dev/null +++ b/cuda-profiling/include/kernels.cuh @@ -0,0 +1,52 @@ +#pragma once + +#include + + +__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 +__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 +__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 +__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) diff --git a/cuda-profiling/include/profiling.h b/cuda-profiling/include/profiling.h new file mode 100644 index 0000000..da1e567 --- /dev/null +++ b/cuda-profiling/include/profiling.h @@ -0,0 +1,31 @@ +#pragma once + +#include + +// RAII NVTX range. Pushed/popped around the region you want visible as a +// named block in the Nsight Systems timeline, and selectable on its own +// with `ncu --nvtx-include "/"` (see Makefile targets `nsys`/`ncu`). +struct NvtxRange { + explicit NvtxRange(const char* name) { nvtxRangePushA(name); } + ~NvtxRange() { nvtxRangePop(); } +}; + +// One profiling entry point per kernel implementation. Each function mallocs +// its own device buffers, runs a few untimed warmup launches, then wraps the +// launch(es) to actually profile in an NvtxRange named after `name` below. +void profile_naive(int S, const float* h_A, const float* h_B); +void profile_coalesced(int S, const float* h_A, const float* h_B); +void profile_tiled(int S, const float* h_A, const float* h_B); +void profile_coarsened(int S, const float* h_A, const float* h_B); +void profile_cublas(int S, const float* h_A, const float* h_B); +void profile_cutlass_fp32(int S, const float* h_A, const float* h_B); +void profile_cutlass_tf32(int S, const float* h_A, const float* h_B); +void profile_cutlass_fp16(int S, const float* h_A, const float* h_B); + +struct ProfileEntry { + const char* name; // also the NVTX range label used by --nvtx-include + void (*fn)(int, const float*, const float*); +}; + +extern const ProfileEntry PROFILE_TABLE[]; +extern const int PROFILE_TABLE_SIZE; diff --git a/cuda-profiling/include/timer.h b/cuda-profiling/include/timer.h new file mode 100644 index 0000000..89e3a62 --- /dev/null +++ b/cuda-profiling/include/timer.h @@ -0,0 +1,29 @@ +#pragma once + +#include + +// --------------------------------------------------------------------------- +// 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; + } +}; diff --git a/cuda-profiling/include/utils.h b/cuda-profiling/include/utils.h new file mode 100644 index 0000000..e214a12 --- /dev/null +++ b/cuda-profiling/include/utils.h @@ -0,0 +1,6 @@ +#pragma once + +#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); diff --git a/cuda-profiling/main.cu b/cuda-profiling/main.cu new file mode 100644 index 0000000..5698404 --- /dev/null +++ b/cuda-profiling/main.cu @@ -0,0 +1,51 @@ +#include +#include +#include + +#include "profiling.h" +#include "utils.h" + +// Usage: cuda_profiling [SIZE] [KERNEL] +// SIZE square matrix dimension (default 4096) +// KERNEL one of the names in PROFILE_TABLE; if omitted, all run in sequence. +// +// Each kernel's profiled launches are wrapped in an NVTX range named after +// it, so a single run can be: +// - traced end-to-end with `nsys profile` (every kernel + memcpy on one +// timeline), or +// - narrowed to one kernel with `ncu --nvtx-include "/"`. +int main(int argc, char** argv) { + srand(42); + + int S = (argc > 1) ? atoi(argv[1]) : 4096; + const char* kernel = (argc > 2) ? argv[2] : nullptr; + + float* h_A = new float[(size_t)S * S]; + float* h_B = new float[(size_t)S * S]; + fill_random(h_A, S * S); + fill_random(h_B, S * S); + + bool ran_any = false; + for (int i = 0; i < PROFILE_TABLE_SIZE; ++i) { + const ProfileEntry& entry = PROFILE_TABLE[i]; + if (kernel && strcmp(kernel, entry.name) != 0) + continue; + printf("Profiling %-14s S=%d\n", entry.name, S); + entry.fn(S, h_A, h_B); + ran_any = true; + } + + if (!ran_any) { + fprintf(stderr, "Unknown kernel '%s'. Available: ", kernel); + for (int i = 0; i < PROFILE_TABLE_SIZE; ++i) + fprintf(stderr, "%s%s", PROFILE_TABLE[i].name, + i + 1 < PROFILE_TABLE_SIZE ? ", " : "\n"); + delete[] h_A; + delete[] h_B; + return EXIT_FAILURE; + } + + delete[] h_A; + delete[] h_B; + return 0; +} diff --git a/cuda-profiling/scripts/download_vendor.sh b/cuda-profiling/scripts/download_vendor.sh new file mode 100755 index 0000000..773bd02 --- /dev/null +++ b/cuda-profiling/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." diff --git a/cuda-profiling/src/kernels.cu b/cuda-profiling/src/kernels.cu new file mode 100644 index 0000000..d05e262 --- /dev/null +++ b/cuda-profiling/src/kernels.cu @@ -0,0 +1,193 @@ +#include + +#include "kernels.cuh" + +// Single Precision Matrix Multiplication Kernels SGEMM: C = alpha * A * B + beta * C + + +// Naive implementation: 1 thread per output element, no shared memory, non-coalesced accesses. +__global__ void sgemm_naive( + 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) + { + + const uint x = blockIdx.x * blockDim.x + threadIdx.x; // 0 .. M-1 (Rows of C) + const uint y = blockIdx.y * blockDim.y + threadIdx.y; // 0 .. N-1 (Columns of C) + + // Warp access e.g. warp 0: + // x is consecutive for threads 0..31, y is the same for threads 0..31 + // each warp accesses a column of C, with non-coalesced accesses to A and C, but broadcast access to B + // Ref: https://siboehm.com/articles/22/CUDA-MMM + + if (x >= M || y >= N) return; + + float acc = 0.0f; + for (size_t i = 0; i < K; ++i){ + acc += A[x * K + i] * B[i * N + y]; + } + + C[x * N + y] = alpha * acc + beta * C[x * N + y]; +} + + +// Coalesced access version: 1 thread per output element, no shared memory, but coalesced accesses to A and B. +template +__global__ void sgemm_coalesced( + size_t M, + size_t N, + size_t K, + float alpha, + const float *A, + const float *B, + float beta, + float *C) { + + const uint x = blockIdx.x * BLOCKSIZE + (threadIdx.x / BLOCKSIZE); // 0 .. M-1 (Rows of C) + const uint y = blockIdx.y * BLOCKSIZE + (threadIdx.x % BLOCKSIZE); // 0 .. N-1 (Columns of C) + + // Warp access e.g. warp 0: + // x is the same for threads 0..31, y is consecutive for threads 0..31 + // each warp accesses a BLOCKSIZE (32) tile of C, with coalesced accesses to A and C but broadcast acces to B + // Ref: https://siboehm.com/articles/22/CUDA-MMM + + if (x >= M || y >= N) return; + + float acc = 0.0f; + for (size_t i = 0; i < K; ++i){ + acc += A[x * K + i] * B[i * N + y]; + } + + C[x * N + y] = alpha * acc + beta * C[x * N + y]; +} + +// Tiled version: 1 thread per output element, shared memory for tiles of A and B. +template +__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) +{ + // Declaring the tiles + __shared__ float A_tile[TILE_SIZE * TILE_SIZE]; + __shared__ float B_tile[TILE_SIZE * TILE_SIZE]; + + // the element of C for this specific thread + size_t y = blockIdx.y * TILE_SIZE + threadIdx.y; // Rows of C + size_t x = blockIdx.x * TILE_SIZE + threadIdx.x; // Columns of C + + // we multiply along the K dimension, iterate over this + size_t n_tiles = (K + TILE_SIZE -1)/TILE_SIZE; + + float acc = .0f; + for(size_t tile=0; tile(size_t M, size_t N, size_t K, float alpha, + const float *A, const float *B, + float beta, float *C); +template __global__ void sgemm_tiled<32>(size_t M, size_t N, size_t K, float alpha, + const float *A, const float *B, + float beta, float *C); +template __global__ void sgemm_coalesced<16>(size_t M, size_t N, size_t K, float alpha, + const float *A, const float *B, + float beta, float *C); +template __global__ void sgemm_coalesced<32>(size_t M, size_t N, size_t K, float alpha, + const float *A, const float *B, + float beta, float *C); + + +// Thread-coarsened tiled SGEMM. +// Block tile BM×BN, K-step BK, each thread handles TM rows. +// Constraint: BM*BK == blockDim.x (512 with defaults BM=64,BN=64,BK=8,TM=8) +// i.e. BM*BN/TM threads, each loading one element of As and one of Bs. +template +__global__ void sgemm_coarsened( + size_t M, size_t N, size_t K, + float alpha, const float *A, const float *B, + float beta, float *C) +{ + const int blockRow = blockIdx.y; + const int blockCol = blockIdx.x; + + // Each thread's output position within the block tile + const int threadRow = threadIdx.x / BN; // row group [0, BM/TM) + const int threadCol = threadIdx.x % BN; // column [0, BN) + + // Loading indices: each thread loads one element into As and one into Bs + const int innerRowA = threadIdx.x / BK; // [0, BM) + const int innerColA = threadIdx.x % BK; // [0, BK) + const int innerRowB = threadIdx.x / BN; // [0, BK) + const int innerColB = threadIdx.x % BN; // [0, BN) + + __shared__ float As[BM * BK]; + __shared__ float Bs[BK * BN]; + + float threadResults[TM] = {}; + + for (int tileIdx = 0; tileIdx < (int)((K + BK - 1) / BK); ++tileIdx) { + // Load one element of A tile + size_t gRowA = (size_t)blockRow * BM + innerRowA; + size_t gColA = (size_t)tileIdx * BK + innerColA; + As[innerRowA * BK + innerColA] = (gRowA < M && gColA < K) ? A[gRowA * K + gColA] : 0.0f; + + // Load one element of B tile + size_t gRowB = (size_t)tileIdx * BK + innerRowB; + size_t gColB = (size_t)blockCol * BN + innerColB; + Bs[innerRowB * BN + innerColB] = (gRowB < K && gColB < N) ? B[gRowB * N + gColB] : 0.0f; + + __syncthreads(); + + // Accumulate: load Bs value once, reuse across all TM rows + #pragma unroll + for (int dotIdx = 0; dotIdx < BK; ++dotIdx) { + float bVal = Bs[dotIdx * BN + threadCol]; + #pragma unroll + for (int tm = 0; tm < TM; ++tm) + threadResults[tm] += As[(threadRow * TM + tm) * BK + dotIdx] * bVal; + } + + __syncthreads(); + } + + // Write TM results to global memory + #pragma unroll + for (int tm = 0; tm < TM; ++tm) { + size_t gRow = (size_t)blockRow * BM + threadRow * TM + tm; + size_t gCol = (size_t)blockCol * BN + threadCol; + if (gRow < M && gCol < N) + C[gRow * N + gCol] = alpha * threadResults[tm] + beta * C[gRow * N + gCol]; + } +} + +template __global__ void sgemm_coarsened<64, 64, 8, 8>(size_t M, size_t N, size_t K, + float alpha, const float *A, + const float *B, float beta, float *C); diff --git a/cuda-profiling/src/profile_cutlass.cu b/cuda-profiling/src/profile_cutlass.cu new file mode 100644 index 0000000..4b67234 --- /dev/null +++ b/cuda-profiling/src/profile_cutlass.cu @@ -0,0 +1,219 @@ +#include +#include +#include + +#include +#include + +#include "cuda_check.cuh" +#include "profiling.h" + +#define WARMUP 3 +#define REPEAT 3 + +// CUTLASS version using FP16 Tensor Cores; A/B are narrowed to half precision on host before upload. +void profile_cutlass_fp16(int S, const float* h_A, const float* h_B) +{ + using CutlassGemm = cutlass::gemm::device::Gemm< + cutlass::half_t, cutlass::layout::RowMajor, // A (fp16 half precision) and row major + cutlass::half_t, cutlass::layout::RowMajor, // B (fp16 half precision) and row major + float, cutlass::layout::RowMajor, // C / D (full float output) + float, // accumulator + cutlass::arch::OpClassTensorOp, // FP16 tensor cores + cutlass::arch::Sm80 // Ampere architecture for our A100 GPU + >; + + size_t n = (size_t)S * S; + size_t szAB = n * sizeof(cutlass::half_t); + size_t szC = n * sizeof(float); + + // Narrow the float operands to half on the host (outside the profiled region). + std::vector h_A_half(n), h_B_half(n); + for (size_t i = 0; i < n; ++i) { + h_A_half[i] = cutlass::half_t(h_A[i]); + h_B_half[i] = cutlass::half_t(h_B[i]); + } + + cutlass::half_t *d_A, *d_B; + float *d_C; + CUDA_CHECK(cudaMalloc(&d_A, szAB)); + CUDA_CHECK(cudaMalloc(&d_B, szAB)); + CUDA_CHECK(cudaMalloc(&d_C, szC)); + + CUDA_CHECK(cudaMemcpy(d_A, h_A_half.data(), szAB, cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemcpy(d_B, h_B_half.data(), szAB, cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemset(d_C, 0, szC)); + + const float alpha = 1.0f, beta = 0.0f; + + CutlassGemm gemm_op; + CutlassGemm::Arguments args( + {S, S, S}, // problem size M, N, K + {d_A, S}, // A ref: pointer + leading dim + {d_B, S}, // B ref + {d_C, S}, // C ref (source for beta*C) + {d_C, S}, // D ref (output destination) + {alpha, beta} // epilogue: alpha*A*B + beta*C + ); + + size_t workspace_bytes = gemm_op.get_workspace_size(args); + void* workspace = nullptr; + if (workspace_bytes) + CUDA_CHECK(cudaMalloc(&workspace, workspace_bytes)); + + for (int i = 0; i < WARMUP; ++i) { + cutlass::Status s = gemm_op(args, workspace); + if (s != cutlass::Status::kSuccess) { + fprintf(stderr, "CUTLASS error %d\n", (int)s); + exit(EXIT_FAILURE); + } + } + CUDA_CHECK(cudaDeviceSynchronize()); + + { + NvtxRange range("cutlass_fp16"); + for (int i = 0; i < REPEAT; ++i) + gemm_op(args, workspace); + CUDA_CHECK(cudaDeviceSynchronize()); + } + + if (workspace) CUDA_CHECK(cudaFree(workspace)); + CUDA_CHECK(cudaFree(d_A)); + CUDA_CHECK(cudaFree(d_B)); + CUDA_CHECK(cudaFree(d_C)); +} + +// CUTLASS version using TF32 Tensor Cores; A/B are bitcast (still 4 bytes), not narrowed. +void profile_cutlass_tf32(int S, const float* h_A, const float* h_B) +{ + // tfloat32_t is bit-identical width to float (4 bytes); float inputs are + // bitcast, not copied, which truncates the mantissa to 10 bits for the MMA. + // Explicit tile shapes are required: CUTLASS's generic Sm80/tf32 default + // config hits a broken epilogue dispatch (FragmentIteratorComplexTensorOp) + // in this CUTLASS version. + using CutlassGemmTF32 = cutlass::gemm::device::Gemm< + cutlass::tfloat32_t, cutlass::layout::RowMajor, // A (tf32 ~ float with 10-bit mantissa) + cutlass::tfloat32_t, cutlass::layout::RowMajor, // B + float, cutlass::layout::RowMajor, // C / D (full float output) + float, // accumulator + cutlass::arch::OpClassTensorOp, // TF32 tensor cores + cutlass::arch::Sm80, // Ampere + cutlass::gemm::GemmShape<128, 128, 32>, // threadblock tile + cutlass::gemm::GemmShape<64, 64, 32>, // warp tile + cutlass::gemm::GemmShape<16, 8, 8> // TF32 MMA instruction shape + >; + + size_t n = (size_t)S * S; + size_t sz = n * sizeof(float); + + float *d_A, *d_B, *d_C; + CUDA_CHECK(cudaMalloc(&d_A, sz)); + CUDA_CHECK(cudaMalloc(&d_B, sz)); + CUDA_CHECK(cudaMalloc(&d_C, sz)); + CUDA_CHECK(cudaMemcpy(d_A, h_A, sz, cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemcpy(d_B, h_B, sz, cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemset(d_C, 0, sz)); + + const float alpha = 1.0f, beta = 0.0f; + + // Reinterpret float* as tfloat32_t* — same bit width, bitcast is intentional. + auto* tf32_A = reinterpret_cast(d_A); + auto* tf32_B = reinterpret_cast(d_B); + + CutlassGemmTF32 gemm_op; + CutlassGemmTF32::Arguments args( + {S, S, S}, + {tf32_A, S}, + {tf32_B, S}, + {d_C, S}, + {d_C, S}, + {alpha, beta} + ); + + size_t workspace_bytes = gemm_op.get_workspace_size(args); + void* workspace = nullptr; + if (workspace_bytes) + CUDA_CHECK(cudaMalloc(&workspace, workspace_bytes)); + + for (int i = 0; i < WARMUP; ++i) { + cutlass::Status s = gemm_op(args, workspace); + if (s != cutlass::Status::kSuccess) { + fprintf(stderr, "CUTLASS error %d\n", (int)s); + exit(EXIT_FAILURE); + } + } + CUDA_CHECK(cudaDeviceSynchronize()); + + { + NvtxRange range("cutlass_tf32"); + for (int i = 0; i < REPEAT; ++i) + gemm_op(args, workspace); + CUDA_CHECK(cudaDeviceSynchronize()); + } + + if (workspace) CUDA_CHECK(cudaFree(workspace)); + CUDA_CHECK(cudaFree(d_A)); + CUDA_CHECK(cudaFree(d_B)); + CUDA_CHECK(cudaFree(d_C)); +} + +// CUTLASS version using regular fp32 CUDA cores (no precision narrowing, no tensor cores). +void profile_cutlass_fp32(int S, const float* h_A, const float* h_B) +{ + using CutlassGemmFP32 = cutlass::gemm::device::Gemm< + float, cutlass::layout::RowMajor, // A and row major + float, cutlass::layout::RowMajor, // B and row major + float, cutlass::layout::RowMajor, // C / D (full float output) + float, // accumulator + cutlass::arch::OpClassSimt, // CUDA cores, not tensor cores + cutlass::arch::Sm80 // Ampere architecture for our A100 GPU + >; + size_t n = (size_t)S * S; + size_t sz = n * sizeof(float); + + float *d_A, *d_B, *d_C; + CUDA_CHECK(cudaMalloc(&d_A, sz)); + CUDA_CHECK(cudaMalloc(&d_B, sz)); + CUDA_CHECK(cudaMalloc(&d_C, sz)); + CUDA_CHECK(cudaMemcpy(d_A, h_A, sz, cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemcpy(d_B, h_B, sz, cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemset(d_C, 0, sz)); + + const float alpha = 1.0f, beta = 0.0f; + + CutlassGemmFP32 gemm_op; + CutlassGemmFP32::Arguments args( + {S, S, S}, + {d_A, S}, + {d_B, S}, + {d_C, S}, + {d_C, S}, + {alpha, beta} + ); + + size_t workspace_bytes = gemm_op.get_workspace_size(args); + void* workspace = nullptr; + if (workspace_bytes) + CUDA_CHECK(cudaMalloc(&workspace, workspace_bytes)); + + for (int i = 0; i < WARMUP; ++i) { + cutlass::Status s = gemm_op(args, workspace); + if (s != cutlass::Status::kSuccess) { + fprintf(stderr, "CUTLASS error %d\n", (int)s); + exit(EXIT_FAILURE); + } + } + CUDA_CHECK(cudaDeviceSynchronize()); + + { + NvtxRange range("cutlass_fp32"); + for (int i = 0; i < REPEAT; ++i) + gemm_op(args, workspace); + CUDA_CHECK(cudaDeviceSynchronize()); + } + + if (workspace) CUDA_CHECK(cudaFree(workspace)); + CUDA_CHECK(cudaFree(d_A)); + CUDA_CHECK(cudaFree(d_B)); + CUDA_CHECK(cudaFree(d_C)); +} diff --git a/cuda-profiling/src/profile_kernels.cu b/cuda-profiling/src/profile_kernels.cu new file mode 100644 index 0000000..d53a105 --- /dev/null +++ b/cuda-profiling/src/profile_kernels.cu @@ -0,0 +1,224 @@ +#include +#include + +#include + +#include "cuda_check.cuh" +#include "kernels.cuh" +#include "profiling.h" +#include "utils.h" + +#define CUBLAS_CHECK(call) \ + do { \ + cublasStatus_t _s = (call); \ + if (_s != CUBLAS_STATUS_SUCCESS) { \ + fprintf(stderr, "cuBLAS error %d at %s:%d\n", _s, \ + __FILE__, __LINE__); \ + exit(EXIT_FAILURE); \ + } \ + } while (0) + +// Untimed launches to warm up clocks/caches before the profiled region. +#define WARMUP 3 +// Launches inside the NVTX range, so the timeline/report shows more than +// one occurrence. `ncu` defaults to profiling all of them; pass +// `--launch-count 1` (done by the Makefile) to only profile the first. +#define REPEAT 3 + +void profile_naive(int S, const float* h_A, const float* h_B) +{ + size_t sz = (size_t)S * S * sizeof(float); + + float *d_A, *d_B, *d_C; + CUDA_CHECK(cudaMalloc(&d_A, sz)); + CUDA_CHECK(cudaMalloc(&d_B, sz)); + CUDA_CHECK(cudaMalloc(&d_C, sz)); + CUDA_CHECK(cudaMemcpy(d_A, h_A, sz, cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemcpy(d_B, h_B, sz, cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemset(d_C, 0, sz)); + + dim3 gridDim(CEIL_DIV(S, 32), CEIL_DIV(S, 32)); + dim3 blockDim(32, 32, 1); + + for (int i = 0; i < WARMUP; ++i) + sgemm_naive<<>>(S, S, S, 1.0f, d_A, d_B, 0.0f, d_C); + CUDA_CHECK(cudaDeviceSynchronize()); + + { + NvtxRange range("naive"); + for (int i = 0; i < REPEAT; ++i) + sgemm_naive<<>>(S, S, S, 1.0f, d_A, d_B, 0.0f, d_C); + CUDA_CHECK(cudaDeviceSynchronize()); + } + + CUDA_CHECK(cudaFree(d_A)); + CUDA_CHECK(cudaFree(d_B)); + CUDA_CHECK(cudaFree(d_C)); +} + +void profile_coalesced(int S, const float* h_A, const float* h_B) +{ + size_t sz = (size_t)S * S * sizeof(float); + + float *d_A, *d_B, *d_C; + CUDA_CHECK(cudaMalloc(&d_A, sz)); + CUDA_CHECK(cudaMalloc(&d_B, sz)); + CUDA_CHECK(cudaMalloc(&d_C, sz)); + CUDA_CHECK(cudaMemcpy(d_A, h_A, sz, cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemcpy(d_B, h_B, sz, cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemset(d_C, 0, sz)); + + dim3 gridDim(CEIL_DIV(S, 32), CEIL_DIV(S, 32)); + dim3 blockDim(32 * 32); + + for (int i = 0; i < WARMUP; ++i) + sgemm_coalesced<32><<>>(S, S, S, 1.0f, d_A, d_B, 0.0f, d_C); + CUDA_CHECK(cudaDeviceSynchronize()); + + { + NvtxRange range("coalesced"); + for (int i = 0; i < REPEAT; ++i) + sgemm_coalesced<32><<>>(S, S, S, 1.0f, d_A, d_B, 0.0f, d_C); + CUDA_CHECK(cudaDeviceSynchronize()); + } + + CUDA_CHECK(cudaFree(d_A)); + CUDA_CHECK(cudaFree(d_B)); + CUDA_CHECK(cudaFree(d_C)); +} + +void profile_tiled(int S, const float* h_A, const float* h_B) +{ + size_t sz = (size_t)S * S * sizeof(float); + + float *d_A, *d_B, *d_C; + CUDA_CHECK(cudaMalloc(&d_A, sz)); + CUDA_CHECK(cudaMalloc(&d_B, sz)); + CUDA_CHECK(cudaMalloc(&d_C, sz)); + CUDA_CHECK(cudaMemcpy(d_A, h_A, sz, cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemcpy(d_B, h_B, sz, cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemset(d_C, 0, sz)); + + dim3 blockDim(32, 32, 1); + dim3 gridDim(CEIL_DIV(S, blockDim.x), CEIL_DIV(S, blockDim.y)); + + for (int i = 0; i < WARMUP; ++i) + sgemm_tiled<32><<>>(S, S, S, 1.0f, d_A, d_B, 0.0f, d_C); + CUDA_CHECK(cudaDeviceSynchronize()); + + { + NvtxRange range("tiled"); + for (int i = 0; i < REPEAT; ++i) + sgemm_tiled<32><<>>(S, S, S, 1.0f, d_A, d_B, 0.0f, d_C); + CUDA_CHECK(cudaDeviceSynchronize()); + } + + CUDA_CHECK(cudaFree(d_A)); + CUDA_CHECK(cudaFree(d_B)); + CUDA_CHECK(cudaFree(d_C)); +} + +void profile_coarsened(int S, const float* h_A, const float* h_B) +{ + constexpr int BM = 64, BN = 64, BK = 8, TM = 8; + + size_t sz = (size_t)S * S * sizeof(float); + + float *d_A, *d_B, *d_C; + CUDA_CHECK(cudaMalloc(&d_A, sz)); + CUDA_CHECK(cudaMalloc(&d_B, sz)); + CUDA_CHECK(cudaMalloc(&d_C, sz)); + CUDA_CHECK(cudaMemcpy(d_A, h_A, sz, cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemcpy(d_B, h_B, sz, cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemset(d_C, 0, sz)); + + dim3 gridDim(CEIL_DIV(S, BN), CEIL_DIV(S, BM)); + dim3 blockDim(BM * BN / TM); // 512 threads + + for (int i = 0; i < WARMUP; ++i) + sgemm_coarsened<<>>(S, S, S, 1.0f, d_A, d_B, 0.0f, d_C); + CUDA_CHECK(cudaDeviceSynchronize()); + + { + NvtxRange range("coarsened"); + for (int i = 0; i < REPEAT; ++i) + sgemm_coarsened<<>>(S, S, S, 1.0f, d_A, d_B, 0.0f, d_C); + CUDA_CHECK(cudaDeviceSynchronize()); + } + + CUDA_CHECK(cudaFree(d_A)); + CUDA_CHECK(cudaFree(d_B)); + CUDA_CHECK(cudaFree(d_C)); +} + +// Persistent handle for cuBLAS +static cublasHandle_t cublas_handle() +{ + static cublasHandle_t h = []() { + cublasHandle_t handle; + CUBLAS_CHECK(cublasCreate(&handle)); + CUBLAS_CHECK(cublasSetMathMode(handle, CUBLAS_TF32_TENSOR_OP_MATH)); + return handle; + }(); + return h; +} + +void profile_cublas(int S, const float* h_A, const float* h_B) +{ + size_t sz = (size_t)S * S * sizeof(float); + + float *d_A, *d_B, *d_C; + CUDA_CHECK(cudaMalloc(&d_A, sz)); + CUDA_CHECK(cudaMalloc(&d_B, sz)); + CUDA_CHECK(cudaMalloc(&d_C, sz)); + CUDA_CHECK(cudaMemcpy(d_A, h_A, sz, cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemcpy(d_B, h_B, sz, cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemset(d_C, 0, sz)); + + cublasHandle_t handle = cublas_handle(); + + // cuBLAS is column-major. For row-major C=A*B, use the identity + // C^T = B^T * A^T, so pass B first with leading dimension S. + const float alpha = 1.0f, beta = 0.0f; + + for (int i = 0; i < WARMUP; ++i) + CUBLAS_CHECK(cublasSgemm(handle, + CUBLAS_OP_N, CUBLAS_OP_N, + S, S, S, + &alpha, + d_B, S, + d_A, S, + &beta, + d_C, S)); + CUDA_CHECK(cudaDeviceSynchronize()); + + { + NvtxRange range("cublas"); + for (int i = 0; i < REPEAT; ++i) + CUBLAS_CHECK(cublasSgemm(handle, + CUBLAS_OP_N, CUBLAS_OP_N, + S, S, S, + &alpha, + d_B, S, + d_A, S, + &beta, + d_C, S)); + CUDA_CHECK(cudaDeviceSynchronize()); + } + + CUDA_CHECK(cudaFree(d_A)); + CUDA_CHECK(cudaFree(d_B)); + CUDA_CHECK(cudaFree(d_C)); +} + +const ProfileEntry PROFILE_TABLE[] = { + {"naive", profile_naive}, + {"coalesced", profile_coalesced}, + {"tiled", profile_tiled}, + {"coarsened", profile_coarsened}, + {"cublas", profile_cublas}, + {"cutlass_fp32", profile_cutlass_fp32}, + {"cutlass_tf32", profile_cutlass_tf32}, + {"cutlass_fp16", profile_cutlass_fp16}, +}; +const int PROFILE_TABLE_SIZE = sizeof(PROFILE_TABLE) / sizeof(PROFILE_TABLE[0]); diff --git a/cuda-profiling/src/utils.cu b/cuda-profiling/src/utils.cu new file mode 100644 index 0000000..e132663 --- /dev/null +++ b/cuda-profiling/src/utils.cu @@ -0,0 +1,6 @@ +#include "utils.h" + +void fill_random(float* p, int n) { + for (int i = 0; i < n; ++i) + p[i] = (float)rand() / RAND_MAX * 2.0f - 1.0f; +} From 2489ad9108475eefe7116a82277fc6fe9eb86814 Mon Sep 17 00:00:00 2001 From: Sebas <12949040+SebastiaAgramunt@users.noreply.github.com> Date: Thu, 9 Jul 2026 14:46:48 +0000 Subject: [PATCH 19/21] Untrack cuda-performance/ Stop version-controlling this directory; files remain on disk. Co-Authored-By: Claude Sonnet 5 --- .gitignore | 3 + cuda-performance/README.md | 31 --- cuda-performance/include/utils.h | 34 --- cuda-performance/include/vector_add.h | 17 -- cuda-performance/scripts/analyze.py | 174 -------------- cuda-performance/scripts/compile.sh | 64 ----- cuda-performance/scripts/execute.sh | 26 --- .../scripts/install_python_env.sh | 10 - cuda-performance/src/main.cu | 220 ------------------ cuda-performance/src/vector_add.cu | 63 ----- 10 files changed, 3 insertions(+), 639 deletions(-) delete mode 100644 cuda-performance/README.md delete mode 100644 cuda-performance/include/utils.h delete mode 100644 cuda-performance/include/vector_add.h delete mode 100644 cuda-performance/scripts/analyze.py delete mode 100755 cuda-performance/scripts/compile.sh delete mode 100755 cuda-performance/scripts/execute.sh delete mode 100755 cuda-performance/scripts/install_python_env.sh delete mode 100644 cuda-performance/src/main.cu delete mode 100644 cuda-performance/src/vector_add.cu diff --git a/.gitignore b/.gitignore index 0d66ca8..28b1a29 100644 --- a/.gitignore +++ b/.gitignore @@ -188,3 +188,6 @@ external # results and analysis *.deb + +# untracked local project +cuda-performance/ diff --git a/cuda-performance/README.md b/cuda-performance/README.md deleted file mode 100644 index b6622cf..0000000 --- a/cuda-performance/README.md +++ /dev/null @@ -1,31 +0,0 @@ -# CUDA performance on vector addition and matrix multiplication - -Supporting material for blogpost [agramunt.me/posts/cuda-performance/](https://agramunt.me/posts/cuda-performance/). - -## Compile & Run - -Compile the code (all intermediate objects will be placed in `build` directory) - -```bash -./scripts/compile.sh -``` - -Execute the benchmark - -```bash -./scripts/execute.sh -``` - -## Analyze data and generate plots - -First create a new python environment and install depencencies, you can do that with - -```bash -./scripts/install_python_env.sh -``` - -That will create an environment in the root directory `.venv`. Then just execute the script with - -```bash -.venv/bin/python scripts/analyze.py -``` diff --git a/cuda-performance/include/utils.h b/cuda-performance/include/utils.h deleted file mode 100644 index e2596c5..0000000 --- a/cuda-performance/include/utils.h +++ /dev/null @@ -1,34 +0,0 @@ -#ifndef UTILS_H -#define UTILS_H - -#include -#include - -inline void check(cudaError_t err, const char* const func, const char* const file, - const int line) -{ - if (err != cudaSuccess) - { - std::cerr << "CUDA Runtime Error at: " << file << ":" << line - << std::endl; - std::cerr << cudaGetErrorString(err) << " " << func << std::endl; - std::exit(EXIT_FAILURE); - } -} - -inline void checkLast(const char* const file, const int line) -{ - cudaError_t const err{cudaGetLastError()}; - if (err != cudaSuccess) - { - std::cerr << "CUDA Runtime Error at: " << file << ":" << line - << std::endl; - std::cerr << cudaGetErrorString(err) << std::endl; - std::exit(EXIT_FAILURE); - } -} - -#define CHECK_CUDA_ERROR(expr) (check((expr), #expr, __FILE__, __LINE__)) -#define CHECK_LAST_CUDA_ERROR() (checkLast(__FILE__, __LINE__)) - -#endif \ No newline at end of file diff --git a/cuda-performance/include/vector_add.h b/cuda-performance/include/vector_add.h deleted file mode 100644 index 09337e3..0000000 --- a/cuda-performance/include/vector_add.h +++ /dev/null @@ -1,17 +0,0 @@ -#ifndef VECTOR_ADD_H -#define VECTOR_ADD_H - -#include "utils.h" - -template -void vectorAdd_wrapper(const T* __restrict__ d_a, - const T* __restrict__ d_b, - T* __restrict__ d_c, - int N, - int m = 1, - int ThreadsPerBlock = 256); - -template -void dummyCalculation(); - -#endif \ No newline at end of file diff --git a/cuda-performance/scripts/analyze.py b/cuda-performance/scripts/analyze.py deleted file mode 100644 index af4f1fa..0000000 --- a/cuda-performance/scripts/analyze.py +++ /dev/null @@ -1,174 +0,0 @@ -#!/usr/bin/env python - -import pandas as pd -import numpy as np -import matplotlib.pyplot as plt -import pathlib - -THIS_DIR = pathlib.Path(__file__).parent.resolve() -ROOT_DIR = THIS_DIR.parent.resolve() - - -def plot_cpu_vs_gpu(df: pd.DataFrame): - fig, ax = plt.subplots(1, 1, dpi=100, layout='constrained', figsize=(15, 8)) - - additions = [1, 256, 4096, 16384] - tags = ["1", "$2^8$", "$2^{11}$", "$2^{14}$"] - colors = plt.cm.viridis(np.linspace(0, 1, len(additions))) - for i, key in enumerate(additions): - df_tmp = df.loc[df["number_of_additions"] == key, :] - plt.plot(np.log10(df_tmp["N"]), df_tmp["total_GPU_time"], marker='o', label=f"GPU aditions={tags[i]}", color=colors[i], lw=3) - - number_of_additions = 1 - df_tmp = df.loc[df["number_of_additions"]==number_of_additions, :] - plt.plot(np.log10(df_tmp["N"]), df_tmp["total_CPU_time"], marker='o', label=f"CPU aditions={number_of_additions}", color="black", lw=5) - - number_of_additions = 256 - df_tmp = df.loc[df["number_of_additions"]==number_of_additions, :] - plt.plot(np.log10(df_tmp["N"]), df_tmp["total_CPU_time"], marker='o', label=f"CPU aditions=$2^8$", color="gray", lw=5) - - plt.xlabel("$\log_{10}(N)$", fontsize=16) - plt.ylabel("Time ($ms$)", fontsize=16) - plt.title("Computation time as function of number of elements in the vectors", fontsize=18) - plt.grid(True) - - ax.set_xlim(4, 9) - ax.set_ylim(-110, 1500) - ax.tick_params(axis="both", which="major", labelsize=14) - ax.legend(fontsize=20) - ax.grid(True, which="both", ls="--", lw=2.0, alpha=0.5) - fig.savefig(ROOT_DIR / "results" / "gpu_cpu_performance.png") - plt.close(fig) - - -def plot_cpu_vs_gpu_1(df: pd.DataFrame): - fig, ax = plt.subplots(1, 1, dpi=100, layout='constrained', figsize=(15, 8)) - - number_of_additions = 1 - df_tmp = df.loc[df["number_of_additions"]==number_of_additions, :] - - plt.plot(np.log10(df_tmp["N"]), df_tmp["total_GPU_time"], marker='o', label=f"GPU", color="green", lw=3) - plt.plot(np.log10(df_tmp["N"]), df_tmp["total_CPU_time"], marker='o', label=f"CPU", color="black", lw=3) - plt.vlines(x=6.5, ymin=-1.0, ymax=25, color="gray", lw=4) - ax.text(x=6.15, y=6, s="vector length\n~3.1M floats", fontsize=14, alpha=1) - ax.text(x=6.15, y=22.5, s="CPU is faster\nthan GPU", fontsize=14, alpha=1) - ax.text(x=6.65, y=22.5, s="GPU is faster\nthan CPU", fontsize=14, alpha=1) - - plt.xlabel("$\log_{10}(N)$", fontsize=16) - plt.ylabel("Time ($ms$)", fontsize=16) - plt.title("Total time for one vector adition", fontsize=18) - plt.grid(True) - # plt.tight_layout() - - ax.set_xlim(5, 8) - ax.set_ylim(-1, 25) - ax.tick_params(axis="both", which="major", labelsize=14) - ax.legend(fontsize=20) - ax.grid(True, which="both", ls="--", lw=2.0, alpha=0.5) - fig.savefig(ROOT_DIR / "results" / "gpu_cpu_performance_1.png") - plt.close(fig) - -def plot_gpu_times(df, mult_times=2048): - - # df_tmp = df.loc[df["number_of_additions"] == 1, :] - df_tmp = df.loc[df["number_of_additions"] == mult_times, :] - # df_tmp = df.loc[df["number_of_additions"] == 2048, :] - x = np.log10(df_tmp["N"].to_numpy()) - alloc = df_tmp["allocate_time"].to_numpy() - copy_h2d = df_tmp["copy_H2D_time"].to_numpy() - compute = df_tmp["compute_time"].to_numpy() - copy_d2h = df_tmp["copy_D2H_time"].to_numpy() - free_ = df_tmp["free_time"].to_numpy() - - series = [alloc, copy_h2d, compute, copy_d2h, free_] - labels = ["Allocation", "Copy H→D", "Compute", "Copy D→H", "Free"] - - fig, ax = plt.subplots(1, 1, dpi=100, layout='constrained', figsize=(15, 8)) - colors = plt.cm.viridis(np.linspace(0, 1, len(series))) - - ax.stackplot( - x, - series, - labels=labels, - colors=colors, - linewidth=0.0, - alpha=0.9, - ) - - # # overlay total time as a line (use your column if you have it) - total = df_tmp["total_time_gpu"].to_numpy() if "total_time_gpu" in df_tmp.columns else np.sum(series, axis=0) - ax.plot(x, total, label="Total GPU Time", color="red", lw=4) - - # # plot cpu time but filter by small times - # cpu_time = df_tmp.loc[df["N"] < 1024, "compute_time"].to_numpy() - # cpu_time_x = df_tmp.loc[df["N"] < 1024, "N"].to_numpy() - # ax.plot(cpu_time_x, cpu_time, label="CPU Time", lw=4, color="black") - - # ax.set_xlim(xrange[0], xrange[1]) - # ax.set_ylim(yrange[0], yrange[1]) - # ax.set_xlim(0, 1.0e9) - ax.set_ylim(0, 1000) - ax.set_xlim(6, 9) - ax.set_title(r"vector addition " + f"{mult_times} times", fontsize=18) - ax.tick_params(axis="both", which="major", labelsize=14) - ax.set_xlabel("$\log_{10} N$", fontsize=16) - ax.set_ylabel("time (ms)", fontsize=16) - ax.legend(fontsize=16, loc="best", ncol=2, frameon=True) - ax.grid(True, which="both", ls="--", lw=2.0, alpha=0.5) - - filename=f"gpu_times_stacked_{mult_times}.png" - fig.savefig(ROOT_DIR / "results" / filename) - plt.close(fig) - -def plot_percentages_bar(df, number_of_additions=1, N=512): - metrics = ["allocate_time", "copy_H2D_time", "compute_time", "copy_D2H_time", "free_time"] - times = df.loc[(df["N"] == N) & (df["number_of_additions"]==number_of_additions), metrics].values.flatten() - - total = sum(times) - - fig, ax = plt.subplots(figsize=(5,8)) - - total = sum(times) - bottom = 0 - colors = plt.cm.viridis(np.linspace(0, 1, len(metrics))) - - for metric, time, color in zip(metrics, times, colors): - ax.bar(0, time, bottom=bottom, color=color, width=0.6, - label=f"{metric}: {(time/total)*100:.1f}%") - bottom += time - - leg = ax.legend( - title="Times", - loc="upper left", - bbox_to_anchor=(1.1, 1.0), # just outside the axes - fontsize=9, - title_fontsize=10, - frameon=False, - ) - - # Make room on the right for the legend - fig.subplots_adjust(right=0.5) # tweak 0.72–0.85 as needed - - # Optional clean-up - ax.set_xticks([0]) - ax.set_xticklabels(["Total Time"]) - ax.set_ylabel("Time (ms)") - ax.set_title(r"GPU Time Percentages $log_{10}$" +f"\nN={int(np.floor(np.log10(N)))}, additions {number_of_additions}") - - plt.tight_layout() - filename=f"percentages_performance_N_{N}additions_{number_of_additions}.png" - fig.savefig(ROOT_DIR / "results" / filename) - plt.close(fig) - -def main(): - df = pd.read_csv(ROOT_DIR / "results" / "results.csv") - plot_cpu_vs_gpu_1(df) - plot_cpu_vs_gpu(df) - plot_gpu_times(df, mult_times=1) - plot_gpu_times(df, mult_times=1024) - plot_gpu_times(df, mult_times=2048) - plot_percentages_bar(df, number_of_additions=2048, N=1073741824) - plot_percentages_bar(df, number_of_additions=2048, N=512) - -if __name__ == "__main__": - main() \ No newline at end of file diff --git a/cuda-performance/scripts/compile.sh b/cuda-performance/scripts/compile.sh deleted file mode 100755 index 3ea3cf4..0000000 --- a/cuda-performance/scripts/compile.sh +++ /dev/null @@ -1,64 +0,0 @@ -#!/bin/bash - -THIS_DIR=$(dirname "$(realpath "$0")") -ROOT_DIR=$(dirname ${THIS_DIR}) - -recreate_dirs(){ - # removing build directory - echo "Removing ${ROOT_DIR}/build and recreating..." - rm -rf ${ROOT_DIR}/build - mkdir ${ROOT_DIR}/build - - # creating directories for the build - mkdir ${ROOT_DIR}/build/obj - mkdir ${ROOT_DIR}/build/bin - mkdir ${ROOT_DIR}/build/lib -} - -compile_exec(){ - recreate_dirs - echo "Compiling..." - - # project includes - INCLUDES="-I${ROOT_DIR}/include" - - # cuda specific includes (modify to your cuda installation) - CUDA_INCLUDES="-I/usr/include" - CUDA_LIB_DIRS="-L/usr/lib/x86_64-linux-gnu/" - - # cuda libraries to link - CUDA_LIB="-lcudart" - - # flags. Using sm_80 for A100 GPU (compute cabability 8.0) - FLAGS="-O3 -arch=sm_80" - - # compile to objects - nvcc ${FLAGS} \ - -c ${ROOT_DIR}/src/vector_add.cu \ - ${INCLUDES} \ - ${CUDA_INCLUDES} \ - -o ${ROOT_DIR}/build/obj/vector_add.o - - nvcc ${FLAGS} \ - -c ${ROOT_DIR}/src/vector_add.cu \ - ${INCLUDES} \ - ${CUDA_INCLUDES} \ - -o ${ROOT_DIR}/build/obj/vector_add.o - - nvcc ${FLAGS} \ - -c ${ROOT_DIR}/src/main.cu \ - ${INCLUDES} \ - ${CUDA_INCLUDES} \ - -o ${ROOT_DIR}/build/obj/main.o - - - # link all the objects - nvcc ${FLAGS} \ - ${ROOT_DIR}/build/obj/vector_add.o \ - ${ROOT_DIR}/build/obj/main.o \ - ${CUDA_LIB_DIRS} \ - ${CUDA_LIB} \ - -o ${ROOT_DIR}/build/bin/main -} - -compile_exec \ No newline at end of file diff --git a/cuda-performance/scripts/execute.sh b/cuda-performance/scripts/execute.sh deleted file mode 100755 index 979798c..0000000 --- a/cuda-performance/scripts/execute.sh +++ /dev/null @@ -1,26 +0,0 @@ -#!/bin/bash - -THIS_DIR=$(dirname "$(realpath "$0")") -ROOT_DIR=$(dirname ${THIS_DIR}) - -DIRECTORY=${ROOT_DIR}/results -if [ ! -d ${DIRECTORY} ]; then - mkdir ${DIRECTORY} -fi - -FILENAME=${DIRECTORY}/results.csv -rm ${FILENAME} - -"${ROOT_DIR}/build/bin/main" --threads_per_block 512 --number_of_additions 1 --use_cpu 1 --output_file "${FILENAME}" -"${ROOT_DIR}/build/bin/main" --threads_per_block 512 --number_of_additions 256 --use_cpu 1 --output_file "${FILENAME}" - -# calculate n_additions without CPU compute (takes too long) -n_additions=(512 1024 2048 4096 8192 16384) - -for n in "${n_additions[@]}"; do - "${ROOT_DIR}/build/bin/main" \ - --threads_per_block 512 \ - --number_of_additions "$n" \ - --use_cpu 0 \ - --output_file "${FILENAME}" -done diff --git a/cuda-performance/scripts/install_python_env.sh b/cuda-performance/scripts/install_python_env.sh deleted file mode 100755 index 804b4e8..0000000 --- a/cuda-performance/scripts/install_python_env.sh +++ /dev/null @@ -1,10 +0,0 @@ -#!/bin/bash - -# not overcomplicating python environment. Use default python in the system -THIS_DIR=$(dirname "$(realpath "$0")") -ROOT_DIR=$(dirname ${THIS_DIR}) - -rm -rf ${ROOT_DIR}/.venv -python -m venv ${ROOT_DIR}/.venv -${ROOT_DIR}/.venv/bin/pip install --upgrade pip -${ROOT_DIR}/.venv/bin/pip install numpy matplotlib pandas seaborn nvitop diff --git a/cuda-performance/src/main.cu b/cuda-performance/src/main.cu deleted file mode 100644 index 0ac9510..0000000 --- a/cuda-performance/src/main.cu +++ /dev/null @@ -1,220 +0,0 @@ -#include "utils.h" -#include "vector_add.h" - -#include -#include -#include -#include - -void open_new_csv_file(const std::string& filename) { - bool new_file = !std::filesystem::exists(filename); - std::ofstream file(filename, std::ios::app); - if (!file.is_open()) { - std::cerr << "Error: could not open " << filename << " for writing.\n"; - return; - } - - if (new_file) { - file << "threads_per_block,number_of_additions,N,sizeMB,allocate_time,copy_H2D_time,compute_time,copy_D2H_time,free_time,total_GPU_time,total_CPU_time\n"; - } - - file.close(); -} - -// Append a row with timing results -void append_to_csv(const std::string& filename, - int threads_per_block, - int number_of_additions, - size_t N, - double sizeMB, - float allocateTime, - float loadH2D, - float calcTime, - float loadD2H, - float freeTime, - float totalGPUTime, - float totalCPUTime) -{ - std::ofstream file(filename, std::ios::app); - if (!file.is_open()) { - std::cerr << "Error: could not open " << filename << " for appending.\n"; - return; - } - - file << threads_per_block << "," - << number_of_additions << "," - << N << "," - << sizeMB << "," - << allocateTime << "," - << loadH2D << "," - << calcTime << "," - << loadD2H << "," - << freeTime << "," - << totalGPUTime << "," - << totalCPUTime - << "\n"; - - file.close(); -} - -int main(int argc, char **argv) { - - int threads_per_block; // threads per block - int m; // number of additions to perform - std::string filename; // file to save results - bool use_cpu = true; // use cpu by default - - using T = float; - - // arg parsing - for (int i = 1; i < argc; ++i) { - std::string arg = argv[i]; - if (arg == "--threads_per_block" && i + 1 < argc) { - threads_per_block = std::stoi(argv[++i]); - } - if (arg == "--number_of_additions" && i + 1 < argc) { - m = std::stoi(argv[++i]); - } - if (arg == "--output_file" && i + 1 < argc) { - filename = argv[++i]; - } - if (arg == "--use_cpu" && i + 1 < argc){ - std::string val = argv[++i]; - use_cpu = (val == "1" || val == "true" || val == "True"); - } - if (arg == "--help") { - std::cerr << "Usage: " << argv[0] << " --threads_per_block --vector_size --number_of_additions --use_cpu --output_file \n"; - return 1; - } - } - - std::cout << "threads_per_block: " << threads_per_block << ", number_of_additions: " << m << ", use_cpu:: " << use_cpu << std::endl; - - // vector sizes to consider - long int sizes[] = { - 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096, 8192, 16384, 32768, 65536, 131072, 262144, 524288, - 1048576, // 4 MB - 2097152, // 8 MB - 4194304, // 16 MB - 8388608, // 32 MB - 16777216, // 64 MB - 33554432, // 128 MB - 67108864, // 256 MB - 134217728, // 512 MB - 268435456, // 1 GB - 536870912, // 2 GB - 1073741824 // 4 GB - }; - - // create file if it doesn't exist - open_new_csv_file(filename); - - // Create events to log CUDA times - cudaEvent_t startAllocate, stopAllocate, startH2D, stopH2D, startCalc, stopCalc, startD2H, stopD2H, startFree, stopFree; - cudaEventCreate(&startAllocate); - cudaEventCreate(&stopAllocate); - cudaEventCreate(&startH2D); - cudaEventCreate(&stopH2D); - cudaEventCreate(&startCalc); - cudaEventCreate(&stopCalc); - cudaEventCreate(&startD2H); - cudaEventCreate(&stopD2H); - cudaEventCreate(&startFree); - cudaEventCreate(&stopFree); - - float allocateTime = 0; - float loadH2D = 0; - float calcTime = 0; - float loadD2H = 0; - float freeTime = 0; - - // dummy calculation to warm up the GPU - dummyCalculation(); - CHECK_LAST_CUDA_ERROR(); - CHECK_CUDA_ERROR(cudaDeviceSynchronize()); - - dummyCalculation(); - CHECK_LAST_CUDA_ERROR(); - CHECK_CUDA_ERROR(cudaDeviceSynchronize()); - - dummyCalculation(); - CHECK_LAST_CUDA_ERROR(); - CHECK_CUDA_ERROR(cudaDeviceSynchronize()); - - // for every size calculate the addition and record the results - for(int s=0; s < sizeof(sizes)/sizeof(sizes[0]); s++){ - long int N = sizes[s]; - - // ceate host vectors and randomize - std::vector h_a(N), h_b(N), h_c(N), h_c_cpu(N); - for (int i = 0; i < N; ++i) { - h_a[i] = static_cast(std::rand()) / RAND_MAX; - h_b[i] = static_cast(std::rand()) / RAND_MAX; - } - - // CPU time - auto cpu_start = std::chrono::high_resolution_clock::now(); - if(use_cpu){ - for (int i = 0; i < N; ++i) { - double acc = 0; - double s = h_a[i] + h_b[i]; - for (int j = 0; j < m-1; ++j) { - acc = acc + s; - } - h_c_cpu[i] = acc; - } - } - auto cpu_end = std::chrono::high_resolution_clock::now(); - double cpuTime = std::chrono::duration(cpu_end - cpu_start).count(); - - // CUDA times - // Allocate device memory - T *d_a, *d_b, *d_c; - CHECK_CUDA_ERROR(cudaEventRecord(startAllocate)); - CHECK_CUDA_ERROR(cudaMalloc(&d_a, N * sizeof(T))); - CHECK_CUDA_ERROR(cudaMalloc(&d_b, N * sizeof(T))); - CHECK_CUDA_ERROR(cudaMalloc(&d_c, N * sizeof(T))); - CHECK_CUDA_ERROR(cudaEventRecord(stopAllocate)); - CHECK_CUDA_ERROR(cudaEventSynchronize(stopAllocate)); - CHECK_CUDA_ERROR(cudaEventElapsedTime(&allocateTime, startAllocate, stopAllocate)); - - // loading data from host to device - CHECK_CUDA_ERROR(cudaEventRecord(startH2D)); - CHECK_CUDA_ERROR(cudaMemcpy(d_a, h_a.data(), N * sizeof(T), cudaMemcpyHostToDevice)); - CHECK_CUDA_ERROR(cudaMemcpy(d_b, h_b.data(), N * sizeof(T), cudaMemcpyHostToDevice)); - CHECK_CUDA_ERROR(cudaEventRecord(stopH2D)); - CHECK_CUDA_ERROR(cudaEventSynchronize(stopH2D)); - CHECK_CUDA_ERROR(cudaEventElapsedTime(&loadH2D, startH2D, stopH2D)); - - // perform calculation on GPU - CHECK_CUDA_ERROR(cudaEventRecord(startCalc)); - vectorAdd_wrapper(d_a, d_b, d_c, N, m, threads_per_block); - CHECK_LAST_CUDA_ERROR(); - CHECK_CUDA_ERROR(cudaDeviceSynchronize()); - CHECK_CUDA_ERROR(cudaEventRecord(stopCalc)); - CHECK_CUDA_ERROR(cudaEventSynchronize(stopCalc)); - CHECK_CUDA_ERROR(cudaEventElapsedTime(&calcTime, startCalc, stopCalc)); - - // load back data from device to host - CHECK_CUDA_ERROR(cudaEventRecord(startD2H)); - CHECK_CUDA_ERROR(cudaMemcpy(h_c.data(), d_c, N * sizeof(float), cudaMemcpyDeviceToHost)); - CHECK_CUDA_ERROR(cudaEventRecord(stopD2H)); - CHECK_CUDA_ERROR(cudaEventSynchronize(stopD2H)); - CHECK_CUDA_ERROR(cudaEventElapsedTime(&loadD2H, startD2H, stopD2H)); - - // free cuda memory - CHECK_CUDA_ERROR(cudaEventRecord(startFree)); - cudaFree(d_a); - cudaFree(d_b); - cudaFree(d_c); - CHECK_CUDA_ERROR(cudaEventRecord(stopFree)); - CHECK_CUDA_ERROR(cudaEventSynchronize(stopFree)); - CHECK_CUDA_ERROR(cudaEventElapsedTime(&freeTime, startFree, stopFree)); - - float totalTime = allocateTime + loadH2D + calcTime + loadD2H + freeTime; - float sizeMB = N * sizeof(T) / (1024.0 * 1024.0); - - append_to_csv(filename, threads_per_block, m, N, sizeMB, allocateTime, loadH2D, calcTime, loadD2H, freeTime, totalTime, cpuTime); - } - return 0; -} \ No newline at end of file diff --git a/cuda-performance/src/vector_add.cu b/cuda-performance/src/vector_add.cu deleted file mode 100644 index 709b30c..0000000 --- a/cuda-performance/src/vector_add.cu +++ /dev/null @@ -1,63 +0,0 @@ - -#include "vector_add.h" - -// Addition kernel: Adds two vectors a, b and stores the result in c -// The parameter m is used to increase the computation time per thread -// by repeating the addition operation m times. This is useful for performance -// testing with different computation loads. -template -__global__ void vectorAdd( - const T* __restrict__ d_a, - const T* __restrict__ d_b, - T* __restrict__ d_c, - int N, - int m=1) { - - int idx = blockDim.x * blockIdx.x + threadIdx.x; - if (idx >= N) return; - T s = d_a[idx] + d_b[idx]; - T acc = T(0); - for (int j = 0; j < m -1; ++j) { - acc = acc + s; - } - d_c[idx] = acc; -} - -// Wrapper to the cuda kernel call -template -void vectorAdd_wrapper( - const T* __restrict__ d_a, - const T* __restrict__ d_b, - T* __restrict__ d_c, - int N, - int m, - int ThreadsPerBlock){ - - int blocksPerGrid = (N + ThreadsPerBlock - 1) / ThreadsPerBlock; - vectorAdd<<>>(d_a, d_b, d_c, N, m); - CHECK_LAST_CUDA_ERROR(); - CHECK_CUDA_ERROR(cudaDeviceSynchronize()); -} - -template void vectorAdd_wrapper(const float*, const float*, float*, int, int, int); -template void vectorAdd_wrapper(const double*, const double*, double*, int, int, int); - - -// Function to perform a dummy calculation to warm up the GPU -// This is only used to boot up the GPU and avoid measuring initialization overhead in timing -template -void dummyCalculation(){ - T *dummy_a, *dummy_b, *dummy_c; - CHECK_CUDA_ERROR(cudaMalloc(&dummy_a, sizeof(T))); - CHECK_CUDA_ERROR(cudaMalloc(&dummy_b, sizeof(T))); - CHECK_CUDA_ERROR(cudaMalloc(&dummy_c, sizeof(T))); - vectorAdd<<<1, 1>>>(dummy_a, dummy_b, dummy_c, 1); - CHECK_LAST_CUDA_ERROR(); - CHECK_CUDA_ERROR(cudaFree(dummy_a)); - CHECK_CUDA_ERROR(cudaFree(dummy_b)); - CHECK_CUDA_ERROR(cudaFree(dummy_c)); -} - -template void dummyCalculation(); -template void dummyCalculation(); -template void dummyCalculation(); \ No newline at end of file From 30ad28de4a3b6bfecae2ec492846d9bcae0530f4 Mon Sep 17 00:00:00 2001 From: Sebas <12949040+SebastiaAgramunt@users.noreply.github.com> Date: Thu, 9 Jul 2026 14:48:37 +0000 Subject: [PATCH 20/21] gitignore --- .gitignore | 3 --- 1 file changed, 3 deletions(-) diff --git a/.gitignore b/.gitignore index 28b1a29..0d66ca8 100644 --- a/.gitignore +++ b/.gitignore @@ -188,6 +188,3 @@ external # results and analysis *.deb - -# untracked local project -cuda-performance/ From 7ed4d7b89084bc4406e3fe9450ee8896d6bd927b Mon Sep 17 00:00:00 2001 From: Sebas <12949040+SebastiaAgramunt@users.noreply.github.com> Date: Sat, 11 Jul 2026 09:17:37 +0000 Subject: [PATCH 21/21] modif --- vLLM-profiling/README.md | 36 ++++++++++++++++++++++++++++++++++++ vLLM-profiling/benchmark.sh | 12 ++++++++++-- vLLM-profiling/install.sh | 17 ++--------------- 3 files changed, 48 insertions(+), 17 deletions(-) create mode 100644 vLLM-profiling/README.md diff --git a/vLLM-profiling/README.md b/vLLM-profiling/README.md new file mode 100644 index 0000000..3a1a259 --- /dev/null +++ b/vLLM-profiling/README.md @@ -0,0 +1,36 @@ +# vLLM Profiling + +Profiles a vLLM batch-inference run with Nsight Systems to inspect CUDA kernel activity, NVTX ranges, and OS-level events during generation. + +## What it does + +`bench.py` loads `meta-llama/Meta-Llama-3-8B-Instruct` with vLLM and runs a batch of 32 prompts (4 topics x 8 repeats) through `llm.generate`. `benchmark.sh` wraps that script with `nsys profile`, capturing a trace you can open in the Nsight Systems GUI. + +## Files + +- `install.sh` — creates a `.venv`, installs `vllm` + `huggingface_hub`, and checks the GPU driver via `nvidia-smi`. Nsight Systems CLI installation is included but commented out (assumes it's already installed at `/opt/nvidia/nsight-systems-cli`). +- `benchmark.sh` — activates the venv, loads `.env`, and runs `bench.py` under `nsys profile` (trace: `cuda,nvtx,osrt`, with CUDA graph node tracing). +- `bench.py` — the actual vLLM workload being profiled. +- `.env` — local secrets (gitignored). Must define `HF_TOKEN` for the gated Llama 3 model. +- `*.nsys-rep` — Nsight Systems trace output (gitignored). + +## Prerequisites + +- NVIDIA GPU + driver (verified by `nvidia-smi` in `install.sh`) +- Nsight Systems CLI (`nsys`) installed and on `PATH`, or available at `/opt/nvidia/nsight-systems-cli` +- A Hugging Face token with access to `meta-llama/Meta-Llama-3-8B-Instruct` + +## Usage + +```bash +# one-time setup +./install.sh + +# create .env with your HF token +echo 'export HF_TOKEN=hf_...' > .env + +# run the profiled benchmark +./benchmark.sh +``` + +This produces `vllm_llama3_8b_bs32.nsys-rep`, which can be opened with `nsys-ui` or the Nsight Systems desktop app. diff --git a/vLLM-profiling/benchmark.sh b/vLLM-profiling/benchmark.sh index 8feb026..0ca8e73 100755 --- a/vLLM-profiling/benchmark.sh +++ b/vLLM-profiling/benchmark.sh @@ -10,8 +10,16 @@ source ${THIS_DIR}/.venv/bin/activate python ${THIS_DIR}/bench.py -# now the real profiled run -${NSIGHT_SYSTEMS_CLI_PATH}/${CAPPED_NSIGHT_SYS_VERSION}/bin/nsys profile \ +# # now the real profiled run +# ${NSIGHT_SYSTEMS_CLI_PATH}/${CAPPED_NSIGHT_SYS_VERSION}/bin/nsys profile \ +# --trace=cuda,nvtx,osrt \ +# --cuda-graph-trace=node \ +# --output=vllm_llama3_8b_bs32 \ +# --force-overwrite=true \ +# python ${THIS_DIR}/bench.py + + +nsys profile \ --trace=cuda,nvtx,osrt \ --cuda-graph-trace=node \ --output=vllm_llama3_8b_bs32 \ diff --git a/vLLM-profiling/install.sh b/vLLM-profiling/install.sh index 2d0c72d..f6cfb8a 100755 --- a/vLLM-profiling/install.sh +++ b/vLLM-profiling/install.sh @@ -1,9 +1,6 @@ #!/bin/bash THIS_DIR=$(dirname "$(realpath "$0")") -NSIGHT_SYS_VERSION=2026.3.1.157-3804839 -CAPPED_NSIGHT_SYS_VERSION=$(echo ${NSIGHT_SYS_VERSION} | cut -d. -f1-3) -NSIGHT_SYSTEMS_CLI_PATH=/opt/nvidia/nsight-systems-cli # create python virtual env rm -rf ${THIS_DIR}/.venv @@ -16,15 +13,5 @@ ${THIS_DIR}/.venv/bin/python -m pip install vllm huggingface_hub # Confirm the driver itself is healthy before touching packages nvidia-smi || { echo "GPU/driver not responding — stop here, don't proceed"; exit 1; } -wget https://developer.nvidia.com/downloads/assets/tools/secure/nsight-systems/2026_3/NsightSystems-linux-cli-public-${NSIGHT_SYS_VERSION}.deb -O ${THIS_DIR}/nsight-systems.deb -sudo apt install ${THIS_DIR}/nsight-systems.deb - -# # What did dpkg actually complain about? (rerun to see the real error, not swallowed by the script) -sudo dpkg -i ${THIS_DIR}/nsight-systems.deb - -# # If that reports missing deps, resolve them narrowly: -sudo apt-get install -f -y --no-install-recommends - -${NSIGHT_SYSTEMS_CLI_PATH}/${CAPPED_NSIGHT_SYS_VERSION}/bin/nsys --version - - +# confirm nsys is installed and working +nsys --version || { echo "nsys not found or not working — stop here, don't proceed"; exit 1; }