Incorrect `dgemm_` results when called concurrently with `dsyevd_`

[andreasskeidsvoll]

Description

I am seeing incorrect dgemm_ results when one application thread is inside dsyevd_ and another application thread calls dgemm_, using the same loaded OpenBLAS library in the same process.

This originally surfaced in jax-ml/jax#38056, where JAX async dispatch could leave a CPU LAPACK eigensolver call in flight while a later NumPy/OpenBLAS matmul started. The C reproducer below removes JAX, NumPy, and SciPy. It links directly against OpenBLAS and calls the Fortran BLAS/LAPACK symbols.

Environment

OpenBLAS config: OpenBLAS 0.3.33 NO_AFFINITY USE_OPENMP VORTEX MAX_THREADS=128
Platform: macOS / Apple ARM64
Install source: conda-forge

Summary

The reproducer does the following:

background pthread:  dsyevd_ on its own matrix/work arrays
main thread:         repeatedly calls dgemm_ while dsyevd_ is in flight

The dgemm_ inputs are all ones, so every output entry should equal K.

The user arrays for the eigensolver and GEMM are separate. The only intended shared object is the OpenBLAS library/runtime itself.

I used OMP_NUM_THREADS to control the thread count because this OpenBLAS build reports USE_OPENMP.

Observed behavior

This fails for me:

OMP_NUM_THREADS=8 ./openblas_dsyevd_dgemm_repro 0 64 64 64 128 50

Output from one run:

OpenBLAS config: OpenBLAS 0.3.33 NO_AFFINITY USE_OPENMP VORTEX MAX_THREADS=128
OpenBLAS num_threads: 8
use_mutex=0 eig_n=64 M=64 N=64 K=128 iterations=50 GEMM_work=524288
FAIL iter=1 gemm_iter=1 max_abs_err=96.4034 bad=128 first_bad_index=1056 first_bad_value=31.7662 expected=128
Observed failure after 10 GEMM calls.

With a mutex around both OpenBLAS calls, the same case passes:

OMP_NUM_THREADS=8 ./openblas_dsyevd_dgemm_repro 1 64 64 64 128 50

Output:

OpenBLAS config: OpenBLAS 0.3.33 NO_AFFINITY USE_OPENMP VORTEX MAX_THREADS=128
OpenBLAS num_threads: 8
use_mutex=1 eig_n=64 M=64 N=64 K=128 iterations=50 GEMM_work=524288
No incorrect result observed after 50 GEMM calls.

Starting the process with one OpenMP thread also passes:

OMP_NUM_THREADS=1 ./openblas_dsyevd_dgemm_repro 0 64 64 64 128 50

Output:

OpenBLAS config: OpenBLAS 0.3.33 NO_AFFINITY USE_OPENMP VORTEX MAX_THREADS=128
OpenBLAS num_threads: 1
use_mutex=0 eig_n=64 M=64 N=64 K=128 iterations=50 GEMM_work=524288
No incorrect result observed after 314 GEMM calls.

I also see the same threshold behavior for several OpenMP thread counts. With M=N=64, K=128 fails for OMP_NUM_THREADS=2,3,4,8, while the adjacent K=127 case passes for all tested thread counts. With OMP_NUM_THREADS=1, both cases pass.

Threshold-like behavior

The failing case above has:

eig_n = 64
M * N * K = 64 * 64 * 128 = 524288 = 2^19

I also see failures for other factorizations of the same GEMM work:

OMP_NUM_THREADS=8 ./openblas_dsyevd_dgemm_repro 0 64 64 64 128 50
OMP_NUM_THREADS=8 ./openblas_dsyevd_dgemm_repro 0 64 64 128 64 50
OMP_NUM_THREADS=8 ./openblas_dsyevd_dgemm_repro 0 64 128 64 64 50

Outputs from one run:

use_mutex=0 eig_n=64 M=64 N=64 K=128 iterations=50 GEMM_work=524288
FAIL iter=1 gemm_iter=1 max_abs_err=96.4034 bad=128 first_bad_index=1056 first_bad_value=31.7662 expected=128
Observed failure after 10 GEMM calls.

