OpenBLAS RISC-V verification suite + TRSM VLEN-fix sweep

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<cblas_include> -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.

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<cblas_include> -lopenblas -lm
+ * Run:    ./verify_ctrsm            # exit 0 = all pass, prints "N cases, F fails"
+ * SPDX-License-Identifier: MIT
+ */
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <complex.h>
+#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;j<n;j++) for(int i=0;i<n;i++){
+        cd v;
+        if(i==j) v = unit ? (1.0+0.0*I) : ((2.0+(double)(i%3)) + 0.5*I);
+        else if((uplo_upper && i<j) || (!uplo_upper && i>j)) 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 && sr<sc)||(uplo==CblasLower && sr>sc)) 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;i<m;i++) for(int jj=0;jj<n;jj++){
+            cd s=0.0;
+            for(int k=0;k<m;k++) s += opA(trans,uplo,is_unit,A,lda,i,k) * X[k+jj*ldb];
+            double d=cabs(s - B0[i+jj*ldb]); if(d>maxr) maxr=d;
+        }
+    } else {
+        for(int i=0;i<m;i++) for(int jj=0;jj<n;jj++){
+            cd s=0.0;
+            for(int k=0;k<n;k++) s += X[i+k*ldb]*opA(trans,uplo,is_unit,A,lda,k,jj);
+            double d=cabs(s - B0[i+jj*ldb]); if(d>maxr) 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;im<ng;im++) for(int in=0;in<ng;in++){
+        int m=grid[im], n=grid[in];
+        int side=sides[si],uplo=uplos[ui],trans=transs[ti],unit=diags[di];
+        int ao=(side==CblasLeft)?m:n; int lda=ao, ldb=m;
+        cd *A=malloc(sizeof(cd)*ao*ao);
+        cd *B=malloc(sizeof(cd)*m*n), *B0=malloc(sizeof(cd)*m*n);
+        make_tri(A,ao, uplo==CblasUpper, unit==CblasUnit);
+        for(int x=0;x<m*n;x++){ B[x]=crand(); B0[x]=B[x]; }
+        cblas_ctrsm(CblasColMajor,side,uplo,trans,unit,m,n,&alpha,A,lda,B,ldb);
+        double r=resid(side,uplo,trans,(unit==CblasUnit),m,n,A,lda,B,B0,ldb);
+        double tol=1e-4*(ao+1);
+        total++; if(r>worst) 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;
+}

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**,

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).