Complex SYR2/SPR2 negative increments use half the required starting offset and can segfault

[BenKnill]

The Fortran-ABI exports csyr2_, zsyr2_, cspr2_, and zspr2_ can
access memory before their input vectors when given valid negative increments.
With UPLO='U', n=4, INCX=INCY=-1, four complex elements in each of X and Y,
alpha=1+0i, and valid zero-initialized output storage, the supplied guard-page
reproducer terminates each call with SIGSEGV.

The corresponding calls with INCX=INCY=1 return normally. X and Y start at
the beginning of readable mappings, with an inaccessible page immediately
before each mapping. This makes the pre-input access visible as a fault; it
does not imply every ordinary allocation would segfault.

Observed output

The unchanged reproducer below prints this on the independently checked
0.3.21 binary:

csyr2_ positive=PASS negative=SIGSEGV
cspr2_ positive=PASS negative=SIGSEGV
zsyr2_ positive=PASS negative=SIGSEGV
zspr2_ positive=PASS negative=SIGSEGV
SUMMARY: BUG reproduced=4/4 controls=4/4

The independently checked 0.3.30 binary gives the same outcome for all four
exports, with scipy_ prefixed symbol names. Here positive=PASS means only
normal return, not a verified numerical matrix result. The program is a
crash detector; it is not a full correctness or access-bounds regression test.

Source diagnosis and input convention

interface/zsyr2.c and
interface/zspr2.c are shared by the
single- and double-complex builds. Their negative-increment setup advances
FLOAT *x and FLOAT *y by (n-1)*abs(inc) scalar slots. Each complex element
occupies two scalar slots. For this example the required starting offset is
six scalar slots; the wrapper uses three.

The complex reference implementation bundled in OpenBLAS
declares X and Y as complex arrays and uses KX = 1 - (N-1)*INCX and the
corresponding Y expression for negative increments. The caller supplies the
base of the storage array, not an already-adjusted pointer to its last element.
The reproducer follows that convention. This n=4, inc=-1 witness does not
involve a large-integer overflow.

The apparent fix is to apply the complex-element-to-scalar factor consistently
to both X and Y offsets, using appropriately wide arithmetic. This report does
not claim a patch has been tested or that removing the crash alone establishes
numerical correctness. The scope here is these four rank-2 Fortran-ABI exports;
no CBLAS, Hermitian, rank-1, or exploitability claim is being made.

Runtime evidence and source scope

An independent x86-64 Linux check reproduced the result in both of these
existing bundled libraries, using GCC 14.2.0 and OPENBLAS_NUM_THREADS=1:

Library	openblas_get_config()	Reproduction
CasADi-bundled libcasadi-tp-openblas.so	OpenBLAS 0.3.21 NO_AFFINITY CORE2 MAX_THREADS=16	Four of four exports
SciPy-bundled libscipy_openblas-6cdc3b4a.so	OpenBLAS 0.3.30 DYNAMIC_ARCH NO_AFFINITY SkylakeX MAX_THREADS=64	Four of four exports; only the lookup names were changed to add scipy_

These are configuration strings from particular distributed binaries, not
claims about every build of those releases. The original investigation also
reported the same result in AArch64 builds of v0.3.34 (e0166008be8e466242aa76b2ff75ce3f0fbf574a) and
31e82fa8c509e6f0d96288de3c20d8916d894e72. Those AArch64 runs were not independently repeated in the x86-64
check; the physical runner/operating-system details are not inferred here.

The relevant source was inspected at v0.3.34 and 31e82fa8c509e6f0d96288de3c20d8916d894e72. At the time of
review, develop was 277ef7801416412ee970873701fdf14d603d2c3f; its intervening WASM packing commit does not
change these interface files. The current-head statement is a source check,
not a claim of running a binary built from that commit.

Reproduce on Linux

Use an LP64 OpenBLAS shared library with 32-bit BLAS integers and unprefixed
symbols, compiled for the same architecture as the test. Do not use this
hard-coded int call interface against an ILP64 / INTERFACE64=1 library.
No Fortran compiler is needed to compile this C caller of the library's
Fortran-ABI exports.

cc -O2 -Wall -Wextra -Werror -o repro repro.c -ldl
OPENBLAS_NUM_THREADS=1 ./repro /absolute/path/to/libopenblas.so

The program forks one child per call. On affected builds the negative-stride
children intentionally fault. Exit 0 means all four target crashes reproduced
and all four positive controls returned normally. Exit 1 is incomplete
reproduction, not a verified fix; exit 2 is a detected setup error.

Complete repro.c

/* Isolated direct-call reproduction for the four complex SYR2/SPR2 exports.
 * A guard before each complex input makes a wrong negative-stride adjustment
 * fail deterministically while the positive-stride control stays valid. */
#define _GNU_SOURCE
#include <dlfcn.h>
#include <signal.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/mman.h>
#include <sys/wait.h>
#include <unistd.h>

typedef void (*csyr2_fn)(char *, int *, float *, float *, int *, float *, int *, float *, int *);
typedef void (*cspr2_fn)(char *, int *, float *, float *, int *, float *, int *, float *);
typedef void (*zsyr2_fn)(char *, int *, double *, double *, int *, double *, int *, double *, int *);
typedef void (*zspr2_fn)(char *, int *, double *, double *, int *, double *, int *, double *);

static void *guarded_input(size_t bytes) {
    size_t page = (size_t)sysconf(_SC_PAGESIZE);
    if (bytes > page) exit(2);
    unsigned char *map = mmap(NULL, 2 * page, PROT_NONE,
                              MAP_PRIVATE | MAP_ANONYMOUS, -1, 0);
    if (map == MAP_FAILED || mprotect(map + page, page, PROT_READ | PROT_WRITE)) exit(2);
    return map + page;
}

static void invoke(void *symbol, int variant, int inc) {
    char upper = 'U'; int n = 4, lda = 4;
    if (variant < 2) {
        float alpha[2] = {1.0f, 0.0f};
        float *x = guarded_input(8 * sizeof(float));
        float *y = guarded_input(8 * sizeof(float));
        float a[32] = {0}, ap[20] = {0};
        for (int i = 0; i < 4; i++) { x[2*i] = i + 1; y[2*i] = i + 2; }
        if (variant == 0) ((csyr2_fn)symbol)(&upper, &n, alpha, x, &inc, y, &inc, a, &lda);
        else ((cspr2_fn)symbol)(&upper, &n, alpha, x, &inc, y, &inc, ap);
    } else {
        double alpha[2] = {1.0, 0.0};
        double *x = guarded_input(8 * sizeof(double));
        double *y = guarded_input(8 * sizeof(double));
        double a[32] = {0}, ap[20] = {0};
        for (int i = 0; i < 4; i++) { x[2*i] = i + 1; y[2*i] = i + 2; }
        if (variant == 2) ((zsyr2_fn)symbol)(&upper, &n, alpha, x, &inc, y, &inc, a, &lda);
        else ((zspr2_fn)symbol)(&upper, &n, alpha, x, &inc, y, &inc, ap);
    }
}

static int run_case(void *symbol, int variant, int inc) {
    pid_t child = fork();
    if (child < 0) return -1;
    if (!child) { invoke(symbol, variant, inc); _exit(0); }
    int status;
    if (waitpid(child, &status, 0) != child) return -1;
    if (WIFSIGNALED(status)) return -WTERMSIG(status);
    return WIFEXITED(status) ? WEXITSTATUS(status) : -1;
}

int main(int argc, char **argv) {
    if (argc != 2) { fprintf(stderr, "usage: repro AARCH64_OPENBLAS_SO\n"); return 2; }
    setenv("OPENBLAS_NUM_THREADS", "1", 1);
    void *library = dlopen(argv[1], RTLD_NOW | RTLD_LOCAL);
    if (!library) { fprintf(stderr, "dlopen: %s\n", dlerror()); return 2; }
    const char *names[] = {"csyr2_", "cspr2_", "zsyr2_", "zspr2_"};
    int controls = 0, reproduced = 0;
    for (int i = 0; i < 4; i++) {
        void *symbol = dlsym(library, names[i]);
        if (!symbol) { fprintf(stderr, "missing: %s\n", names[i]); return 2; }
        int positive = run_case(symbol, i, 1);
        int negative = run_case(symbol, i, -1);
        printf("%s positive=%s negative=%s\n", names[i],
               positive == 0 ? "PASS" : "ERROR",
               negative == -SIGSEGV ? "SIGSEGV" : "NO_SIGSEGV");
        controls += positive == 0;
        reproduced += positive == 0 && negative == -SIGSEGV;
    }
    printf("SUMMARY: %s reproduced=%d/4 controls=%d/4\n",
           reproduced == 4 && controls == 4 ? "BUG" : "INCOMPLETE",
           reproduced, controls);
    return reproduced == 4 && controls == 4 ? 0 : 1;
}

[martin-frbg]

Interesting, this fairly trivial bug (missing the factor of 2 in the recalculation of x or y) appears to have gone unnoticed since the days of GotoBLAS. I'll push a fix tomorrow if I can find the time

[BenKnill]

Thanks! Confirmed on develop commit 277ef7801 (x86-64, TARGET=HASWELL): the guard-page reproducer reports SIGSEGV for all four exports with inc = -1. A numerical result check against a naive reference (n=5, both uplo values, seven (incx, incy) pairs) also fails every tested negative-increment case. Adding the factor of 2 to both adjustments in interface/zsyr2.c and interface/zspr2.c locally removes all four SIGSEGV reports and makes the result check pass 56/56, so your diagnosis checks out here. I'll rerun them against your fix when it's up and report back.

[BenKnill]

Tested d1c068401 (PR #6110) on x86-64 (TARGET=HASWELL, GCC 15.2): the guard-page reproducer no longer reports SIGSEGV for the four exports with inc = -1 (reproduced=0/4, positive controls 4/4), and the numerical result check matches the reference in 56/56 cases, including all 40 tested cases with a negative increment. Thanks for the fix!