use_mutex=0 eig_n=64 M=64 N=128 K=64 iterations=50 GEMM_work=524288
FAIL iter=16 gemm_iter=0 max_abs_err=64.1916 bad=192 first_bad_index=7424 first_bad_value=47.9644 expected=64
Observed failure after 104 GEMM calls.

use_mutex=0 eig_n=64 M=128 N=64 K=64 iterations=50 GEMM_work=524288
FAIL iter=1 gemm_iter=1 max_abs_err=64.1916 bad=1152 first_bad_index=64 first_bad_value=-0.137265 expected=64
Observed failure after 10 GEMM calls.

Adjacent cases just below the apparent thresholds pass in my tests:

OMP_NUM_THREADS=8 ./openblas_dsyevd_dgemm_repro 0 63 64 64 128 50
OMP_NUM_THREADS=8 ./openblas_dsyevd_dgemm_repro 0 64 63 64 128 50
OMP_NUM_THREADS=8 ./openblas_dsyevd_dgemm_repro 0 64 64 63 128 50
OMP_NUM_THREADS=8 ./openblas_dsyevd_dgemm_repro 0 64 64 64 127 50
OMP_NUM_THREADS=8 ./openblas_dsyevd_dgemm_repro 0 64 64 127 64 50
OMP_NUM_THREADS=8 ./openblas_dsyevd_dgemm_repro 0 64 127 64 64 50

Outputs from one run:

use_mutex=0 eig_n=63 M=64 N=64 K=128 iterations=50 GEMM_work=524288
No incorrect result observed after 380 GEMM calls.

use_mutex=0 eig_n=64 M=63 N=64 K=128 iterations=50 GEMM_work=516096
No incorrect result observed after 465 GEMM calls.

use_mutex=0 eig_n=64 M=64 N=63 K=128 iterations=50 GEMM_work=516096
No incorrect result observed after 516 GEMM calls.

use_mutex=0 eig_n=64 M=64 N=64 K=127 iterations=50 GEMM_work=520192
No incorrect result observed after 539 GEMM calls.

use_mutex=0 eig_n=64 M=64 N=127 K=64 iterations=50 GEMM_work=520192
No incorrect result observed after 368 GEMM calls.

use_mutex=0 eig_n=64 M=127 N=64 K=64 iterations=50 GEMM_work=520192
No incorrect result observed after 495 GEMM calls.

On this machine/build, the issue appears to show up when both of these are true:

eig_n >= 64
M * N * K >= 2^19

under concurrent dsyevd_ / dgemm_ entry with multiple OpenMP threads.

Expected behavior

I would not expect dgemm_ to silently return incorrect values when it is operating on independent user arrays.

I realize concurrent calls into a multithreaded BLAS library from an already-multithreaded application might be discouraged or unsupported. Still, the failure mode here is silent incorrect output from dgemm_, not just oversubscription or poor performance. If this usage pattern is unsupported for USE_OPENMP builds, it would be useful to know what the intended safe usage pattern is.

Reproducer

// openblas_dsyevd_dgemm_repro.c
//
// Reproducer for incorrect dgemm_ results when dgemm_ is called while
// another application thread is inside dsyevd_.
//
//   thread 1: dsyevd_
//   thread 0: dgemm_ while dsyevd_ is in flight
//
// GEMM computes ones(M,K) @ ones(K,N), so every C entry should equal K.
//
// Usage:
//   ./openblas_dsyevd_dgemm_repro [use_mutex] [eig_n] [M] [N] [K] [iterations]
//
// Default:
//   use_mutex = 0
//   eig_n     = 64
//   M         = 64
//   N         = 64
//   K         = 128
//   iterations = 50
//
// For the default GEMM, M*N*K = 64*64*128 = 524288 = 2^19.
//
// Example failing case:
//   OMP_NUM_THREADS=8 ./openblas_dsyevd_dgemm_repro 0
//
// Mutex control:
//   OMP_NUM_THREADS=8 ./openblas_dsyevd_dgemm_repro 1
//
// Single-thread OpenMP control:
//   OMP_NUM_THREADS=1 ./openblas_dsyevd_dgemm_repro 0
//
// Adjacent below-threshold case:
//   OMP_NUM_THREADS=8 ./openblas_dsyevd_dgemm_repro 0 64 64 64 127 50

