[cakedev0]

In threaded builds without OpenMP, blas_memory_alloc and blas_memory_free take the global alloc_lock on every call, so all concurrent callers of single-threaded BLAS serialize on it twice per call. This is the case when an application calls small dgemm from the threads of its own OpenMP parallel region (e.g. scikit-learn's KMeans, see scikit-learn/scikit-learn#17334). OpenMP builds don't have this problem: they claim a slot of memory[] with its own lock and release it without any lock.

This PR makes non-OpenMP builds use the same per-slot path for memory[]. alloc_lock is now only taken for the overflow area (newmemory) and on the first call. This covers pthreads builds, Windows builds (where alloc_lock is already a blas_lock spinlock) and USE_THREAD=0 USE_LOCKING=1 builds. Single-threaded builds without USE_LOCKING still take no lock at all.

I think this is the same contention as OpenMathLib#5589 (symv/trsv from many Julia threads, which uses a pthreads build), but I haven't measured symv or trsv.

The global mutex replaced the per-slot lock in OpenMathLib/OpenBLAS@e882f3d6f (2017, OpenMathLib#1299), for data races found by TSan in OpenMathLib#888. At the time, the fields of memory[] were plain volatile, blas_lock was inline assembly on x86 (xchgl), which TSan can't see, and the diagnosis in OpenMathLib#888 was that mixing pthread mutexes with blas_lock doesn't work. Since OpenMathLib#5876 (in 0.3.34), the fields are _Atomic and, with HAVE_C11, blas_lock is a C11 compare-and-swap that TSan understands. With this PR, outside of blas_shutdown, no field of memory[] is protected by alloc_lock anymore: a slot is claimed with its blas_lock, and used and addr are only read and written atomically. It's the code OpenMP builds already run on every architecture, and the linux_thread_sanitizer CI job will check it for pthreads builds. Without HAVE_C11, the slot is released with WMB and a volatile store, as in OpenMP builds. The same lock was discussed in OpenMathLib#1782.

With 64 threads each calling a small dgemm (Lloyd layout below, 102400 calls), pthread_mutex_lock is called 204933 times on develop (31% contended; 64% with the scipy-openblas wheel) and 132 times with this PR.

Benchmarks

Xeon 6787P, one socket (86 cores), TARGET=SAPPHIRERAPIDS USE_OPENMP=0 NUM_THREADS=64, GCC 16.1 libgomp with GOMP_SPINCOUNT=300000, OPENBLAS_NUM_THREADS=1. Each parallel region does a fixed number of dgemm("T", "N", k, 256, d, ...) calls (one per chunk of 256 samples, as in k-means). I picked the TN layout on purpose: with M * N > 1200, the Sapphire Rapids small-matrix kernel doesn't take it, so the call goes through the buffer allocator. Min-max of 3 interleaved runs, ms per region:

threads	k=10, d=3, develop	k=10, d=3, this PR	k=10, d=16, develop	k=10, d=16, this PR
1	2.01-2.02	2.03-2.10	4.25	4.32-4.45
8	1.87-1.89	0.38-0.41	1.21-1.53	0.67-0.69
16	3.12-3.19	0.42	3.18-3.42	0.46
32	3.14-3.25	0.42-0.47	3.26-3.34	0.52-0.53
64	3.40-3.46	0.52-0.55	3.35-3.48	0.63-0.64

On develop, these regions get slower than single-threaded beyond 8 threads. With bigger calls the lock matters less: 4.81-4.98 -> 4.00-4.07 ms at 64 threads for k=100, d=50, and 4.08-4.16 -> 3.95-3.96 ms for k=10, d=128 (memory-bound). With 1 thread, this PR is 1-4% slower on the small calls above. I'm not sure why. It could be the seq_cst stores to used.

On hybrid CPUs, the lock also causes stalls when the OpenMP runtime spins (scikit-learn/scikit-learn#17334): a contended mutex inside a parallel region stalls the whole team. Core Ultra 5 225H, 14 threads, TARGET=HASWELL, 20 calls of dgemm("N", "N", 256, 8, 50, ...) per region, GOMP_SPINCOUNT=300000: 4.2-5.0 ms -> 0.07-0.31 ms per region (laptop on battery, so only indicative).

End to end, with the reproducer of scikit-learn/scikit-learn#17334 on a Core Ultra X7 358H (16 threads, GCC 13 libgomp) and a scikit-learn branch where OpenBLAS is the last blocking call in KMeans' parallel regions, I swapped the OpenBLAS bundled with scipy: 7.4-21.6 s with scipy's wheel, 8.3-14.0 s with develop, 1.11-1.16 s with this PR.

Benchmark

// gcc -O2 -fopenmp blas_scaling.c -o blas_scaling -ldl
// OPENBLAS_NUM_THREADS=1 OMP_NUM_THREADS=64 ./blas_scaling libopenblas.so dgemm_ 1360 200 10 3
#include <dlfcn.h>
#include <omp.h>
#include <stdio.h>
#include <stdlib.h>

typedef void (*dgemm_t)(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 *);
#define CHUNK 256

int main(int argc, char **argv) {
    void *lib = dlopen(argv[1], RTLD_NOW | RTLD_GLOBAL);
    if (!lib) { fprintf(stderr, "%s\n", dlerror()); return 1; }
    dgemm_t dgemm = (dgemm_t)dlsym(lib, argv[2]);
    int n_chunks = atoi(argv[3]), n_rep = atoi(argv[4]);
    int k = atoi(argv[5]), d = atoi(argv[6]), m = CHUNK;
    double *X = malloc((size_t)n_chunks * CHUNK * d * sizeof(double));
    double *centers = malloc((size_t)d * k * sizeof(double));
    double *C = malloc((size_t)n_chunks * CHUNK * k * sizeof(double));
    for (size_t i = 0; i < (size_t)n_chunks * CHUNK * d; i++) X[i] = (i % 97) * 1e-2;
    for (int i = 0; i < d * k; i++) centers[i] = (i % 13) * 1e-2;
    double alpha = -2.0, beta = 0.0, t0 = 0;
    for (int rep = 0; rep < 2 * n_rep; rep++) {
        if (rep == n_rep) t0 = omp_get_wtime();
        #pragma omp parallel for schedule(static)
        for (int c = 0; c < n_chunks; c++)
            dgemm("T", "N", &k, &m, &d, &alpha, centers, &d,
                  X + (size_t)c * CHUNK * d, &d, &beta, C + (size_t)c * CHUNK * k, &k);
    }
    printf("%9.1f us/region\n", (omp_get_wtime() - t0) / n_rep * 1e6);
    return 0;
}

Tests

On the Xeon with the first commit: cpp_thread_test (dgemm_tester with 64 and 96 concurrent threads, the latter going through the overflow path, dgemv_tester, dgemm_mixed_tester), openblas_utest and openblas_utest_ext pass. Measured on develop at 3374642 and rebased since (no change to memory.c or common.h in between). The benchmarks above are also from the first commit: the second one only changes the first mapping of a slot and builds without SMP.

With both commits, under TSan (GCC 16.1 with the options of the linux_thread_sanitizer job, NUM_THREADS=32 NUM_PARALLEL=2, so 128 buffers):

• pthreads build: no report on the three tests of the CI job, on dgemm_thread_safety and dgemv_thread_safety with 64 and 200 concurrent callers and OPENBLAS_NUM_THREADS=1, and on dgemm_thread_safety 1024 200 1, which goes through the overflow path.
• USE_THREAD=0 USE_LOCKING=1: one report, on memory_initialized in blas_memory_alloc, whose init isn't locked without SMP. It's the same on develop, so I left it out of this PR.
• USE_THREAD=0 (no TSan): builds and ctest passes.

I couldn't test on ARM or POWER, nor with a compiler without C11 atomics, nor on Windows.

Any other comments?

USE_TLS=1 doesn't avoid the problem: it takes 3 mutexes per call (alloc_lock to read memory_initialized, key_lock twice in get_memory_table). I left it out of this PR, but I can look at it in a follow-up if you think it's worth it.

AI usage disclosure

I used Claude Code to read memory.c and its history, write the patch and the benchmark, and run the measurements.

🤖 Generated with Claude Code

[cakedev0]

Closed in favor of OpenMathLib#6105