Reproducer with SciPy's test suite
SciPy's sparse.linalg iterative-solver tests (lgmres, gcrotmk in single precision) fail on linux-aarch64 in some environments: lgmres
makes zero progress and gcrotmk diverges on a well-conditioned (cond≈48) Poisson matrix, because every sdot-based orthogonalization returns garbage. In an AddressSanitizer-instrumented CPython, v0 reliably contains poison-pattern bytes at the call site, making the failure deterministic (4 × float32(0xcbcbcbcb) ≈ -1.07e8 added to every dot product). In ordinary processes v0 often happens to be zero, which is why the bug usually stays hidden — but any caller leaving FP junk in v0 (a caller-saved register!) gets silently wrong results.
Reproducer with a SciPy ASan build from the branch at scipy/scipy#24066, and with pixi installed on a Linux aarch64 machine (for me a container on an M5 Macbook):
# fails reproducibly:
pixi run test-asan -s sparse.linalg -- -k gcrotmk
# edit pixi.toml to add `env.OPENBLAS_CORETYPE = "ARMv8"`,
# then the same test command will pass
ddot does not seem to be affected (the allocator happened to pick a compatible register there), but it uses the same macro pattern and should be audited along with the dsdot/sdsdot variants and the SVE kernels.
OpenBLAS 0.3.33, conda-forge aarch64 build (libopenblas 0.3.33 pthreads_h9d3fd7e_0), DYNAMIC_ARCH, core detected: NEOVERSEN1
kernel/arm64/dot_kernel_asimd.c hasn't been touched since the 0.3.33 release, so the bug still exists on develop.
Reproducing the issue was one thing, but understanding the root cause was quite another, that was much more Claude Fable 5 than me - the analysis below is what it produced. I don't really understand the assembly code that seems to be the problem here, so I won't open a PR with the suggested fix (I can help test it of course). @Mousius would you be able to help with this?
Summary of the bug
On arm64, sdot/cblas_sdot silently return wrong results for unit-stride inputs with n >= 64 whenever the caller happens to leave non-zero data in SIMD register v0. The result returned is dot(x, y) + horizontal_sum_f32(v0_at_entry).
The cause is in kernel/arm64/dot_kernel_asimd.c (still present on develop): the KERNEL_F inline-asm body hardcodes v0–v7 as accumulators,
" fmla v0.4s, v16.4s, v24.4s \n"
" fmla v1.4s, v17.4s, v25.4s \n"
...
but the initialization only zeroes d1–d7 plus the compiler-chosen output operand:
" fmov "OUT", "REG0" \n"
" fmov d1, xzr \n"
...
" fmov d7, xzr \n"
where OUT is the asm operand [DOT_] "=&w" (dot). The code implicitly assumes the register allocator assigns DOT_ to v0. In the conda-forge 0.3.33 aarch64 build (GCC), it assigns v31 instead, so v0 is never zeroed, yet the main loop accumulates into it and the final reduction sums v0..v7 into the result. Disassembly of the shipped binary (sdot_k_NEOVERSEN1 → dot_kernel_asimd):
597f8c: fmov s31, wzr ; OUT allocated to v31, zeroed
597f90: fmov d1, xzr ; v1..v7 zeroed
...
597fa8: fmov d7, xzr
; v0 never initialized!
597fe0: ldp q16, q17, [x1]
597fe4: ldp q24, q25, [x3]
597ff0: fmla v0.4s, v16.4s, v24.4s ; accumulates into garbage
...
598098: fadd v0.4s, v0.4s, v1.4s ; reduction sums the garbage
Independently of the above, the clobber list
: "cc", "memory", "d1", "d2", "d3", "d4", "d5", "d6", "d7"
is missing the hardcoded v0 and the scratch registers v16–v31 that the asm body writes. That is undefined behavior on its own: the compiler is free to keep live values in those registers across the asm block, so surrounding code can also be miscompiled. Whether the bug manifests in any given binary depends on register allocation, i.e. on compiler version and flags.
Reproducer (no sanitizers, plain C)
#include <stdio.h>
extern float sdot_(int*, float*, int*, float*, int*);
int main(void) {
static float x[64], y[64];
for (int i = 0; i < 64; i++) { x[i] = 1.0f; y[i] = 1.0f; }
int n = 64, inc = 1;
float junk[4] = {1e6f, 1e6f, 1e6f, 1e6f};
__asm__ volatile("ld1 {v0.4s}, [%0]" :: "r"(junk) : "v0");
float r = sdot_(&n, x, &inc, y, &inc);
printf("sdot(ones, ones, n=64) = %f (correct: 64.0)\n", r);
return r == 64.0f ? 0 : 1;
}$ gcc repro.c -lopenblas && LD_BIND_NOW=1 ./a.out
sdot(ones, ones, n=64) = 4000064.000000 (correct: 64.0)
(LD_BIND_NOW=1 just keeps the lazy-binding trampoline from touching v0; the bug does not depend on it.)
Results are correct for n <= 63 (scalar path) and for non-unit strides, which use a different path. OPENBLAS_CORETYPE=ARMV8 avoids the affected kernel entirely and is a usable workaround.
Suggested fix
In dot_kernel_asimd.c:
- Zero the hardcoded accumulator explicitly (e.g.
movi v0.16b, #0, or fmov d0, xzr) instead of relying on OUT being allocated to v0, or consistently use the %[DOT_] operand rather than hardcoded registers.
- Add all registers the asm writes to the clobber list:
v0 (or d0) and v16–v31 (d16–d31).
Reproducer with SciPy's test suite
SciPy's
sparse.linalgiterative-solver tests (lgmres,gcrotmkin single precision) fail on linux-aarch64 in some environments:lgmresmakes zero progress and
gcrotmkdiverges on a well-conditioned (cond≈48) Poisson matrix, because everysdot-based orthogonalization returns garbage. In an AddressSanitizer-instrumented CPython,v0reliably contains poison-pattern bytes at the call site, making the failure deterministic (4 × float32(0xcbcbcbcb) ≈ -1.07e8added to every dot product). In ordinary processesv0often happens to be zero, which is why the bug usually stays hidden — but any caller leaving FP junk inv0(a caller-saved register!) gets silently wrong results.Reproducer with a SciPy ASan build from the branch at scipy/scipy#24066, and with
pixiinstalled on a Linux aarch64 machine (for me a container on an M5 Macbook):ddotdoes not seem to be affected (the allocator happened to pick a compatible register there), but it uses the same macro pattern and should be audited along with thedsdot/sdsdotvariants and the SVE kernels.OpenBLAS 0.3.33, conda-forge aarch64 build (
libopenblas 0.3.33 pthreads_h9d3fd7e_0),DYNAMIC_ARCH, core detected:NEOVERSEN1kernel/arm64/dot_kernel_asimd.chasn't been touched since the 0.3.33 release, so the bug still exists ondevelop.Reproducing the issue was one thing, but understanding the root cause was quite another, that was much more Claude Fable 5 than me - the analysis below is what it produced. I don't really understand the assembly code that seems to be the problem here, so I won't open a PR with the suggested fix (I can help test it of course). @Mousius would you be able to help with this?
Summary of the bug
On arm64,
sdot/cblas_sdotsilently return wrong results for unit-stride inputs withn >= 64whenever the caller happens to leave non-zero data in SIMD registerv0. The result returned isdot(x, y) + horizontal_sum_f32(v0_at_entry).The cause is in
kernel/arm64/dot_kernel_asimd.c(still present ondevelop): theKERNEL_Finline-asm body hardcodesv0–v7as accumulators,but the initialization only zeroes
d1–d7plus the compiler-chosen output operand:where
OUTis the asm operand[DOT_] "=&w" (dot). The code implicitly assumes the register allocator assignsDOT_tov0. In the conda-forge 0.3.33 aarch64 build (GCC), it assignsv31instead, sov0is never zeroed, yet the main loop accumulates into it and the final reduction sumsv0..v7into the result. Disassembly of the shipped binary (sdot_k_NEOVERSEN1→dot_kernel_asimd):Independently of the above, the clobber list
is missing the hardcoded
v0and the scratch registersv16–v31that the asm body writes. That is undefined behavior on its own: the compiler is free to keep live values in those registers across the asm block, so surrounding code can also be miscompiled. Whether the bug manifests in any given binary depends on register allocation, i.e. on compiler version and flags.Reproducer (no sanitizers, plain C)
(
LD_BIND_NOW=1just keeps the lazy-binding trampoline from touchingv0; the bug does not depend on it.)Results are correct for
n <= 63(scalar path) and for non-unit strides, which use a different path.OPENBLAS_CORETYPE=ARMV8avoids the affected kernel entirely and is a usable workaround.Suggested fix
In
dot_kernel_asimd.c:movi v0.16b, #0, orfmov d0, xzr) instead of relying onOUTbeing allocated tov0, or consistently use the%[DOT_]operand rather than hardcoded registers.v0(ord0) andv16–v31(d16–d31).