#define _GNU_SOURCE

#include <math.h>
#include <pthread.h>
#include <stdatomic.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

extern char *openblas_get_config(void);
extern int openblas_get_num_threads(void);

extern void dgemm_(
    const char *TRANSA, const char *TRANSB,
    const int *M, const int *N, const int *K,
    const double *ALPHA,
    const double *A, const int *LDA,
    const double *B, const int *LDB,
    const double *BETA,
    double *C, const int *LDC
);

extern void dsyevd_(
    const char *JOBZ, const char *UPLO,
    const int *N,
    double *A, const int *LDA,
    double *W,
    double *WORK, const int *LWORK,
    int *IWORK, const int *LIWORK,
    int *INFO
);

static pthread_mutex_t blas_mutex = PTHREAD_MUTEX_INITIALIZER;
static atomic_int eig_entered = 0;
static atomic_int eig_done = 0;
static int use_mutex = 0;

static void *xalloc(size_t nbytes) {
    void *p = NULL;
    if (posix_memalign(&p, 64, nbytes) != 0 || p == NULL) {
        fprintf(stderr, "allocation failed for %zu bytes\n", nbytes);
        exit(1);
    }
    return p;
}

static void fill(double *x, size_t n, double v) {
    for (size_t i = 0; i < n; ++i) {
        x[i] = v;
    }
}

static void make_symmetric(double *A, int n) {
    for (int j = 0; j < n; ++j) {
        for (int i = 0; i < n; ++i) {
            A[i + (size_t)j * n] =
                (i == j) ? (double)n + 1.0
                         : 1.0 / (1.0 + fabs((double)i - (double)j));
        }
    }
}

typedef struct {
    int n;
    double *A0;
    double *A;
    double *W;
    double *WORK;
    int *IWORK;
    int LWORK;
    int LIWORK;
    int INFO;
} EigCtx;

static void *eig_thread(void *arg) {
    EigCtx *e = (EigCtx *)arg;

    char jobz = 'V';
    char uplo = 'U';
    int lda = e->n;

    memcpy(e->A, e->A0, sizeof(double) * (size_t)e->n * (size_t)e->n);
    e->INFO = 0;

    if (use_mutex) {
        pthread_mutex_lock(&blas_mutex);
    }

    atomic_store_explicit(&eig_entered, 1, memory_order_release);

    dsyevd_(
        &jobz, &uplo,
        &e->n,
        e->A, &lda,
        e->W,
        e->WORK, &e->LWORK,
        e->IWORK, &e->LIWORK,
        &e->INFO
    );

    if (use_mutex) {
        pthread_mutex_unlock(&blas_mutex);
    }

    atomic_store_explicit(&eig_done, 1, memory_order_release);

    return NULL;
}

static void call_dgemm(double *A, double *B, double *C, int M, int N, int K) {
    char trans = 'N';
    double alpha = 1.0;
    double beta = 0.0;

    if (use_mutex) {
        pthread_mutex_lock(&blas_mutex);
    }

    dgemm_(
        &trans, &trans,
        &M, &N, &K,
        &alpha,
        A, &M,
        B, &K,
        &beta,
        C, &M
    );

    if (use_mutex) {
        pthread_mutex_unlock(&blas_mutex);
    }
}

