diff --git a/OpenBLAS/README.md b/OpenBLAS/README.md new file mode 100644 index 0000000..8ce8dbe --- /dev/null +++ b/OpenBLAS/README.md @@ -0,0 +1,167 @@ +# OpenBLAS on RISC-V - verification & microbenchmarks + +Small, self-contained programs used to check the **performance** and +**numerical correctness** of a BLAS backend, and - when something is wrong - to +localize *which* kernel is broken. Written while bringing up the RISC-V RVV +(SpaceMiT X60 / K1) OpenBLAS build, which surfaced two distinct kernel bugs: a +`gemv_n` vector kernel that silently returned `NaN`, and a set of `_rvv_v1` +TRSM kernels that were not VLEN-agnostic and corrupted results on non-`x280` +targets. + +| File | What it does | +|---|---| +| `bench_dgemm.c` | Times square `C = A*B` (level-3 `dgemm`, 3 reps), reports GFLOP/s, prints `C[0]` so garbage is visible. | +| `difftest.c` | `dlopen`s a BLAS `.so` and runs level-1/2/3 routines (incl. `dgemv`, `dgemm`, `dtrsm`) on fixed inputs, printing `sum / sumsq / nan / inf` per routine for differential comparison. | +| `verify_ctrsm.c` | Full-parameter TRSM correctness sweep ({Left,Right} x {Upper,Lower} x {No,T,ConjT} x {Non,Unit} x size grid), recomputes `op(A)*X` and checks max residual. Shown for the complex-single (`ctrsm`) case; S/D/Z are structurally identical (swap type + `cblas_?trsm` + element generator). | + +All are single C files, MIT-licensed. `bench_dgemm` / `difftest` select the +backend at runtime via FlexiBLAS (`FLEXIBLAS=...`) or OpenBLAS +(`OPENBLAS_CORETYPE`, `OPENBLAS_NUM_THREADS`) - no recompilation needed to A/B +two backends. + +## Build + +```bash +gcc -O2 bench_dgemm.c -o bench_dgemm -lflexiblas # or -lopenblas +gcc -O2 difftest.c -o difftest -ldl -lm +gcc -O2 verify_ctrsm.c -o verify_ctrsm -I -lopenblas -lm +``` + +## Run + +```bash +# performance (pick threads + backend via env) +OPENBLAS_NUM_THREADS=8 ./bench_dgemm 4096 + +# correctness of a specific backend .so +./difftest /path/to/libopenblas.so + +# TRSM correctness sweep (exit 0 = all pass) +./verify_ctrsm +``` + +## A/B a backend with FlexiBLAS (no rebuild) + +```bash +# scalar baseline +OPENBLAS_CORETYPE=RISCV64_GENERIC ./bench_dgemm 2048 +OPENBLAS_CORETYPE=RISCV64_GENERIC ./difftest /path/to/libopenblas.so +# vector / alternative backend +FLEXIBLAS=/path/to/libopenblas.so ./bench_dgemm 2048 +./difftest /path/to/vector-libopenblas.so +``` + +## Example - the OpenBLAS 0.3.30 RISC-V `gemv_n` NaN bug (SpaceMiT X60) + +`difftest` run against the **stock** EESSI-CVMFS OpenBLAS 0.3.30 (a +`DYNAMIC_ARCH` build that dispatches the `ZVL256B` RVV kernels on the X60), the +same lib forced to scalar, and a **patched** RVV build: + +| backend / dispatch | `dgemv` nan | `dgemm` nan | `dtrsm` nan | `dgemv` sum | +|---|--:|--:|--:|---| +| stock CVMFS, **default RVV** (`ZVL256B`) | **192** | 0 | 0 | 198.94 (wrong) | +| stock CVMFS, forced scalar (`RISCV64_GENERIC`) | 0 | 0 | 0 | 42.06549 (reference) | +| **patched RVV** (`gemv_n` fix) | 0 | 0 | 0 | 42.06549 (matches reference) | + +The differential test pins the fault to `dgemv` alone - `dgemm` and `dtrsm` are +correct even on the broken RVV build. That is exactly why a plain `dgemm` +benchmark looks fine while **HPL fails with a NaN residual**: HPL's panel +factorization leans on `dgemv`. The patched vector build reproduces the scalar +reference bit-for-bit. + +Performance on the same machine (`bench_dgemm`, 1 core, N=2048, via FlexiBLAS): + +| backend | GFLOP/s | `C[0]` | +|---|--:|--:| +| scalar (`RISCV64_GENERIC`) | 1.16 | 245.24 | +| patched RVV (`ZVL256B`) | 2.62 | 245.24 | + +i.e. ~2.3x faster and numerically identical. (Root cause: OpenBLAS 0.3.30 +`kernel/riscv64/gemv_n_vector.c` zeroes an *uninitialized* vector accumulator; +fixed upstream after 0.3.30.) + +## Cross-board confirmation - Banana Pi BPI-F3 (same K1 / X60 SoC) + +The [Banana Pi BPI-F3](https://www.banana-pi.org/) uses the same SpaceMiT K1 +(8x X60 @ 1.6 GHz, RVV 1.0 VLEN=256) as the Orange Pi RV2 above, so the same +bug and fix reproduce on a second board - here on EESSI `2025.06-001`, patched +OpenBLAS 0.3.30 ([easyconfigs #26444](https://github.com/easybuilders/easybuild-easyconfigs/pull/26444)) +via FlexiBLAS. + +`difftest` - bit-identical to the RV2: + +| backend / dispatch | `dgemv` nan | `dgemm` nan | `dtrsm` nan | `dgemv` sum | +|---|--:|--:|--:|---| +| stock CVMFS, **default RVV** (`ZVL256B`) | **192** | 0 | 0 | 198.94 (wrong) | +| stock CVMFS, forced scalar (`RISCV64_GENERIC`) | 0 | 0 | 0 | 42.06549 (reference) | +| **patched RVV** (`gemv_n` fix) | 0 | 0 | 0 | 42.06549 (matches reference) | + +`bench_dgemm` (1 core, N=2048, via FlexiBLAS): + +| backend | GFLOP/s | `C[0]` | +|---|--:|--:| +| scalar (`RISCV64_GENERIC`) | 1.26 | 245.24 | +| patched RVV (`ZVL256B`) | 2.96 | 245.24 | + +~2.35x faster and numerically identical - confirming the RV2 result on a second +K1 board. (Threaded: **17.71 GFLOP/s** at 8 cores, N=4096.) + +## Example - the RVV `_rvv_v1` TRSM VLEN bug (OpenMathLib/OpenBLAS#5928) + +A second, independent kernel bug on the same X60 build: the RVV triangular-solve +kernels `trsm_kernel_{LN,LT,RN,RT}_rvv_v1.c` are **not VLEN-agnostic**. They tile +the packed `A` panel by the *runtime* `VSETVL_MAX` (16 fp32 elements at LMUL=2 on +the X60's VLEN=256), but the GEMM driver they link against packs `A` by the +*compile-time* `GEMM_UNROLL_M`. Those two agree only when `GEMM_UNROLL_M == 16` +(the `x280` / `ZVL256B` target). On any build where `GEMM_UNROLL_M = 8` +(`ZVL128B`, the default for a plain RVV target) the tile stride `16 != 8` +diverges from how the panel was packed, and TRSM silently returns wrong results. + +`verify_ctrsm` (and its S/D/Z siblings) run the full parameter sweep. On a +`ZVL128B` build (`GEMM_UNROLL_M = 8`): + +| build | STRSM | DTRSM | CTRSM | ZTRSM | total | +|---|--:|--:|--:|--:|--:| +| stock `_rvv_v1` (pre-fix) | fails | **408 / 1600 fail** | fails | fails | broken | +| **fixed** (`GEMM_UNROLL_M`-tiled) | 0 fail | 0 fail | 0 fail | 0 fail | **8000 / 0** | + +worst residuals on the fixed build: STRSM `1.4e-7`, DTRSM `3.3e-16`, +CTRSM `3.0e-7`, ZTRSM `5.5e-16`. The fix also passes `8000 / 0` on a `ZVL256B` +build (`GEMM_UNROLL_M = 16`), confirming it is genuinely VLEN-agnostic rather +than re-hardcoded to the other width. + +**The fix:** replace the `VSETVL_MAX` outer tiling with generic +`GEMM_UNROLL_M`-tiled outer control plus a power-of-2 `M`-remainder path, keeping +the vectorized RVV `solve()` inner kernel unchanged. + +### Performance - no regression + +`benchmark/trsm.c` (from the OpenBLAS tree), single-thread, `L/U/N/N`, +`OPENBLAS_LOOPS=20`, best-of-3 on the X60, both libs `ZVL128B`. Delta = fixed vs +original `_rvv_v1`, MFlops: + +| size | STRSM | DTRSM | CTRSM | ZTRSM | +|---|--:|--:|--:|--:| +| 8 | +4.3% | +4.3% | -3.0% | +0.3% | +| 16 | +32.6% | +0.2% | +31.9% | +2.6% | +| 32 | +40.9% | +0.4% | +29.2% | +1.6% | +| 64 | +35.6% | +1.4% | +21.4% | +1.5% | +| 128 | +26.0% | ~0% | +8.6% | -2.6% | +| 256 | +17.3% | -3.9% | +6.6% | -4.2% | +| 400 | +11.5% | -6.9% | -0.8% | -1.5% | + +- **STRSM / CTRSM**: consistent speedups (+10-40%), largest at small/mid sizes - + `GEMM_UNROLL_M` tiling reuses the packed diagonal block better than per-call + `VSETVL_MAX` retiling. +- **DTRSM / ZTRSM**: within measurement noise (+-~4%) on this loaded board; a + 6-rep re-measure at DTRSM `m=128` gave overlapping distributions (orig + 1273-1412, fixed 1257-1371; equal medians), so the scattered negatives are + jitter, not a systematic slowdown. + +> **Caveat on the baseline:** on `ZVL128B` the original `_rvv_v1` kernels are +> *numerically broken* (that is the bug), so their timings are not a legitimate +> speed baseline - you cannot "regress" against a kernel that returns garbage. +> The fix makes `ZVL128B` correct for the first time while being at least as +> fast, usually faster. On `ZVL256B` / `x280` (`VSETVL_MAX == GEMM_UNROLL_M == +> 16`) the tile width is unchanged, so the rewrite is a structural no-op for +> throughput there. diff --git a/dgemm/bench_dgemm.c b/OpenBLAS/bench_dgemm.c similarity index 100% rename from dgemm/bench_dgemm.c rename to OpenBLAS/bench_dgemm.c diff --git a/dgemm/difftest.c b/OpenBLAS/difftest.c similarity index 100% rename from dgemm/difftest.c rename to OpenBLAS/difftest.c diff --git a/OpenBLAS/verify_ctrsm.c b/OpenBLAS/verify_ctrsm.c new file mode 100644 index 0000000..aeb7d97 --- /dev/null +++ b/OpenBLAS/verify_ctrsm.c @@ -0,0 +1,101 @@ +/* + * verify_ctrsm.c - CTRSM (complex single) numerical-correctness sweep. + * + * Solves op(A) X = alpha B (Left) or X op(A) = alpha B (Right) for a triangular + * A across the full parameter space - {Left,Right} x {Upper,Lower} x + * {NoTrans,Trans,ConjTrans} x {NonUnit,Unit} x a size grid - then recomputes + * op(A)*X (or X*op(A)) and checks the max residual against the original B. + * Used to validate the RISC-V RVV (_rvv_v1) TRSM kernels: the VLEN-agnostic fix + * in OpenMathLib/OpenBLAS#5928 passes 0 fails here on both ZVL128B (GEMM_UNROLL_M + * =8) and ZVL256B (=16); the pre-fix kernels fail on ZVL128B (unroll 8 != the + * VSETVL_MAX=16 the kernel tiled packed-A by). + * + * The S/D/Z variants are structurally identical - swap the type, the cblas_?trsm + * call, and the element generator (real vs complex). + * + * Build: gcc -O2 verify_ctrsm.c -o verify_ctrsm -I -lopenblas -lm + * Run: ./verify_ctrsm # exit 0 = all pass, prints "N cases, F fails" + * SPDX-License-Identifier: MIT + */ +#include +#include +#include +#include +#include "cblas.h" + +typedef float _Complex cd; +static double frand(void){ return (double)rand()/RAND_MAX - 0.5; } +static cd crand(void){ return frand() + frand()*I; } + +static void make_tri(cd *A, int n, int uplo_upper, int unit){ + for(int j=0;jj)) v = 0.3*crand(); + else v = 0.0; + A[i+j*n]=v; + } +} + +static cd opA(int trans,int uplo,int is_unit,const cd*A,int lda,int r,int c){ + /* returns element (r,c) of op(A) given stored triangular A */ + int sr,sc; + if(trans==CblasNoTrans){ sr=r; sc=c; } + else if(trans==CblasTrans){ sr=c; sc=r; } + else { sr=c; sc=r; } /* ConjTrans handled by caller via conj */ + cd a; + if(sr==sc) a = is_unit ? (1.0+0.0*I) : A[sr+sc*lda]; + else if((uplo==CblasUpper && srsc)) a=A[sr+sc*lda]; + else a=0.0; + if(trans==CblasConjTrans) a = conj(a); + return a; +} + +static double resid(int side,int uplo,int trans,int is_unit,int m,int n, + const cd*A,int lda,const cd*X,const cd*B0,int ldb){ + double maxr=0.0; + if(side==CblasLeft){ + for(int i=0;imaxr) maxr=d; + } + } else { + for(int i=0;imaxr) maxr=d; + } + } + return maxr; +} + +int main(void){ + int grid[]={1,2,4,7,8,9,10,12,15,23}; + int ng=sizeof(grid)/sizeof(grid[0]); + int sides[]={CblasLeft,CblasRight}, uplos[]={CblasUpper,CblasLower}; + int transs[]={CblasNoTrans,CblasTrans,CblasConjTrans}, diags[]={CblasNonUnit,CblasUnit}; + int fails=0,total=0; double worst=0.0; srand(999); + cd alpha=1.0+0.0*I; + for(int si=0;si<2;si++) for(int ui=0;ui<2;ui++) for(int ti=0;ti<3;ti++) for(int di=0;di<2;di++) + for(int im=0;imworst) worst=r; + if(!(r<=tol)||isnan(r)){ fails++; + if(fails<=30) + printf("FAIL side=%d uplo=%d trans=%d unit=%d m=%d n=%d resid=%g tol=%g\n", + side,uplo,trans,unit,m,n,r,tol); } + free(A);free(B);free(B0); + } + printf("=== CTRSM: %d cases, %d fails, worst_resid=%g ===\n", total, fails, worst); + return fails?1:0; +} diff --git a/README.md b/README.md index f5c8393..49a0aa4 100644 --- a/README.md +++ b/README.md @@ -20,7 +20,7 @@ correctness (finite / bit-identical results) is checked in every case. | Dir | What it measures | Axis swapped | Kind | |---|---|---|---| -| [`dgemm/`](dgemm) | BLAS DGEMM performance **and** numerical correctness; localizes a broken kernel | BLAS | microbench + verification | +| [`OpenBLAS/`](OpenBLAS) | OpenBLAS on RISC-V: DGEMM performance + differential correctness + TRSM sweep; localizes two broken RVV kernels (`gemv_n` NaN, `_rvv_v1` TRSM VLEN bug) | BLAS | microbench + verification | | [`numpy/`](numpy) | BLAS/LAPACK backend as seen through NumPy/SciPy | BLAS | application proxy | | [`hpl/`](hpl) | High-Performance Linpack, end-to-end | BLAS (FlexiBLAS) | application A/B | | [`elpa/`](elpa) | Dense real-symmetric eigensolver (ELPA, 1-stage) | BLAS | microbench | @@ -38,7 +38,7 @@ Several directories deliberately pair a **microbenchmark** with a **real-application** measurement of the *same* backend, because they often disagree — and that disagreement is the interesting result: -- **BLAS axis:** [`dgemm`](dgemm)/[`numpy`](numpy)/[`elpa`](elpa) (kernel level) +- **BLAS axis:** [`OpenBLAS`](OpenBLAS)/[`numpy`](numpy)/[`elpa`](elpa) (kernel level) → [`hpl`](hpl)/[`qe`](qe) (whole application). - **FFT axis:** [`fftw`](fftw) standalone microbench (RVV wins **1.06–1.60×**) → the same RVV FFTW dropped into a Quantum ESPRESSO SCF (**~0% end-to-end**, diff --git a/dgemm/README.md b/dgemm/README.md deleted file mode 100644 index 14f6d7c..0000000 --- a/dgemm/README.md +++ /dev/null @@ -1,98 +0,0 @@ -# BLAS / DGEMM verification - -Two small programs to check both the **performance** and the **numerical -correctness** of a BLAS backend, and - when something is wrong - to localize -*which* kernel is broken. Written while bringing up a RISC-V RVV (SpaceMiT X60 / -K1) OpenBLAS build, where a vector kernel silently returned `NaN`. - -| File | What it does | -|---|---| -| `bench_dgemm.c` | Times square `C = A*B` (level-3 `dgemm`, 3 reps), reports GFLOP/s, prints `C[0]` so garbage is visible. | -| `difftest.c` | `dlopen`s a BLAS `.so` and runs level-1/2/3 routines on fixed inputs, printing `sum / sumsq / nan / inf` per routine for differential comparison. | - -Both are single C files, MIT-licensed, and select/switch the backend at runtime -via FlexiBLAS (`FLEXIBLAS=...`) or OpenBLAS (`OPENBLAS_CORETYPE`, -`OPENBLAS_NUM_THREADS`) - no recompilation needed to A/B two backends. - -## Build - -```bash -gcc -O2 bench_dgemm.c -o bench_dgemm -lflexiblas # or -lopenblas -gcc -O2 difftest.c -o difftest -ldl -lm -``` - -## Run - -```bash -# performance (pick threads + backend via env) -OPENBLAS_NUM_THREADS=8 ./bench_dgemm 4096 - -# correctness of a specific backend .so -./difftest /path/to/libopenblas.so -``` - -## A/B a backend with FlexiBLAS (no rebuild) - -```bash -# scalar baseline -OPENBLAS_CORETYPE=RISCV64_GENERIC ./bench_dgemm 2048 -OPENBLAS_CORETYPE=RISCV64_GENERIC ./difftest /path/to/libopenblas.so -# vector / alternative backend -FLEXIBLAS=/path/to/libopenblas.so ./bench_dgemm 2048 -./difftest /path/to/vector-libopenblas.so -``` - -## Example - the OpenBLAS 0.3.30 RISC-V `gemv_n` NaN bug (SpaceMiT X60) - -`difftest` run against the **stock** EESSI-CVMFS OpenBLAS 0.3.30 (a -`DYNAMIC_ARCH` build that dispatches the `ZVL256B` RVV kernels on the X60), the -same lib forced to scalar, and a **patched** RVV build: - -| backend / dispatch | `dgemv` nan | `dgemm` nan | `dtrsm` nan | `dgemv` sum | -|---|--:|--:|--:|---| -| stock CVMFS, **default RVV** (`ZVL256B`) | **192** | 0 | 0 | 198.94 (wrong) | -| stock CVMFS, forced scalar (`RISCV64_GENERIC`) | 0 | 0 | 0 | 42.06549 (reference) | -| **patched RVV** (`gemv_n` fix) | 0 | 0 | 0 | 42.06549 (matches reference) | - -The differential test pins the fault to `dgemv` alone - `dgemm` and `dtrsm` are -correct even on the broken RVV build. That is exactly why a plain `dgemm` -benchmark looks fine while **HPL fails with a NaN residual**: HPL's panel -factorization leans on `dgemv`. The patched vector build reproduces the scalar -reference bit-for-bit. - -Performance on the same machine (`bench_dgemm`, 1 core, N=2048, via FlexiBLAS): - -| backend | GFLOP/s | `C[0]` | -|---|--:|--:| -| scalar (`RISCV64_GENERIC`) | 1.16 | 245.24 | -| patched RVV (`ZVL256B`) | 2.62 | 245.24 | - -i.e. ~2.3x faster and numerically identical. (Root cause: OpenBLAS 0.3.30 -`kernel/riscv64/gemv_n_vector.c` zeroes an *uninitialized* vector accumulator; -fixed upstream after 0.3.30.) - -## Cross-board confirmation - Banana Pi BPI-F3 (same K1 / X60 SoC) - -The [Banana Pi BPI-F3](https://www.banana-pi.org/) uses the same SpaceMiT K1 -(8x X60 @ 1.6 GHz, RVV 1.0 VLEN=256) as the Orange Pi RV2 above, so the same -bug and fix reproduce on a second board - here on EESSI `2025.06-001`, patched -OpenBLAS 0.3.30 ([easyconfigs #26444](https://github.com/easybuilders/easybuild-easyconfigs/pull/26444)) -via FlexiBLAS. - -`difftest` - bit-identical to the RV2: - -| backend / dispatch | `dgemv` nan | `dgemm` nan | `dtrsm` nan | `dgemv` sum | -|---|--:|--:|--:|---| -| stock CVMFS, **default RVV** (`ZVL256B`) | **192** | 0 | 0 | 198.94 (wrong) | -| stock CVMFS, forced scalar (`RISCV64_GENERIC`) | 0 | 0 | 0 | 42.06549 (reference) | -| **patched RVV** (`gemv_n` fix) | 0 | 0 | 0 | 42.06549 (matches reference) | - -`bench_dgemm` (1 core, N=2048, via FlexiBLAS): - -| backend | GFLOP/s | `C[0]` | -|---|--:|--:| -| scalar (`RISCV64_GENERIC`) | 1.26 | 245.24 | -| patched RVV (`ZVL256B`) | 2.96 | 245.24 | - -~2.35x faster and numerically identical - confirming the RV2 result on a second -K1 board. (Threaded: **17.71 GFLOP/s** at 8 cores, N=4096.) diff --git a/hpl/README.md b/hpl/README.md index b79bc31..f1eecac 100644 --- a/hpl/README.md +++ b/hpl/README.md @@ -53,7 +53,7 @@ the `ZVL256B` RVV kernels on the X60 - and its RVV `gemv_n` has a NaN bug, so th ``` Swapping in an OpenBLAS with the `gemv_n` fix backported gives residual -4.04e-03, **PASSED**. See [`../dgemm/`](../dgemm) for the differential test that +4.04e-03, **PASSED**. See [`../OpenBLAS/`](../OpenBLAS) for the differential test that isolated the bug to `dgemv` (HPL's panel factorization leans on `dgemv`, which is why HPL NaNs while a plain `dgemm` benchmark looks fine). diff --git a/numpy/README.md b/numpy/README.md index c2d88a0..f44c174 100644 --- a/numpy/README.md +++ b/numpy/README.md @@ -78,4 +78,4 @@ is worse than a silent NaN here: `eigvalsh` aborts with `LinAlgError: Eigenvalues did not converge` (the `gemv_n` NaN propagating through LAPACK `dsyevd`), so only scalar and the patched vector build return an answer at all. (DGEMM stays finite on stock RVV - 10.94 GFLOP/s - because only `gemv`, not -`gemm`, is affected; see [`../dgemm/`](../dgemm).) +`gemm`, is affected; see [`../OpenBLAS/`](../OpenBLAS).)