RISC-V (riscv64_zvl256b, SpaceMiT K1): DGEMM regression between 0.3.31 and 0.3.33 produces non-PSD matrix from A^T @ A

On SpaceMiT K1 hardware (RVA22 + RVV 1.0 (VLEN=256), runtime DYNAMIC_ARCH selects riscv64_zvl256b), DGEMM produces a matrix with significantly negative eigenvalues when computing A.T @ A for a 50x50 random matrix. Mathematically the result must be positive semi-definite. This regressed between OpenBLAS 0.3.31 and 0.3.33.

The downstream symptom is numpy.linalg.cholesky raising LinAlgError: Matrix is not positive definite on a matrix that should have been PSD by construction.

Hardware

BananaPi F3 (SpaceMiT K1 SoC, 8x X60 cores, RVA22 + RVV 1.0, VLEN=256)
Armbian 25.11.2 trixie, Linux 6.6.99-current-spacemit
Also confirmed on Muse Pi Pro (also K1) by another user
NOT reproducible on Pine64 Star64 (SiFive U74, no RVV, dispatches to riscv64_generic), which sharpens the suspicion that the regression is in the riscv64_zvl256b kernels

Build configuration

OpenBLAS built with DYNAMIC_ARCH, USE64BITINT, NO_AFFINITY. Runtime dispatch picks riscv64_zvl256b on the K1 (verified with OPENBLAS_VERBOSE=2).

Bisect

Tested via LD_PRELOAD swap of the OpenBLAS shared library against an unchanged numpy 2.5.0.dev0 wheel. All 0.3.31.x scipy-openblas64 builds produce correct results; 0.3.33.0.0 produces wrong results.

scipy-openblas64	OpenBLAS embedded	min eigval of A^T @ A	cholesky
0.3.31.22.1	0.3.31.dev	+1.405e-02	OK
0.3.31.126.1	0.3.31.dev	+1.405e-02	OK
0.3.31.126.5	0.3.31.dev	+1.405e-02	OK
0.3.31.159.0	0.3.31.159.0	+1.405e-02	OK
0.3.31.188.0	0.3.31.188.0	+1.405e-02	OK
0.3.33.0.0	0.3.33.0.0	-2.779e+00	LinAlgError

The Debian 13 trixie system OpenBLAS 0.3.29 (libopenblas0:riscv64 0.3.29+ds-3) also produces correct results on the same machine.

scipy-openblas64 has no 0.3.32 release on PyPI, so the regression range at OpenBLAS source level is v0.3.31..v0.3.33 (~250 commits).

Minimal reproducer

Requires Python with numpy linked against the OpenBLAS in question:

import numpy as np

np.random.seed(1)
a = np.random.randn(50, 50)
a = (a.T @ a).astype(np.float64)  # mathematically PSD

print("min eigval =", np.linalg.eigvalsh(a).min())
# OpenBLAS 0.3.31: +1.405e-02 (correct)
# OpenBLAS 0.3.33: -2.779e+00 (wrong by ~2 orders of magnitude)

np.linalg.cholesky(a, upper=True)  # raises LinAlgError on 0.3.33

The matrix is dense, symmetric (within FP roundoff), 50x50, float64. The error magnitude (~2.8) is far beyond floating-point precision. DGEMM is producing an actively wrong result, not a tiny precision drift.

Context

Surfaced in numpy/numpy#30995 (riscv64 wheel CI). Confirmed by multiple K1 users. Cross-referencing here so OpenBLAS maintainers can pick up the bisect.

Available to test patches on the F3 if useful.

@threexc

Correction to the original report: the K1 X60 implements RVV 1.0 (VLEN=256), not generic riscv64. Verified with OPENBLAS_VERBOSE=2 against both 0.3.31.188.0 and 0.3.33.0.0 wheels, both dispatch to riscv64_zvl256b on the K1.

The bisect range is unchanged; the regression is in the riscv64_zvl256b path (or something it calls), not the scalar generic kernels. Title and body edited accordingly.