static int check_gemm(double *C, int M, int N, int K, int iter, int gemm_iter) {
    const double expected = (double)K;
    const size_t total = (size_t)M * (size_t)N;

    double max_abs_err = 0.0;
    size_t bad = 0;
    size_t first_bad = 0;
    double first_value = 0.0;

    for (size_t i = 0; i < total; ++i) {
        double err = fabs(C[i] - expected);

        if (err > max_abs_err) {
            max_abs_err = err;
        }

        if (err > 1e-9) {
            if (bad == 0) {
                first_bad = i;
                first_value = C[i];
            }
            ++bad;
        }
    }

    if (bad) {
        fprintf(stderr,
                "FAIL iter=%d gemm_iter=%d max_abs_err=%g bad=%zu "
                "first_bad_index=%zu first_bad_value=%g expected=%g\n",
                iter, gemm_iter, max_abs_err, bad,
                first_bad, first_value, expected);
        return 1;
    }

    return 0;
}

int main(int argc, char **argv) {
    use_mutex = (argc > 1) ? atoi(argv[1]) : 0;

    int eig_n      = (argc > 2) ? atoi(argv[2]) : 64;
    int M          = (argc > 3) ? atoi(argv[3]) : 64;
    int N          = (argc > 4) ? atoi(argv[4]) : 64;
    int K          = (argc > 5) ? atoi(argv[5]) : 128;
    int iterations = (argc > 6) ? atoi(argv[6]) : 50;

    if (eig_n < 1 || M < 1 || N < 1 || K < 1 || iterations < 1) {
        fprintf(stderr,
                "Usage: %s [use_mutex] [eig_n] [M] [N] [K] [iterations]\n",
                argv[0]);
        return 2;
    }

    fprintf(stderr, "OpenBLAS config: %s\n", openblas_get_config());
    fprintf(stderr, "OpenBLAS num_threads: %d\n", openblas_get_num_threads());
    fprintf(stderr,
            "use_mutex=%d eig_n=%d M=%d N=%d K=%d iterations=%d "
            "GEMM_work=%zu\n",
            use_mutex, eig_n, M, N, K, iterations,
            (size_t)M * (size_t)N * (size_t)K);

    EigCtx e;
    e.n = eig_n;
    e.INFO = 0;

    e.A0 = xalloc(sizeof(double) * (size_t)eig_n * (size_t)eig_n);
    e.A  = xalloc(sizeof(double) * (size_t)eig_n * (size_t)eig_n);
    e.W  = xalloc(sizeof(double) * (size_t)eig_n);

    make_symmetric(e.A0, eig_n);

    char jobz = 'V';
    char uplo = 'U';
    int lda = eig_n;

    int lwork_query = -1;
    int liwork_query = -1;
    double work_query = 0.0;
    int iwork_query = 0;

    memcpy(e.A, e.A0, sizeof(double) * (size_t)eig_n * (size_t)eig_n);

    dsyevd_(
        &jobz, &uplo,
        &eig_n,
        e.A, &lda,
        e.W,
        &work_query, &lwork_query,
        &iwork_query, &liwork_query,
        &e.INFO
    );

    if (e.INFO != 0) {
        fprintf(stderr, "dsyevd workspace query failed: info=%d\n", e.INFO);
        return 1;
    }

    e.LWORK = (int)work_query;
    e.LIWORK = iwork_query;

    e.WORK = xalloc(sizeof(double) * (size_t)e.LWORK);
    e.IWORK = xalloc(sizeof(int) * (size_t)e.LIWORK);

    double *GA = xalloc(sizeof(double) * (size_t)M * (size_t)K);
    double *GB = xalloc(sizeof(double) * (size_t)K * (size_t)N);
    double *GC = xalloc(sizeof(double) * (size_t)M * (size_t)N);

    fill(GA, (size_t)M * (size_t)K, 1.0);
    fill(GB, (size_t)K * (size_t)N, 1.0);

    long total_gemm = 0;

    for (int iter = 0; iter < iterations; ++iter) {
        pthread_t th;

        atomic_store_explicit(&eig_entered, 0, memory_order_release);
        atomic_store_explicit(&eig_done, 0, memory_order_release);

        if (pthread_create(&th, NULL, eig_thread, &e) != 0) {
            fprintf(stderr, "pthread_create failed\n");
            return 1;
        }

        while (atomic_load_explicit(&eig_entered, memory_order_acquire) == 0) {
        }

        int gemm_iter = 0;

        while (atomic_load_explicit(&eig_done, memory_order_acquire) == 0) {
            memset(GC, 0, sizeof(double) * (size_t)M * (size_t)N);

            call_dgemm(GA, GB, GC, M, N, K);

            ++total_gemm;

            if (check_gemm(GC, M, N, K, iter, gemm_iter)) {
                pthread_join(th, NULL);
                fprintf(stderr, "Observed failure after %ld GEMM calls.\n", total_gemm);

                free(e.A0);
                free(e.A);
                free(e.W);
                free(e.WORK);
                free(e.IWORK);
                free(GA);
                free(GB);
                free(GC);

                return 1;
            }

            ++gemm_iter;
        }

        pthread_join(th, NULL);

        if (e.INFO != 0) {
            fprintf(stderr, "dsyevd failed: info=%d\n", e.INFO);
            return 1;
        }
    }

    fprintf(stderr, "No incorrect result observed after %ld GEMM calls.\n", total_gemm);

    free(e.A0);
    free(e.A);
    free(e.W);
    free(e.WORK);
    free(e.IWORK);
    free(GA);
    free(GB);
    free(GC);

    return 0;
}

