opensolvers · hugomeiland · Jul 19, 2026 · Jul 19, 2026
diff --git a/BLIS/README.md b/BLIS/README.md
@@ -0,0 +1,147 @@
+# BLIS on RISC-V - build, validation & performance vs OpenBLAS
+
+Builds [BLIS](https://github.com/flame/blis) (BLAS-like Library Instantiation
+Software, from the FLAME group) for the **SpaceMiT X60 / K1** (Orange Pi RV2 /
+BPI-F3, `rv64gcv`, RVV 1.0, VLEN=256) under the [EESSI](https://www.eessi.io/)
+GCC 14.3.0 toolchain, verifies it is numerically correct, and A/Bs its BLAS
+throughput against OpenBLAS.
+
+BLIS ships **hand-written RVV assembly** level-3 microkernels under
+[`kernels/rviv/3/`](https://github.com/flame/blis/tree/master/kernels/rviv/3)
+(dynamic-VLEN via `get_vlenb()`), selected by the `rv64iv` config target - so
+this is a real vector-vs-vector comparison, not vector-vs-scalar.
+
+| File | What it does |
+|---|---|
+| `build-blis.sh` | Clones + configures (`rv64iv`, `--enable-cblas --enable-blas`) + builds BLIS with EESSI GCC 14.3.0. Emits `libblis.a` + CBLAS/BLAS symbols. |
+| `run-ab.sh` | Links the backend-agnostic `bench_dgemm.c` against BLIS **and** OpenBLAS, sweeps sizes × threads, prints GFLOP/s side-by-side + a ratio. |
+| `bench_dgemm.c` | Square `C=A*B` (level-3 `dgemm`, 3 reps), reports GFLOP/s, prints `C[0]` so a NaN/garbage backend is visible. *(shared with [`../OpenBLAS/`](../OpenBLAS))* |
+| `difftest.c` | `dlopen`s a BLAS `.so` (or link a `.a`) and runs level-1/2/3 routines on fixed inputs, printing `sum / sumsq / nan / inf` per routine for differential comparison against a reference. |
+| `verify_ctrsm.c` | Full-parameter TRSM correctness sweep (residual-checked); confirms BLIS's triangular solve matches the reference. |
+
+## Why link, not FlexiBLAS
+
+The other BLAS-axis dirs ([`../OpenBLAS`](../OpenBLAS), [`../hpl`](../hpl),
+[`../numpy`](../numpy)) A/B by swapping the backend at runtime via FlexiBLAS.
+**FlexiBLAS is not installed on the RV2**, so here the A/B is done by *linking*
+the same unchanged `bench_dgemm.c` against each library in turn. This still
+honours the repo principle - change exactly one variable (the BLAS
+implementation), hold source and compiler flags constant.
+
+## Build
+
+```bash
+# on the Orange Pi RV2 (native build); ~10-20 min for a full BLIS build
+./build-blis.sh                       # installs to $HOME/blis-install
+# or:  ./build-blis.sh /custom/prefix
+```
+
+`build-blis.sh` forces EESSI **GCC 14.3.0** ahead of the EESSI compat GCC 13.4.0
+(see the gotcha below) and configures the RVV 1.0 target with OpenMP threading:
+
+```bash
+./configure --prefix=$PREFIX --enable-cblas --enable-blas \
+            --enable-threading=openmp rv64iv
+make -j$(nproc) && make install
+```
+
+> `rv64iv` = RVV 1.0 with **dynamic VLEN detection** (correct for the X60's
+> VLEN=256). Do **not** use `sifive_rvv`, which defaults to VLEN=128.
+> `--enable-threading=openmp` is required for BLIS to use all 8 cores; without
+> it DGEMM stays single-threaded (~2.7 GFLOP/s regardless of thread count).
+
+## Run
+
+```bash
+# performance A/B (BLIS vs OpenBLAS)
+BLIS_PREFIX=$HOME/blis-install \
+OPENBLAS_LIB=$HOME/trsm-pr5830/libopenblas.a \
+  ./run-ab.sh
+# tune the sweep:  SIZES="1024 2048 4096" THREADS="1 8" ./run-ab.sh
+
+# correctness of the BLIS build (threaded lib needs -fopenmp + the GCC14 flags)
+GCC14=/cvmfs/dev.eessi.io/riscv/versions/2025.06-001/software/linux/riscv64/generic/software/GCCcore/14.3.0
+gcc -O3 -march=rv64imafdcv_zvl256b -fopenmp verify_ctrsm.c -o verify_ctrsm \
+    -I$HOME/blis-install/include/blis -L$GCC14/lib -B$GCC14/lib \
+    $HOME/blis-install/lib/libblis.a -lm -lpthread
+OMP_NUM_THREADS=8 ./verify_ctrsm    # exit 0 = all pass
+```
+
+## A/B methodology
+
+- **One variable**: identical `bench_dgemm.c`, identical `-O3
+  -march=rv64imafdcv_zvl256b`, only the linked BLAS differs.
+- **OpenBLAS baseline**: the patched RVV OpenBLAS from this repo's work
+  (`$HOME/trsm-pr5830/libopenblas.a`, `ZVL128B`/`ZVL256B` RVV kernels).
+- **Threads**: 1 (single-core kernel efficiency) and 8 (full X60 parallel).
+- **Correctness first**: `verify_ctrsm` / `difftest` must pass before any speed
+  number is trusted - a fast NaN is worthless.
+
+## Results
+
+Measured on the Orange Pi RV2 (SpaceMiT X60, 8 cores, VLEN=256), BLIS
+`061c2eb` (`rv64iv`, OpenMP threading) vs the patched RVV OpenBLAS
+`0.3.33.dev` (`zvl128bp`), both `-O3 -march=rv64imafdcv_zvl256b`, EESSI
+GCC 14.3.0. DGEMM GFLOP/s (square, 3 reps, best):
+
+| threads | N | BLIS | OpenBLAS | BLIS/OpenBLAS |
+|--:|--:|--:|--:|--:|
+| 1 | 1024 | 1.99 | 2.13 | 0.93x |
+| 1 | 2048 | 2.73 | 2.25 | 1.21x |
+| 1 | 4096 | 2.95 | 2.28 | 1.29x |
+| 8 | 1024 | 8.93 | 10.13 | 0.88x |
+| 8 | 2048 | 9.60 | 10.83 | 0.89x |
+| 8 | 4096 | 9.55 | 11.94 | 0.80x |
+
+Takeaways:
+
+- **Single thread, large N**: BLIS's hand-written RVV assembly gemm microkernel
+  (`kernels/rviv/3/`) **beats OpenBLAS by ~20-30%** once the problem is big
+  enough to amortize packing (N>=2048). At N=1024 the two are within noise.
+- **8 threads**: OpenBLAS scales slightly better (BLIS reaches 0.80-0.89x). Both
+  get ~3.5-5x from 8 cores; BLIS's OpenMP path leaves some headroom vs
+  OpenBLAS's threading.
+- `C[0]=245.24` was identical for both backends at every size - no NaN/garbage.
+
+Correctness (`verify_ctrsm`, BLIS `rv64iv`): **2400 cases, 0 fails,
+worst_resid=2.55e-07** (single- and 8-thread). BLIS's triangular solve is
+numerically correct on the X60 - unlike the stock OpenBLAS `_rvv_v1` TRSM kernel
+this repo's [`../OpenBLAS`](../OpenBLAS) work localized.
+
+## Gotcha - `module load` does not repath `gcc` on the Orange Pi RV2
+
+Same trap as [`../fftw/README.md`](../fftw/README.md): `module load
+GCCcore/14.3.0` returns rc=0 but does **not** put GCC 14 first on `PATH` - the
+EESSI compat-layer GCC 13.4.0 keeps winning. Both scripts here force it
+explicitly, *before* any `set -e` (the lmod `module` function reads unbound vars
+and would kill a strict shell):
+
+```bash
+GCC14=/cvmfs/dev.eessi.io/riscv/versions/2025.06-001/software/linux/riscv64/generic/software/GCCcore/14.3.0
+export PATH="$GCC14/bin:$PATH"
+export LD_LIBRARY_PATH="$GCC14/lib64:$GCC14/lib:$LD_LIBRARY_PATH"
+```
+
+## Gotcha - GCC14 runtime libs are off the default linker search path
+
+The shared `libblis.so` link (and any `-fopenmp` link against the static lib)
+fails with the EESSI GCC 14.3.0 unless you point the driver at `$GCC14/lib`,
+where its `libgcc_s`, `libatomic`, and `libgomp.spec` actually live:
+
+```
+/usr/bin/ld: cannot find -lgcc_s / -latomic          # → needs -L$GCC14/lib
+gcc: fatal error: cannot read spec file 'libgomp.spec'  # → needs -B$GCC14/lib
+```
+
+Both scripts pass `-L$GCC14/lib -B$GCC14/lib -Wl,-rpath,$GCC14/lib` on every
+link. `-L` resolves the libs, `-B` lets the gcc driver read `libgomp.spec` for
+`-fopenmp`, and the rpath lets the resulting binary find `libgomp.so.1` at run
+time. (The static-linked `bench_dgemm` still needs the rpath because BLIS's
+OpenMP calls pull in `libgomp` dynamically.)
+
+## References
+
+- BLIS RISC-V config: [`config_registry`](https://github.com/flame/blis/blob/master/config_registry) (`rv64iv`, `sifive_rvv`)
+- RVV assembly gemm kernels: [`kernels/rviv/3/`](https://github.com/flame/blis/tree/master/kernels/rviv/3)
+- BLIS RISC-V CI: [`ci/do_riscv.sh`](https://github.com/flame/blis/blob/master/ci/do_riscv.sh)
+- Prior RISC-V work: PR [#737](https://github.com/flame/blis/pull/737) (SiFive X280), [#832](https://github.com/flame/blis/pull/832) (generic RVV, VLEN-tunable), [#868](https://github.com/flame/blis/pull/868) (SG2042, +49% sdgemm)
diff --git a/BLIS/bench_dgemm.c b/BLIS/bench_dgemm.c
@@ -0,0 +1,36 @@
+/*
+ * bench_dgemm.c - minimal DGEMM throughput benchmark.
+ *
+ * Times a square double-precision C = A*B (level-3 BLAS dgemm, 3 reps) and
+ * reports GFLOP/s. Links against any BLAS that provides the Fortran dgemm_
+ * symbol; use FlexiBLAS (or OPENBLAS_CORETYPE / OPENBLAS_NUM_THREADS) to select
+ * and A/B the backend without recompiling. Prints C[0] so a NaN/garbage backend
+ * is immediately visible.
+ *
+ * Build:  gcc -O2 bench_dgemm.c -o bench_dgemm -lflexiblas   (or -lopenblas)
+ * Run:    OPENBLAS_NUM_THREADS=8 ./bench_dgemm [N]           # default N=4096
+ * A/B:    scalar : OPENBLAS_CORETYPE=RISCV64_GENERIC ./bench_dgemm
+ *         vector : FLEXIBLAS=/path/to/libopenblas.so  ./bench_dgemm
+ * SPDX-License-Identifier: MIT
+ */
+#include <stdio.h>
+#include <stdlib.h>
+#include <time.h>
+extern void dgemm_(const char*,const char*,const int*,const int*,const int*,
+  const double*,const double*,const int*,const double*,const int*,
+  const double*,double*,const int*);
+static double now(){struct timespec t;clock_gettime(CLOCK_MONOTONIC,&t);return t.tv_sec+t.tv_nsec*1e-9;}
+int main(int argc,char**argv){
+  int n=argc>1?atoi(argv[1]):4096;
+  double *A=malloc(8L*n*n),*B=malloc(8L*n*n),*C=malloc(8L*n*n);
+  if(!A||!B||!C){printf("alloc fail\n");return 1;}
+  for(long i=0;i<(long)n*n;i++){A[i]=(i%7)*0.1;B[i]=(i%5)*0.2;C[i]=0;}
+  double al=1.0,be=0.0;
+  dgemm_("N","N",&n,&n,&n,&al,A,&n,B,&n,&be,C,&n);
+  int reps=3; double t0=now();
+  for(int r=0;r<reps;r++) dgemm_("N","N",&n,&n,&n,&al,A,&n,B,&n,&be,C,&n);
+  double t=now()-t0; double gf=2.0*(double)n*n*n*reps/t/1e9;
+  printf("N=%d threads=%s DGEMM=%.2f GFLOP/s (t=%.2fs C[0]=%g)\n",
+    n,getenv("OPENBLAS_NUM_THREADS")?getenv("OPENBLAS_NUM_THREADS"):"?",gf,t,C[0]);
+  return 0;
+}
diff --git a/BLIS/build-blis.sh b/BLIS/build-blis.sh
@@ -0,0 +1,63 @@
+#!/bin/bash
+# build-blis.sh - build BLIS (flame/blis) for the SpaceMiT X60 / K1 (rv64gcv,
+# RVV 1.0, VLEN=256) under the EESSI GCC 14.3.0 toolchain, with the CBLAS/BLAS
+# compatibility layer so it links against the repo's backend-agnostic harnesses.
+#
+# BLIS target `rv64iv` = 64-bit RVV 1.0 with dynamic VLEN detection (correct for
+# the X60's VLEN=256); it selects the hand-written RVV assembly gemm microkernels
+# under kernels/rviv/3/. (The `sifive_rvv` config defaults to VLEN=128 and is the
+# wrong choice here.)
+#
+# Usage:  ./build-blis.sh [install_prefix]     # default prefix: $HOME/blis-install
+# Output: $PREFIX/lib/libblis.a  (+ CBLAS/BLAS symbols) , $PREFIX/include/blis/
+#
+# NOTE: run ON the Orange Pi RV2 (native build). Modules are loaded BEFORE
+# `set -e` because on the RV2 the lmod `module` function returns nonzero / reads
+# unbound vars and would kill a strict shell (see ../README.md gotcha).
+
+# --- EESSI GCC 14.3.0: module load does NOT repath gcc on the RV2, so force it ---
+# (identical workaround to fftw/README.md - EESSI compat GCC 13.4.0 otherwise wins)
+GCC14=/cvmfs/dev.eessi.io/riscv/versions/2025.06-001/software/linux/riscv64/generic/software/GCCcore/14.3.0
+export PATH="$GCC14/bin:$PATH"
+export LD_LIBRARY_PATH="$GCC14/lib64:$GCC14/lib:$LD_LIBRARY_PATH"
+# GCC14's libgcc_s / libatomic live in $GCC14/lib, off the default linker search
+# path; without this the shared libblis.so link fails on -lgcc_s / -latomic.
+export LDFLAGS="-L$GCC14/lib -B$GCC14/lib -Wl,-rpath,$GCC14/lib ${LDFLAGS:-}"
+
+set -euo pipefail
+
+PREFIX=${1:-$HOME/blis-install}
+SRC=${BLIS_SRC:-$HOME/blis}
+JOBS=${JOBS:-$(nproc)}
+
+echo "=== toolchain ==="
+gcc --version | head -1
+gcc -dumpmachine
+echo "prefix : $PREFIX"
+echo "src    : $SRC"
+echo "jobs   : $JOBS"
+
+# --- clone if needed ---
+if [ ! -d "$SRC/.git" ]; then
+  echo "=== cloning flame/blis into $SRC ==="
+  git clone --depth 1 https://github.com/flame/blis.git "$SRC"
+fi
+cd "$SRC"
+echo "=== BLIS commit ==="; git rev-parse --short HEAD
+
+# --- configure for RVV 1.0, dynamic VLEN, with CBLAS + BLAS compat ---
+echo "=== configure rv64iv ==="
+./configure --prefix="$PREFIX" --enable-cblas --enable-blas --enable-threading=openmp rv64iv
+
+echo "=== make -j$JOBS ==="
+make -j"$JOBS"
+
+echo "=== make install ==="
+make install
+
+echo "=== result ==="
+ls -la "$PREFIX"/lib/libblis.* 2>/dev/null || true
+echo "cblas symbols present?"
+nm -g "$PREFIX"/lib/libblis.a 2>/dev/null | grep -c -E "cblas_dgemm|cblas_dtrsm" || true
+echo
+echo "BLIS built. Next: ./run-ab.sh (bench vs OpenBLAS) and ./verify (correctness)."
diff --git a/BLIS/difftest.c b/BLIS/difftest.c
@@ -0,0 +1,111 @@
+/*
+ * difftest.c - differential BLAS correctness test (single-threaded).
+ *
+ * dlopen()s a given BLAS .so and calls a spread of level-1/2/3 routines on
+ * fixed pseudo-random inputs, printing sum / sumsq / NaN-count / Inf-count for
+ * each. Point it at two backends (or one backend under two OPENBLAS_CORETYPE
+ * settings) and diff the output to localize a numerically broken kernel.
+ *
+ * This is the test that isolated the OpenBLAS 0.3.30 RISC-V vector gemv_n bug:
+ * on a SpaceMiT X60 the stock DYNAMIC_ARCH build's RVV (ZVL256B) cblas_dgemv
+ * returns NaN (VFILL_ZERO zeroes an uninitialized vector register) while
+ * dgemm/dtrsm are correct - which is exactly why HPL (leans on dgemv in panel
+ * factorization) NaNs while a standalone dgemm benchmark looks fine. Forcing
+ * OPENBLAS_CORETYPE=RISCV64_GENERIC, or using a fixed/patched RVV build, makes
+ * dgemv match the scalar reference bit-for-bit (nan=0).
+ *
+ * Build:  gcc -O2 difftest.c -o difftest -ldl -lm
+ * Run:    ./difftest /path/to/libopenblas.so
+ *         OPENBLAS_CORETYPE=RISCV64_GENERIC ./difftest /path/to/libopenblas.so
+ * Requires a backend exposing the cblas_* symbols (OpenBLAS does).
+ * SPDX-License-Identifier: MIT
+ */
+#define _GNU_SOURCE
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <dlfcn.h>
+
+enum {ColMajor=102,RowMajor=101,NoTrans=111,Trans=112,Upper=121,Lower=122,
+      Left=141,Right=142,NonUnit=131,Unit=132};
+
+typedef void   (*dgemm_t)(int,int,int,int,int,int,double,const double*,int,const double*,int,double,double*,int);
+typedef void   (*dtrsm_t)(int,int,int,int,int,int,int,double,const double*,int,double*,int);
+typedef void   (*dgemv_t)(int,int,int,int,double,const double*,int,const double*,int,double,double*,int);
+typedef void   (*dscal_t)(int,double,double*,int);
+typedef void   (*daxpy_t)(int,double,const double*,int,double*,int);
+typedef size_t (*idamax_t)(int,const double*,int);
+typedef double (*dnrm2_t)(int,const double*,int);
+typedef double (*ddot_t)(int,const double*,int,const double*,int);
+typedef void   (*setnt_t)(int);
+
+static double rnd(long i){
+  unsigned long x=(unsigned long)i*2654435761UL+1013904223UL;
+  x^=x>>13; x*=0xD1B54A32D192ED03UL; x^=x>>29;
+  return ((double)(x>>40)/(double)(1UL<<24))*2.0-1.0; /* [-1,1) */
+}
+
+static void report(const char*name,const double*v,long n){
+  double s=0,s2=0; long nan=0,inf=0;
+  for(long i=0;i<n;i++){double x=v[i];
+    if(isnan(x))nan++; else if(isinf(x))inf++; else {s+=x; s2+=x*x;}}
+  printf("%-8s sum=%.17g sumsq=%.17g nan=%ld inf=%ld\n",name,s,s2,nan,inf);
+}
+
+int main(int argc,char**argv){
+  if(argc<2){fprintf(stderr,"usage: %s lib.so\n",argv[0]);return 2;}
+  void*h=dlopen(argv[1],RTLD_NOW|RTLD_LOCAL);
+  if(!h){fprintf(stderr,"dlopen: %s\n",dlerror());return 3;}
+  #define SYM(t,n) t n=(t)dlsym(h,#n); if(!n){fprintf(stderr,"missing %s\n",#n);return 4;}
+  SYM(setnt_t,openblas_set_num_threads); openblas_set_num_threads(1);
+  SYM(dgemm_t,cblas_dgemm)
+  SYM(dtrsm_t,cblas_dtrsm)
+  SYM(dgemv_t,cblas_dgemv)
+  SYM(dscal_t,cblas_dscal)
+  SYM(daxpy_t,cblas_daxpy)
+  SYM(idamax_t,cblas_idamax)
+  SYM(dnrm2_t,cblas_dnrm2)
+  SYM(ddot_t,cblas_ddot)
+
+  /* ---- level-1 vector ops ---- */
+  long nv=4096;
+  double*x=malloc(nv*sizeof(double)),*y=malloc(nv*sizeof(double));
+  for(long i=0;i<nv;i++){x[i]=rnd(i); y[i]=rnd(i+100000);}
+  size_t im=cblas_idamax(nv,x,1);
+  printf("idamax   idx=%zu val=%.17g\n",im,(im<(size_t)nv?x[im]:0.0/0.0));
+  printf("dnrm2    val=%.17g\n",cblas_dnrm2(nv,x,1));
+  printf("ddot     val=%.17g\n",cblas_ddot(nv,x,1,y,1));
+  double*xs=malloc(nv*sizeof(double)); memcpy(xs,x,nv*sizeof(double));
+  cblas_dscal(nv,3.14159,xs,1); report("dscal",xs,nv);
+  double*ya=malloc(nv*sizeof(double)); memcpy(ya,y,nv*sizeof(double));
+  cblas_daxpy(nv,-2.5,x,1,ya,1); report("daxpy",ya,nv);
+
+  /* ---- dgemv (col-major, NoTrans) ---- */
+  int M=768,N=768;
+  double*A=malloc((size_t)M*N*sizeof(double));
+  for(long i=0;i<(long)M*N;i++)A[i]=rnd(i+1);
+  double*xv=malloc(N*sizeof(double)),*yv=malloc(M*sizeof(double));
+  for(int i=0;i<N;i++)xv[i]=rnd(i+7);
+  for(int i=0;i<M;i++)yv[i]=rnd(i+9);
+  cblas_dgemv(ColMajor,NoTrans,M,N,1.3,A,M,xv,1,0.7,yv,1); report("dgemv",yv,M);
+
+  /* ---- dgemm sanity ---- */
+  int mm=256;
+  double*Ag=malloc((size_t)mm*mm*sizeof(double)),*Bg=malloc((size_t)mm*mm*sizeof(double)),*Cg=malloc((size_t)mm*mm*sizeof(double));
+  for(long i=0;i<(long)mm*mm;i++){Ag[i]=rnd(i+3);Bg[i]=rnd(i+5);Cg[i]=0;}
+  cblas_dgemm(ColMajor,NoTrans,NoTrans,mm,mm,mm,1.0,Ag,mm,Bg,mm,0.0,Cg,mm); report("dgemm",Cg,(long)mm*mm);
+
+  /* ---- dtrsm (Left,Lower,NoTrans,NonUnit) diag-dominant ---- */
+  int tm=512,tn=64;
+  double*T=malloc((size_t)tm*tm*sizeof(double));
+  for(int j=0;j<tm;j++)for(int i=0;i<tm;i++){
+    double v=(i==j)?(tm+1.0):(i>j?rnd((long)j*tm+i)*0.1:0.0);
+    T[(size_t)j*tm+i]=v; }
+  double*B=malloc((size_t)tm*tn*sizeof(double));
+  for(long i=0;i<(long)tm*tn;i++)B[i]=rnd(i+13);
+  cblas_dtrsm(ColMajor,Left,Lower,NoTrans,NonUnit,tm,tn,1.0,T,tm,B,tm); report("dtrsm",B,(long)tm*tn);
+
+  printf("DONE\n");
+  return 0;
+}