Thanks to @threexc on numpy/numpy#30995 for catching this.

More narrowing from runtime checks on the F3, all using the same broken 0.3.33 wheel.

Confirmed: bug is in riscv64_zvl256b, not generic.

OPENBLAS_CORETYPE=riscv64_generic, n=50:    min eigval = +1.405e-02  OK
OPENBLAS_CORETYPE=(default zvl256b), n=50:  min eigval = -2.779e+00  FAIL

Same .so, same hardware, just changing the dispatched kernel. Generic produces correct results; zvl256b does not.

Not concurrency. OPENBLAS_NUM_THREADS=1 with zvl256b still fails. Deterministic kernel bug.

Both precisions broken. At n=50, seed=1: SGEMM min eigval = -2.599e+00, DGEMM min eigval = -2.779e+00.

Size pattern: failures when n % 8 ∈ {1, 2, 3}.

Sweep of n=33..65 with seed=1 on zvl256b:

33-35 FAIL  (n%8 = 1,2,3)
36-40 OK    (n%8 = 4,5,6,7,0)
41-43 FAIL
44-48 OK
49-51 FAIL
52-56 OK
57-59 FAIL
60-64 OK
65    FAIL

Perfectly periodic, period 8. With VLEN=256 and float64, that is 4 doubles per vector register, 8 doubles across two registers. The breakage happens when the trailing M or N past the last full block is 1, 2, or 3 elements, which looks like broken tail/cleanup for output dims smaller than half a register's worth.

The earlier coarser sweep (4, 8, 16, 20, 32, 64, 100, 128, 200, 256, 500) all happened to be n%8 ∈ {0, 4}, which is why the bug seemed to only hit n=50.

Seed-independent. 20/20 random seeds at n=50 produce FAIL.

K dimension gates it too. Holding M=N=50, varying K:

A shape	K	K%8	result
40x50	40	0	FAIL
50x50	50	2	FAIL
80x50	80	0	OK
100x50	100	4	OK

K alone doesn't explain it (40 and 80 are both 0 mod 8 with different outcomes); likely the GEMM blocking driver picks a different code path past some K threshold.

So: regression is in the riscv64_zvl256b GEMM kernel (or its blocking driver), specifically the tail handling for M%8 ∈ {1,2,3}, present in both SGEMM and DGEMM, with additional K-dependent gating.

Happy to run more targeted tests on the F3 if useful.

Right, the RISCV64_GENERIC target would have been very unlikely to be broken (simple 2x2 GEMM kernel shared with all other generic and a few not so generic targets).
All relevant recent work on the ZVL256B DGEMM (and SGEMM) kernel occured in the 0.3.33 timeframe, so 0.3.32 should be unaffected and there are only a handful of PRs between late March and early April that could have had this side effect IMHO.
Unfortunately I am away from my BPI-F3 and the one in the GCC Compile Farm appears to be down. I'll check again later, IIRC that machine already has Numpy installed

@ChipKerchner

Source bisect on a BPI-F3 (K1) points to one commit and one BLAS routine.