Compile command I used:

clang -O2 -std=c11 -pthread openblas_dsyevd_dgemm_repro.c \
  "$CONDA_PREFIX/lib/libopenblas.0.dylib" \
  -Wl,-rpath,"$CONDA_PREFIX/lib" \
  -lm \
  -o openblas_dsyevd_dgemm_repro

[andreasskeidsvoll]

I also reproduced this in a fresh conda environment with only conda-forge OpenBLAS and the required build tools installed:

conda create -n oblas-min -c conda-forge openblas clang make
conda activate oblas-min

I then compiled the reproducer with:

clang -O2 -Wall -Wextra \
  openblas_dsyevd_dgemm_repro.c \
  -I"$CONDA_PREFIX/include" \
  -L"$CONDA_PREFIX/lib" \
  -lopenblas \
  -Wl,-rpath,"$CONDA_PREFIX/lib" \
  -o openblas_dsyevd_dgemm_repro

otool -L gives:

./openblas_dsyevd_dgemm_repro:
        @rpath/libopenblas.0.dylib (compatibility version 0.0.0, current version 0.0.0)
        /usr/lib/libSystem.B.dylib (compatibility version 1.0.0, current version 1356.0.0)

The failing case still fails in this environment:

OMP_NUM_THREADS=8 ./openblas_dsyevd_dgemm_repro 0 64 64 64 128 50

OpenBLAS config: OpenBLAS 0.3.33 NO_AFFINITY USE_OPENMP VORTEX MAX_THREADS=128
OpenBLAS num_threads: 8
use_mutex=0 eig_n=64 M=64 N=64 K=128 iterations=50 GEMM_work=524288
FAIL iter=0 gemm_iter=1 max_abs_err=96.4034 bad=608 first_bad_index=32 first_bad_value=31.7445 expected=128
Observed failure after 2 GEMM calls.

The mutex and single-thread controls pass in the same environment:

OMP_NUM_THREADS=8 ./openblas_dsyevd_dgemm_repro 1 64 64 64 128 50
OMP_NUM_THREADS=1 ./openblas_dsyevd_dgemm_repro 0 64 64 64 128 50

OpenBLAS config: OpenBLAS 0.3.33 NO_AFFINITY USE_OPENMP VORTEX MAX_THREADS=128
OpenBLAS num_threads: 8
use_mutex=1 eig_n=64 M=64 N=64 K=128 iterations=50 GEMM_work=524288
No incorrect result observed after 50 GEMM calls.

OpenBLAS config: OpenBLAS 0.3.33 NO_AFFINITY USE_OPENMP VORTEX MAX_THREADS=128
OpenBLAS num_threads: 1
use_mutex=0 eig_n=64 M=64 N=64 K=128 iterations=50 GEMM_work=524288
No incorrect result observed after 252 GEMM calls.

[martin-frbg]

Not immediately reproducible on M1 (nor AMD Zen5) with OpenBLAS built from source. How many cores does your hardware have ?

[andreasskeidsvoll]

The machine has 14 CPU cores: 10 performance cores and 4 efficiency cores.

One possibly relevant difference is that I am using the conda-forge OpenBLAS build rather than a local source build:

OpenBLAS 0.3.33 NO_AFFINITY USE_OPENMP VORTEX MAX_THREADS=128

I also reduced the reproducer further. DSYEVD/DSYTRD are not required; the issue reproduces for me with concurrent DSYR2K and DGEMM only:

background thread: repeated DSYR2K('U','N'), syr_n=32, syr_k=32
foreground thread: repeated DGEMM('N','N'), M=64, N=64, K=128

The reduced case fails for me with OMP_NUM_THREADS=8, and also starts failing already at OMP_NUM_THREADS=2. The controls pass with OMP_NUM_THREADS=1, with a mutex around the BLAS calls, when DSYR2K is reduced to syr_n=31, syr_k=32, or when DGEMM is reduced to K=127.

I can attach the reduced reproducer if useful.

[andreasskeidsvoll]

@martin-frbg Small update: I also tested local OpenBLAS v0.3.33 builds on the same machine, using the reduced DSYR2K/DGEMM reproducer: openblas_dsyr2k_dgemm_minimal_repro.c

The local pthreads/non-OpenMP builds pass for both targets:

OpenBLAS 0.3.33 NO_AFFINITY VORTEX MAX_THREADS=128
OPENBLAS_NUM_THREADS=8
PASS

OpenBLAS 0.3.33 NO_AFFINITY ARMV8 MAX_THREADS=128
OPENBLAS_NUM_THREADS=8
PASS

The local OpenMP builds fail for both targets:

OpenBLAS 0.3.33 NO_AFFINITY USE_OPENMP VORTEX MAX_THREADS=128
OMP_NUM_THREADS=8
FAIL

OpenBLAS 0.3.33 NO_AFFINITY USE_OPENMP ARMV8 MAX_THREADS=128
OMP_NUM_THREADS=8
FAIL

For the OpenMP builds, the mutex and OMP_NUM_THREADS=1 controls pass. So in this local build matrix, the failure appears to track USE_OPENMP, not the target-specific kernel path.

[martin-frbg]

Still can't reproduce this (with either reproducer, and now also including a Neoverse N1 host for more cores).
Which compiler version(s) are you using ? (Any chance the reproducer on your end might be linking to both GNU libgomp and LLVM libomp ?)

[andreasskeidsvoll]

@martin-frbg I checked the compiler/OpenMP-runtime side more carefully.

The failing local clang/LLVM OpenMP build was made with:

clang version 22.1.6
Target: arm64-apple-darwin25.5.0

GNU Fortran (conda-forge gcc 15.2.0-19) 15.2.0

For the failing repro_openmp_vortex executable, I do not see evidence that both GNU libgomp and LLVM libomp are linked or loaded. otool -L shows only LLVM libomp:

./repro_openmp_vortex:
    @rpath/libomp.dylib
    /usr/lib/libSystem.B.dylib

and DYLD_PRINT_LIBRARIES also shows only:

.../envs/oblas-openmp/lib/libomp.dylib

The final executable has LLVM OpenMP-style symbols (__kmpc_*/omp_*) and no GOMP_* symbols.

One nuance: the full clang+gfortran OpenBLAS archive did contain some GOMP_* symbols, presumably from objects that were not pulled into this minimal DSYR2K/DGEMM reproducer. To remove that ambiguity, I also built a cleaner BLAS-only OpenBLAS v0.3.33 with clang/LLVM OpenMP:

OpenBLAS 0.3.33 NO_LAPACK NO_LAPACKE NO_AFFINITY USE_OPENMP ARMV8 MAX_THREADS=128

This was built with NO_LAPACK=1 and NOFORTRAN=1. In that build, the OpenBLAS archive itself has __kmpc_*/omp_* symbols and no GOMP_* symbols. With dynamic LLVM libomp.dylib, the reduced reproducer still fails.

Finally, I built LLVM libomp.a from source and linked the same clean OpenBLAS archive against it statically. In that case, otool -L shows only:

/usr/lib/libSystem.B.dylib

so no dynamic libomp or libgomp is loaded. This static-libomp executable still reproduces the failure:

OMP_NUM_THREADS=8
OpenBLAS 0.3.33 NO_LAPACK NO_LAPACKE NO_AFFINITY USE_OPENMP ARMV8 MAX_THREADS=128
use_mutex=0 syr_n=32 syr_k=32 M=64 N=64 K=128 calls=50
FAIL iter=0 max_abs_err=0.676 bad=512 first_bad_index=32 first_bad_value=128.116 expected=128
Observed failure after 2 background DSYR2K calls.

The controls still pass with a mutex around the BLAS calls and with OMP_NUM_THREADS=1.

As an additional comparison, I also built OpenBLAS with conda-forge GCC 15.2.0. The archive contains GOMP_* symbols, and when I linked the reproducer with static libgomp.a, otool -L again showed only:

/usr/lib/libSystem.B.dylib

That GCC/static-libgomp executable did not reproduce the failure in the same tests, including longer 500-call runs. So in my tests, the failure seems specific to the clang/LLVM OpenMP path, rather than to all USE_OPENMP builds.

[andreasskeidsvoll]

@martin-frbg I added a public GitHub Actions reproducer with a Linux ARM64 / Linux x86_64 matrix:

https://github.com/andreasskeidsvoll/openblas-bug-repro

Latest run:

https://github.com/andreasskeidsvoll/openblas-bug-repro/actions/runs/26885914875

ARM64 job:

https://github.com/andreasskeidsvoll/openblas-bug-repro/actions/runs/26885914875/job/79297965229

x86_64 job:

https://github.com/andreasskeidsvoll/openblas-bug-repro/actions/runs/26885914875/job/79297965219

Both jobs build OpenBLAS v0.3.33 from source as a clean BLAS-only clang/LLVM OpenMP build, using NO_LAPACK=1 and NOFORTRAN=1. Both final executables link to LLVM libomp.so and show LLVM OpenMP-style symbols (__kmpc_*/omp_*) without GOMP_* symbols.

On the ARM64 runner, OpenBLAS reports:

OpenBLAS 0.3.33 NO_LAPACK NO_LAPACKE NO_AFFINITY USE_OPENMP ARMV8 MAX_THREADS=128

and the threaded no-mutex correctness stress test fails:

OMP_NUM_THREADS=8
OpenBLAS num_threads: 4
use_mutex=0 syr_n=32 syr_k=32 M=64 N=64 K=128 calls=1000
FAIL iter=327 max_abs_err=16.944 bad=608 first_bad_index=32 first_bad_value=138.224 expected=128
Observed failure after 272 background DSYR2K calls.

On the x86_64 runner, OpenBLAS reports:

OpenBLAS 0.3.33 NO_LAPACK NO_LAPACKE NO_AFFINITY USE_OPENMP ZEN MAX_THREADS=128

and the same 10 × 1000-call threaded no-mutex stress test does not reproduce the failure.

The workflow is written as a correctness test: the positive controls are run first, and then the threaded no-mutex case is run as an ordinary correctness stress test. The controls include:

OMP_NUM_THREADS=8 with a mutex around the BLAS calls
OMP_NUM_THREADS=1 without a mutex
DSYR2K below threshold: syr_n=31, syr_k=32
DGEMM below threshold: K=127

Those controls pass.

So this no longer looks specific to my local macOS/Apple Silicon setup. It reproduces on Linux ARM64 with a clean clang/LLVM OpenMP OpenBLAS build, while the corresponding Linux x86_64 job passes in this workflow.

[martin-frbg]

Reproduced on M4 with LLVM 22.1.7, but have no immediate suggestion for a fix. The distinction between x86_64 and arm64 is probably that the latter has weak memory ordering (suggesting that the power platform may also be affected). Perhaps trying builds with earlier releases of libomp will provide a hint at a defect there, but it could also be that the build using the GNU implementation is working more by accident.