[REFAC]: Adding changelog + bench-scripts

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](docs/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 <qblas/qblas.h>`. 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](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](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](docs/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.

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
 

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 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 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.
+  - 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}")

performance_bottlenecks.md → 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.