From add6c6e9bfc82fc5f68d4d8514468f477f302a18 Mon Sep 17 00:00:00 2001 From: SwayamInSync Date: Thu, 14 May 2026 17:24:03 +0000 Subject: [PATCH 1/2] post-release: check in bench harness + CHANGELOG.md bench/bench_quad_vs_numpy.py is the harness that produced the quad-vs-f64 numbers in perf_comparison_with_old.md. It now lives in tree so the perf claims are reproducible from a clean checkout. CHANGELOG.md documents the 1.5.0 release (highlights, added/changed/ removed, known limitations) and seeds the format for future versions. --- CHANGELOG.md | 112 ++++++++++++++++++++++ bench/bench_quad_vs_numpy.py | 181 +++++++++++++++++++++++++++++++++++ 2 files changed, 293 insertions(+) create mode 100644 CHANGELOG.md create mode 100644 bench/bench_quad_vs_numpy.py diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 index 0000000..5c453de --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,112 @@ +# Changelog + +All notable changes to QBLAS are documented in this file. Format follows +[Keep a Changelog](https://keepachangelog.com/en/1.1.0/); QBLAS adheres to +[Semantic Versioning](https://semver.org/spec/v2.0.0.html). + +--- + +## [1.5.0] — 2026-05-14 + +Full rewrite of QBLAS from the legacy header-only C++ template +implementation (1.0) to a compiled shared library with a stable +CBLAS-style C ABI, runtime CPU dispatch, and OpenMP-backed parallelism. + +### Highlights + +Measured on AMD EPYC 7V13 (Zen 3, AVX2 tier), same f64 baseline +(scipy-openblas Haswell tier), same harness on both versions: + +| metric | 1.0 | 1.5 | +|----------------------------------------------|----------|----------| +| single-thread gemm slowdown vs OpenBLAS | ~1 400× | ~800× | +| 16-thread gemm @ n=256 slowdown vs OpenBLAS | 7 072× | 287× | +| single-core kernel speedup | baseline | 1.6-1.7× | +| threading at small N (≤256) | broken* | working | + +\* Old code used a fixed `nc = NC_DEFAULT = 512`, so 128×128 gemm at +16 threads ran on a single thread. New code auto-scales `nc` so each +thread gets ≥ 2 blocks. + +Full before/after tables in +[perf_comparison_with_old.md](perf_comparison_with_old.md); +reproduction harness at [bench/bench_quad_vs_numpy.py](bench/bench_quad_vs_numpy.py). + +### Added + +- CBLAS-style C ABI: `cblas_q{dot,nrm2,asum,axpy,scal,gemv,ger,gemm,syrk,trmm,trsm}` + in `include/qblas/qblas.h`. Mirrors the standard CBLAS surface for the + binary128 `Sleef_quad` type, including transpose / layout / conj flags + (previously rejected by old `QuadBLAS::gemm`). +- Per-ISA OBJECT libraries: `qblas_kernels_{generic,sse2,avx2,avx512,neon}`, + all generated from one parametrised template + ([src/kernels/kernels_template.h](src/kernels/kernels_template.h)). +- Runtime CPU dispatch via CPUID at library init; override with + `QBLAS_DISPATCH={generic,sse2,avx2,avx512,neon}`. Falls back safely if + asked to enable a tier the host doesn't support. +- Goto/van-de-Geijn 5-loop blocked GEMM with dynamic `mc / kc / nc` sized + from detected L1 / L2 / L3 cache. +- Blocked TRSM and TRMM with parallelised diagonal solve (5-15× over the + reference path at all sizes). +- OpenMP threading with thresholds detected at init time (measures actual + parallel-region overhead rather than hardcoding 8 192). +- meson.build for use as a Meson subproject by downstream packages + (numpy-quaddtype integration verified on branch `qblas-rewrite-integration`). +- Numpy-vs-QBLAS correctness test (85 cases, rounded-to-double comparison). +- Google Benchmark suite (`qblas_bench`, `qblas_bench_compare`). +- Cross-dtype throughput harness ([bench/bench_quad_vs_numpy.py](bench/bench_quad_vs_numpy.py)) + for honest reproduction of the perf-comparison numbers. + +### Changed + +- **Breaking**: C++ template API removed. Consumers now link against + `libqblas` and `#include `. The old + `QuadBLAS::gemm`/`dot`/`gemv` C++ free-function surface is gone. +- **Breaking**: function signatures take `Sleef_quad` by value (or by + pointer in shim layers for ctypes/libffi callers — see + `tests/python/qblas_shim.c` for the pattern). +- `qblas_set_num_threads()` no longer mutates global OpenMP state; it + only caps QBLAS-internal scheduling. Callers that want process-wide + thread-count changes should use `OMP_NUM_THREADS` / `omp_set_num_threads`. + +### Removed + +- Orphaned `third_party/OpenBLAS` submodule (was never referenced in any + build). Drops ~200 MB from `git clone --recurse-submodules`. + +### Known limitations + +- **Windows / MSVC is not supported.** The kernel template uses GCC-only + flags (`-march`, `-mavx2`), POSIX-only APIs (`clock_gettime`, `sysconf`, + `_SC_LEVEL2_CACHE_SIZE`), and GCC built-ins for CPUID. Downstream + packages can gate QBLAS off via their build system (numpy-quaddtype + uses `-Ddisable_quadblas=true` automatically on Windows). +- **AVX-512 tier is built but never run on hardware in CI** — the bench + host is Zen 3 (no AVX-512). The dispatcher selects it correctly on + Sapphire Rapids / Zen 4, but those paths are uncovered. Tracking: + [performance_bottlenecks.md §3](performance_bottlenecks.md). +- **Single-thread quad-FMA throughput is capped at SLEEF's TD-FMA + ceiling** (~60 cycles per FMA on Zen 3, ~800× the f64 cost). Not a + blocker for 1.5; the optimisation plan to bring this to ~100× lives in + [performance_bottlenecks.md §1](performance_bottlenecks.md). +- **Level-2 `qger` / `qsymv` / `qtrsv` still drop through to a + scalar+dispatched-axpy path** (qgemv is fully vectorised + threaded). + [performance_bottlenecks.md §4](performance_bottlenecks.md). + +### CI + +Matrix: `ubuntu-latest` × {gcc, clang} × `macos-14` × `macos-15`. +Per-tier dispatch test on Linux (`QBLAS_DISPATCH={generic,sse2,avx2}`). +Numpy comparison test on every push. Benchmark artifact uploaded by the +`bench-baseline` job. + +--- + +## [1.0.0] — pre-2026 + +The original header-only C++ template implementation. Compiled into the +consumer at use-site, x86-64 path hardcoded to SSE2, no transpose +support in `gemm`, no runtime dispatch. + +This entry is here for reference only; the 1.0 line is no longer +maintained. diff --git a/bench/bench_quad_vs_numpy.py b/bench/bench_quad_vs_numpy.py new file mode 100644 index 0000000..5730bc1 --- /dev/null +++ b/bench/bench_quad_vs_numpy.py @@ -0,0 +1,181 @@ +"""End-to-end throughput comparison: numpy float64 (scipy-openblas) vs +numpy-quaddtype (this QBLAS) for the BLAS-1/2/3 hot path. + +This is the harness that produced the numbers in perf_comparison_with_old.md +and in the QBLAS 1.5.0 release notes. Reproducing those numbers from a +clean checkout should give you the same shape (single-thread quad-FMA cost +~800x f64; multi-thread closes that gap to ~50-550x depending on size and +thread count). Exact values depend on host CPU, memory bandwidth, and +OpenBLAS tier. + +Methodology (matching what perf_comparison_with_old.md describes): + - Pin both libraries to BENCH_THREADS via OMP_NUM_THREADS / OPENBLAS_NUM_THREADS + before importing numpy. Same thread count on both sides; no asymmetric + advantage. + - For each (op, dtype, N): generate arrays once, warm up, then time an + inner loop that is auto-calibrated so each timed sample takes >= 50 ms + (or 1 iteration if a single call already exceeds 50 ms). + - Repeat the timed sample N_REPEATS times. Report median (robust to GC + and OS-noise outliers) and stdev/median as a stability indicator. Trim + the fastest + slowest sample (Olympic trimming) when N_REPEATS >= 5. + - Same arrays used for both libraries: no setup-cost asymmetry. + +Requirements: + pip install numpy-quaddtype (which depends on numpy) + +Usage: + BENCH_THREADS=1 python bench/bench_quad_vs_numpy.py + BENCH_THREADS=16 python bench/bench_quad_vs_numpy.py + +Output paths can be redirected with BENCH_OUTPUT_DIR (default: $PWD). +""" +import os +import sys +import gc +import time +import json +from statistics import median, stdev + +THREADS = int(os.environ.get("BENCH_THREADS", "16")) +os.environ.setdefault("OMP_NUM_THREADS", str(THREADS)) +os.environ.setdefault("OPENBLAS_NUM_THREADS", str(THREADS)) +os.environ.setdefault("MKL_NUM_THREADS", str(THREADS)) +os.environ.setdefault("OMP_PROC_BIND", "close") +os.environ.setdefault("OMP_PLACES", "cores") + +OUT_DIR = os.environ.get("BENCH_OUTPUT_DIR", os.getcwd()) + +import numpy as np +try: + import numpy_quaddtype as nq +except ImportError: + sys.exit("error: this harness requires numpy-quaddtype; install with " + "`pip install numpy-quaddtype` or build it from source against this QBLAS.") + +QPRC = nq.QuadPrecDType() +RNG = np.random.default_rng(0) + +TARGET_SAMPLE_S = 0.05 # auto-calibrate inner loop to >= this many seconds +N_WARMUP = 3 +N_REPEATS = 9 # odd so median is well-defined; 9 -> trim top/bottom 2 + + +def timed(fn, *args, target=TARGET_SAMPLE_S): + """Return (per_call_seconds, stdev_seconds, n_inner, n_samples).""" + # Calibrate inner loop count so each timed sample is >= `target` seconds. + inner = 1 + while True: + gc.collect(); gc.disable() + t0 = time.perf_counter() + for _ in range(inner): + fn(*args) + dt = time.perf_counter() - t0 + gc.enable() + if dt >= target or inner >= 2048: + break + inner = max(2 * inner, int(inner * target / max(dt, 1e-9) * 1.2)) + + for _ in range(N_WARMUP): + for _ in range(inner): + fn(*args) + + times = [] + for _ in range(N_REPEATS): + gc.collect(); gc.disable() + t0 = time.perf_counter() + for _ in range(inner): + fn(*args) + times.append((time.perf_counter() - t0) / inner) + gc.enable() + + if len(times) >= 5: + times = sorted(times)[1:-1] + med = median(times) + sd = stdev(times) if len(times) > 1 else 0.0 + return med, sd, inner, len(times) + + +def gops(fmas, t): return fmas / t / 1e9 + +def fmt_t(t): + if t < 1e-3: return f"{t*1e6:7.1f} us" + if t < 1: return f"{t*1e3:7.2f} ms" + return f"{t:7.3f} s " + +def run_op(name, ops_per_call, dtype_label, fn, *args): + med, sd, inner, n = timed(fn, *args) + rate = gops(ops_per_call, med) + rel = sd / med if med > 0 else 0 + unit = "GFLOPS" if dtype_label == "f64" else "GFMA/s" + print(f" {name:<22} {fmt_t(med)} {rate:8.3f} {unit:<8} " + f"+/-{rel*100:5.2f}% (inner={inner}, n={n})") + return med + + +print(f"# Threads: {THREADS}") +print(f"# numpy: {np.__version__}") +print(f"# numpy_quaddtype: {nq.get_quadblas_version()}") +print(f"# QBLAS threads: {nq.get_num_threads()}") +print(f"# Methodology: {N_WARMUP} warmup + {N_REPEATS} timed samples (median, trim hi/lo);") +print(f"# inner loop auto-calibrated to >= {int(TARGET_SAMPLE_S*1000)} ms per sample\n") + +DOT_SIZES = [1 << 16, 1 << 18, 1 << 20, 1 << 22] +GEMV_SIZES = [512, 1024, 2048] +# Skip n=1024 quad gemm at single thread (~4.6s/call x 9 samples = too slow). +GEMM_SIZES = [128, 256, 512, 1024] if THREADS > 1 else [128, 256, 512] + +results = {} + +print("== BLAS-1: dot (vector-vector inner product) ==") +for n in DOT_SIZES: + xd = RNG.random(n); yd = RNG.random(n) + xq = xd.astype(QPRC); yq = yd.astype(QPRC) + fmas = n + print(f" N = {n:>7}") + td = run_op("f64", fmas, "f64", np.matmul, xd, yd) + tq = run_op("quad", fmas, "quad", np.matmul, xq, yq) + results[("dot", n)] = (td, tq) + del xd, yd, xq, yq + +print("\n== BLAS-2: gemv (matrix-vector product) ==") +for n in GEMV_SIZES: + Ad = RNG.random((n, n)); xd = RNG.random(n) + Aq = Ad.astype(QPRC); xq = xd.astype(QPRC) + fmas = 2 * n * n + print(f" N = {n:>5}") + td = run_op("f64", fmas, "f64", np.matmul, Ad, xd) + tq = run_op("quad", fmas, "quad", np.matmul, Aq, xq) + results[("gemv", n)] = (td, tq) + del Ad, xd, Aq, xq + +print("\n== BLAS-3: gemm (matrix-matrix product) ==") +for n in GEMM_SIZES: + Ad = RNG.random((n, n)); Bd = RNG.random((n, n)) + Aq = Ad.astype(QPRC); Bq = Bd.astype(QPRC) + fmas = 2 * n * n * n + print(f" N = {n:>5}") + td = run_op("f64", fmas, "f64", np.matmul, Ad, Bd) + tq = run_op("quad", fmas, "quad", np.matmul, Aq, Bq) + results[("gemm", n)] = (td, tq) + del Ad, Bd, Aq, Bq + +print("\n" + "="*78) +print(f"SUMMARY (threads={THREADS})") +print("="*78) +print(f"{'op':<5} {'N':>6} {'f64 time':>11} {'f64 rate':>13} {'quad time':>11} {'quad rate':>13} slowdown") +print("-"*78) +for (op, n), (td, tq) in results.items(): + fmas = {"dot": n, "gemv": 2*n*n, "gemm": 2*n*n*n}[op] + print(f"{op:<5} {n:>6} {fmt_t(td):>11} {gops(fmas, td):>9.2f} GFLOPS" + f" {fmt_t(tq):>11} {gops(fmas, tq):>8.3f} GFMA/s {tq/td:>6.0f}x") + +out_path = os.path.join(OUT_DIR, f"bench_quad_vs_numpy_t{THREADS}.json") +with open(out_path, "w") as f: + json.dump({ + "threads": THREADS, + "numpy_version": np.__version__, + "qblas_version": nq.get_quadblas_version(), + "results": {f"{op}/n={n}": {"f64_s": td, "quad_s": tq} + for (op, n), (td, tq) in results.items()}, + }, f, indent=2) +print(f"\nWrote {out_path}") From ad3659eda33558eb5146fa550d7ec139bb23604c Mon Sep 17 00:00:00 2001 From: SwayamInSync Date: Thu, 14 May 2026 17:32:44 +0000 Subject: [PATCH 2/2] move perf docs into docs/ to keep top-level tree clean MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Top-level used to have five .md files (README, CHANGELOG, CLAUDE, perf_comparison_with_old, performance_bottlenecks). Reference material that the README links to belongs deeper. Moves both perf docs into docs/; updates the cross-references in README, CHANGELOG, the bench harness's docstring, and the internal repro-instructions inside performance_bottlenecks.md. Top level is now README + CHANGELOG + CLAUDE only — the conventional three for an OSS project. --- CHANGELOG.md | 8 ++++---- README.md | 4 ++-- bench/bench_quad_vs_numpy.py | 4 ++-- .../perf_comparison_with_old.md | 0 .../performance_bottlenecks.md | 2 +- 5 files changed, 9 insertions(+), 9 deletions(-) rename perf_comparison_with_old.md => docs/perf_comparison_with_old.md (100%) rename performance_bottlenecks.md => docs/performance_bottlenecks.md (99%) diff --git a/CHANGELOG.md b/CHANGELOG.md index 5c453de..8d4fb53 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -29,7 +29,7 @@ Measured on AMD EPYC 7V13 (Zen 3, AVX2 tier), same f64 baseline thread gets ≥ 2 blocks. Full before/after tables in -[perf_comparison_with_old.md](perf_comparison_with_old.md); +[perf_comparison_with_old.md](docs/perf_comparison_with_old.md); reproduction harness at [bench/bench_quad_vs_numpy.py](bench/bench_quad_vs_numpy.py). ### Added @@ -84,14 +84,14 @@ reproduction harness at [bench/bench_quad_vs_numpy.py](bench/bench_quad_vs_numpy - **AVX-512 tier is built but never run on hardware in CI** — the bench host is Zen 3 (no AVX-512). The dispatcher selects it correctly on Sapphire Rapids / Zen 4, but those paths are uncovered. Tracking: - [performance_bottlenecks.md §3](performance_bottlenecks.md). + [performance_bottlenecks.md §3](docs/performance_bottlenecks.md). - **Single-thread quad-FMA throughput is capped at SLEEF's TD-FMA ceiling** (~60 cycles per FMA on Zen 3, ~800× the f64 cost). Not a blocker for 1.5; the optimisation plan to bring this to ~100× lives in - [performance_bottlenecks.md §1](performance_bottlenecks.md). + [performance_bottlenecks.md §1](docs/performance_bottlenecks.md). - **Level-2 `qger` / `qsymv` / `qtrsv` still drop through to a scalar+dispatched-axpy path** (qgemv is fully vectorised + threaded). - [performance_bottlenecks.md §4](performance_bottlenecks.md). + [performance_bottlenecks.md §4](docs/performance_bottlenecks.md). ### CI diff --git a/README.md b/README.md index 2cdc621..4be7e80 100644 --- a/README.md +++ b/README.md @@ -254,9 +254,9 @@ On AMD EPYC 7V13 (96 cores, AVX2 tier, OMP_NUM_THREADS=16, median of 3): | qtrmm | 512 x 512 | 350.3 | | qtrsm | 512 x 512 | 335.7 | -For full before/after tables (including single-thread and 96-thread numbers), see [perf_comparison_with_old.md](perf_comparison_with_old.md). +For full before/after tables (including single-thread and 96-thread numbers), see [docs/perf_comparison_with_old.md](docs/perf_comparison_with_old.md). -For remaining performance opportunities not yet taken, see [performance_bottlenecks.md](performance_bottlenecks.md). +For remaining performance opportunities not yet taken, see [docs/performance_bottlenecks.md](docs/performance_bottlenecks.md). ## Repository layout diff --git a/bench/bench_quad_vs_numpy.py b/bench/bench_quad_vs_numpy.py index 5730bc1..fea438a 100644 --- a/bench/bench_quad_vs_numpy.py +++ b/bench/bench_quad_vs_numpy.py @@ -1,14 +1,14 @@ """End-to-end throughput comparison: numpy float64 (scipy-openblas) vs numpy-quaddtype (this QBLAS) for the BLAS-1/2/3 hot path. -This is the harness that produced the numbers in perf_comparison_with_old.md +This is the harness that produced the numbers in docs/perf_comparison_with_old.md and in the QBLAS 1.5.0 release notes. Reproducing those numbers from a clean checkout should give you the same shape (single-thread quad-FMA cost ~800x f64; multi-thread closes that gap to ~50-550x depending on size and thread count). Exact values depend on host CPU, memory bandwidth, and OpenBLAS tier. -Methodology (matching what perf_comparison_with_old.md describes): +Methodology (matching what docs/perf_comparison_with_old.md describes): - Pin both libraries to BENCH_THREADS via OMP_NUM_THREADS / OPENBLAS_NUM_THREADS before importing numpy. Same thread count on both sides; no asymmetric advantage. diff --git a/perf_comparison_with_old.md b/docs/perf_comparison_with_old.md similarity index 100% rename from perf_comparison_with_old.md rename to docs/perf_comparison_with_old.md diff --git a/performance_bottlenecks.md b/docs/performance_bottlenecks.md similarity index 99% rename from performance_bottlenecks.md rename to docs/performance_bottlenecks.md index aa1c8c0..9d0820d 100644 --- a/performance_bottlenecks.md +++ b/docs/performance_bottlenecks.md @@ -199,7 +199,7 @@ ctest # correctness MUST stay green OMP_NUM_THREADS=16 ./bench/qblas_bench_compare \ --benchmark_min_time=0.3s --benchmark_repetitions=3 \ --benchmark_report_aggregates_only=true -# Compare against perf_comparison_with_old.md numbers. +# Compare against perf_comparison_with_old.md numbers (in this docs/ dir). ``` Single change at a time. Keep the bench harness untouched.