First-bad commit: d69be17b6ff7 (first commit of PR #5674, "Convert 2X LMUL1 instructions to 1X LMUL2. Improved FP64 GEMM edges").

The bug shows up in cblas_dsyrk with Upper + Trans. Plain cblas_dgemm at the same size returns correct results. Error magnitude (1.535e+01) matches the original numpy.matmul reproducer.

Caveat before the data: I got the K1's ISA wrong earlier in this thread, so please double-check my work. I am not a BLAS expert.

Method

Built OpenBLAS at each candidate commit on the F3 with TARGET=RISCV64_ZVL256B INTERFACE64=1 SYMBOLPREFIX=scipy_ SYMBOLSUFFIX=64_ NO_LAPACK=1 ONLY_CBLAS=1, then called the resulting scipy_cblas_*64_ symbols directly via Python ctypes against a pure-Python triple-loop reference. No LAPACK, no numpy in the call path, just the BLAS kernel.

Per-commit results

sha          label                                                max_err   verdict
4956446c     baseline pre-LMUL2                                   2.13e-14   GOOD
d69be17b     Convert 2X LMUL1 -> 1X LMUL2 FP64 GEMM edges        1.54e+01   BAD  <<<
daa3215f     Remove shadow variable                               1.54e+01   BAD
3b1aef12     Use LMUL2 loads main block                           1.54e+01   BAD
22b7950b     Use LMUL2 calculations main block                    1.54e+01   BAD
cc1b5794     Reduce vectors 32 -> 24                              1.54e+01   BAD
v0.3.33      release tag                                          1.54e+01   BAD

Why dgemm tests look fine but dsyrk doesn't

At M=N=K=50, cblas_dgemm returns max_err = 2.13e-14 on every commit including v0.3.33. I spent a while confused by this before checking what numpy actually calls.

Numpy's matmul(X.T, X) dispatches to cblas_dsyrk(RowMajor, Upper, Trans, N=50, K=50). Looks like it detects A.T @ A and uses the symmetric rank-K update. On the BAD wheel, dsyrk Upper+Trans returns max_err = 1.535e+01, matching the matmul error exactly. Lower+Trans is correct. Plain dgemm is correct.

Speculation: dsyrk's diagonal blocks call dgemm with smaller M than the main block size, and that hits the new LMUL2 edge-handling path. The full 50x50x50 GEMM stays in the main block code, which is fine.

Reproducer (small)

import ctypes, numpy as np
lib = ctypes.CDLL("/path/to/libscipy_openblas64_.so")
fn = lib.scipy_cblas_dsyrk64_
fn.argtypes = [ctypes.c_int]*3 + [ctypes.c_int64]*2 + [ctypes.c_double] + \
    [ctypes.POINTER(ctypes.c_double), ctypes.c_int64, ctypes.c_double,
     ctypes.POINTER(ctypes.c_double), ctypes.c_int64]
fn.restype = None

np.random.seed(1)
X = np.random.randn(50,50).astype(np.float64, order='C')
A = np.zeros((50,50), order='C')
fn(101, 121, 112, 50, 50, 1.0,
   X.ctypes.data_as(ctypes.POINTER(ctypes.c_double)), 50,
   0.0,
   A.ctypes.data_as(ctypes.POINTER(ctypes.c_double)), 50)

ref = X.T @ X  # any GOOD BLAS will give this
print("max_err:", np.abs(np.triu(A) - np.triu(ref)).max())
# K1 with v0.3.33 zvl256b -> ~1.5e+01
# K1 forced to riscv64_generic -> 2.1e-14

Hardware

BPI-F3 (SpaceMiT K1, 8x X60, RVA22 + RVV 1.0, VLEN=256). Default DYNAMIC_ARCH picks riscv64_zvl256b at runtime. OPENBLAS_CORETYPE=riscv64_generic produces correct results.

If a patch would be useful to try on the F3, I can run it. Build artifacts and bisect logs are still on the machine.

cc @ChipKerchner

Thanks, I'll see if I can reproduce it and fix it

@gounthar

@gounthar Could you list the tools you are using - AKA compilers including version?

@ChipKerchner

@ChipKerchner here's the toolchain on the BananaPi F3 (SpaceMiT K1) where I ran the bisect:

OS: Armbian 25.11.2 (Debian 13 trixie)
Kernel: 6.6.99-current-spacemit (SMP PREEMPT_DYNAMIC, 2025-04-10)
GCC: 14.2.0 (Debian 14.2.0-19)
GNU Make 4.4.1
CPU: Spacemit X60, 8 cores at 1.6 GHz max, rv64gc + RVV 1.0, VLEN=256

Each commit was built with TARGET=RISCV64_ZVL256B INTERFACE64=1 NO_LAPACK=1 NO_LAPACKE=1 ONLY_CBLAS=1 SYMBOLPREFIX=scipy_ SYMBOLSUFFIX=64_ to match the scipy-openblas64 wheel layout. I can rebuild with different flags if that would help narrow it down further.

I can reproduce a similar behavior that you are seeing with scipy.

I have retested the dense GEMM and see no issues (other than some dead code that needs to be removed - probably an issue with merging).

I suspect that it may be that ldc doesn't match M that is causing the differences but I need to dig into it more. I'll keep you updated.

I'm not familiar with building libscipy_openblas64_.so. Could you provide instructions of how I build that library including using my own version of libopenblas?

The canonical recipe lives in MacPython/openblas-libs (the repo that produces the scipy-openblas64 wheels): https://github.com/MacPython/openblas-libs/blob/main/tools/build_steps.sh

For a quick test build, no LAPACK, CBLAS only (what I used for the bisect):

make -j$(nproc) \
  TARGET=RISCV64_ZVL256B \
  INTERFACE64=1 \
  SYMBOLPREFIX=scipy_ SYMBOLSUFFIX=64_ \
  LIBNAMEPREFIX=scipy_ LIBNAMESUFFIX=64_ \
  FIXED_LIBNAME=1 \
  NO_LAPACK=1 NO_LAPACKE=1 ONLY_CBLAS=1 \
  BUFFERSIZE=20

make PREFIX=/some/prefix install

That produces libscipy_openblas64_.so with the scipy-mangled symbols (scipy_cblas_dgemm64_, scipy_cblas_dsyrk64_, etc). TARGET can be whatever your hardware needs; I used RISCV64_ZVL256B on the K1 (vlen 256). For a wheel-equivalent full build, drop NO_LAPACK NO_LAPACKE ONLY_CBLAS and use DYNAMIC_ARCH=1 instead of a fixed TARGET.

To swap your .so into an existing scipy-openblas64 wheel for testing, it goes in scipy_openblas64/lib/ as libscipy_openblas64_.so.0 (versioned soname); rename or symlink to match.

I can run the resulting .so against the cholesky reproducer on the K1 if that helps narrow the bisect further.

Could you also provide a version of your test for float32?

FYI - I have access to a K1 (hopefully soon as K3).

ssyrk version, same shape as the dsyrk reproducer (Upper/Trans, N=K=50). INTERFACE64=1 gives int64 indices; the symbol name comes from SYMBOLPREFIX=scipy_ plus SYMBOLSUFFIX=64_.

import ctypes, numpy as np

lib = ctypes.CDLL("/path/to/libscipy_openblas64_.so")
fn = lib.scipy_cblas_ssyrk64_
fn.restype = None
fn.argtypes = [
    ctypes.c_int, ctypes.c_int, ctypes.c_int,
    ctypes.c_int64, ctypes.c_int64,
    ctypes.c_float,
    ctypes.POINTER(ctypes.c_float), ctypes.c_int64,
    ctypes.c_float,
    ctypes.POINTER(ctypes.c_float), ctypes.c_int64,
]

CblasRowMajor, CblasUpper, CblasTrans = 101, 121, 112
N = K = 50
A = np.random.RandomState(0).randn(K, N).astype(np.float32)
C_blas = np.zeros((N, N), dtype=np.float32)
fn(CblasRowMajor, CblasUpper, CblasTrans, N, K, 1.0,
   A.ctypes.data_as(ctypes.POINTER(ctypes.c_float)), N,
   0.0, C_blas.ctypes.data_as(ctypes.POINTER(ctypes.c_float)), N)

C_ref = A.T @ A
mask = np.triu(np.ones((N, N), dtype=bool))
err = np.max(np.abs(C_blas[mask] - C_ref[mask]))
print(f"max abs error: {err:.3e}")

Direct K1 access will iterate faster than round-tripping through me. I can re-run on the BPI-F3 (K1, 8 cores, 16GB) if you want a second data point, or bisect commits between 0.3.29 and 0.3.33 there.

Your above make command does NOT produce a .so file. Any ideas?

My recipe was wrong; I pasted flags from the canonical MacPython wheel build instead of what I actually ran on the F3. The stripped-down version (no LIBNAMEPREFIX/LIBNAMESUFFIX/FIXED_LIBNAME) is what produced a usable .so:

make -j$(nproc) \
  TARGET=RISCV64_ZVL256B \
  INTERFACE64=1 \
  SYMBOLPREFIX=scipy_ SYMBOLSUFFIX=64_ \
  NO_LAPACK=1 NO_LAPACKE=1 ONLY_CBLAS=1

That dropped libopenblas_riscv64_zvl256b_-r0.3.33.dev.so (or a similarly version-suffixed name) into the build root. The scipy-mangled symbols (scipy_cblas_dsyrk64_ and friends) live inside that .so even though the filename itself doesn't carry the scipy_ prefix.

If that still doesn't reproduce, could you paste the tail of the make output plus find . -name "*.so*" from the build root? My hunch: FIXED_LIBNAME=1 combined with ONLY_CBLAS=1 interacted badly and skipped the shared lib stage entirely; the MacPython recipe pairs FIXED_LIBNAME with a full LAPACK build, not a CBLAS-only one. I can run any variant you want on the F3 today if it helps narrow it down.

@martin-frbg

Still doesn't produce the .so. @martin-frbg Any ideas why the above directions does not produce this file?

You shouldn't need the scipy-openblas .so I'd think. The only thing that matters is the regular OpenBLAS build settings (TARGET, INTERFACE64, etc.); no need to deal with symbol or library names, that's irrelevant to any numerical issues.

@ChipKerchner

Right, the symbol and library naming has no bearing on the numbers. My bisect used the scipy-named build only because I was matching the wheel layout; the wrong results come out the same from a plain TARGET=RISCV64_ZVL256B INTERFACE64=1 build, so that is all @ChipKerchner needs.

The bisect tree is still set up on the BPI-F3 (K1, zvl256b). Point me at a branch or a patch and I can build it there and run the dsyrk/cholesky reproducer, so the kernel work does not have to wait on K1 access. Whatever is least work for you.

Turns out I was in the wrong directory to produce the .so file.

Ok, found the problem. Issue was it was checking in M == ldc but it also needs to be less than a block width or it doesn't work correctly (in a small amount of cases)

SGEMM: max abs error: 2.861e-06
DGEMM: max_err: 7.105427357601002e-15

Navigation Menu

Repository navigation

RISC-V (riscv64_zvl256b, SpaceMiT K1): DGEMM regression between 0.3.31 and 0.3.33 produces non-PSD matrix from A^T @ A #5811

Hardware

Build configuration

Bisect

Minimal reproducer

Context

Method

Per-commit results

Why dgemm tests look fine but dsyrk doesn't

Reproducer (small)

Hardware

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

RISC-V (riscv64_zvl256b, SpaceMiT K1): DGEMM regression between 0.3.31 and 0.3.33 produces non-PSD matrix from A^T @ A #5811

Description

Hardware

Build configuration

Bisect

Minimal reproducer

Context

Activity

gounthar commented on May 14, 2026

gounthar commented on May 14, 2026

martin-frbg commented on May 14, 2026

gounthar commented on May 14, 2026

Method

Per-commit results

Why dgemm tests look fine but dsyrk doesn't

Reproducer (small)

Hardware

ChipKerchner commented on May 14, 2026

ChipKerchner commented on May 15, 2026

gounthar commented on May 15, 2026

ChipKerchner commented on May 15, 2026

ChipKerchner commented on May 15, 2026

gounthar commented on May 15, 2026

ChipKerchner commented on May 15, 2026

gounthar commented on May 15, 2026

ChipKerchner commented on May 15, 2026

gounthar commented on May 15, 2026

ChipKerchner commented on May 18, 2026

rgommers commented on May 18, 2026

gounthar commented on May 18, 2026

ChipKerchner commented on May 18, 2026

ChipKerchner commented on May 18, 2026

ChipKerchner commented on May 18, 2026

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions