Possible regression in NEOVERSEV1 kernel(s) after 0.3.29

[martin-frbg]

Our CI was running with OpenBlas 0.3.29. I skipped updating to 0.3.30 because of build failures (possibly this #5459, and I skipped 0.3.31 because it was not packaged with Spack. I now updated OpenBlas to 0.3.32 and now some regression tests are failing on armv8.4a (numbers do not match reference results). I tried with openblas@develop to include #5710, but the regression tests are still failing.

This is on a c7g AWS instance, with gcc 13 and with llvm 19.

Thank you in advance

Originally posted by @Sbozzolo in #5708

[martin-frbg]

@Sbozzolo 0.3.31 had several updates to SGER/DGER, SAXPY/DAXPY, 0.3.30 added GEMV and SYMV kernels , it would help if you could identify BLAS functions involved in the failing tests, or isolate a failing test from your CI.
(For testing, it should be possible to simply comment out (#) the GER and AXPY entries in kernel/arm64/KERNEL.NEOVERSEV1 and rebuild OpenBLAS, and/or to swap in the SGEMVNKERNEL/SGEMVTKERNEL declarations from KERNEL.NEOVERSEN1 - SYMV kernel entries would need to be changed (i.e. borrowed from NEOVERSEN1) in kernel/arm64/KERNEL.ARMV8SVE )

[martin-frbg]

Sbozzolo commented on Mar 31, 2026
Sbozzolo
on Mar 31, 2026

@martin-frbg, thank you for your suggestions.

I think I tried your suggestions of commenting out/swapping the kernels, but I still found failures. So, I bisected OpenBlas versions and found that 0.3.30 was already leading to failures. I git bisected between 0.3.29 and 0.3.30 and found that f1e628b seems to be the culprit.

I'll try digging further, but I don't know how much I will be able to do.

This was my git bisect script (with Spack and Palace):

# Install directly into the spack openblas prefix
source /data/home/gbozzola/spack/share/spack/setup-env.sh
spack env activate /data/home/gbozzola/prod_palace_neoverse_v1_openblas
INSTALLDIR=$(spack location -i openblas)

echo "=== Building OpenBLAS at $(git rev-parse --short HEAD) ==="
echo "=== Installing to $INSTALLDIR ==="

# Clean and build
make -j$(nproc) TARGET=NEOVERSEV1 clean 2>/dev/null || true
make -j$(nproc) TARGET=NEOVERSEV1 2>&1 | tail -3

if [ ${PIPESTATUS[0]} -ne 0 ]; then
    echo "BUILD FAILED — skipping this commit"
    exit 125  # git bisect skip
fi

make TARGET=NEOVERSEV1 PREFIX=$INSTALLDIR install 2>&1 | tail -3

# Source modules
source /etc/profile.d/modules.sh 2>/dev/null || source /usr/share/modules/init/bash 2>/dev/null || true
module load julia

echo "=== Running Palace dipole test ==="
cd /data/home/gbozzola/palace/test/examples

julia --project runtests.jl --test-cases "antenna/antenna_short_dipole" 2>&1 | tail -20
RESULT=${PIPESTATUS[0]}

if [ $RESULT -eq 0 ]; then
    echo "=== PASS (good commit) ==="
    exit 0
else
    echo "=== FAIL (bad commit) ==="
    exit 1
fi

[Sbozzolo]

I looked into this more carefully and found that the failure goes away if I change gemv_n_sve_v1x3.c from using interleaved operations to a sequential one (on the most recent develop and on 0.3.30). The code below was generated with the help of an LLM. I do not understand why this happens.

/***************************************************************************
Copyright (c) 2025, The OpenBLAS Project
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:

   1. Redistributions of source code must retain the above copyright
      notice, this list of conditions and the following disclaimer.

   2. Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in
      the documentation and/or other materials provided with the
      distribution.
   3. Neither the name of the OpenBLAS project nor the names of
      its contributors may be used to endorse or promote products
      derived from this software without specific prior written
      permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*****************************************************************************/

#include <arm_sve.h>

#include "common.h"

#ifdef DOUBLE
#define SV_COUNT svcntd
#define SV_TYPE svfloat64_t
#define SV_TRUE svptrue_b64
#define SV_WHILE svwhilelt_b64_s64
#define SV_DUP svdup_f64
#else
#define SV_COUNT svcntw
#define SV_TYPE svfloat32_t
#define SV_TRUE svptrue_b32
#define SV_WHILE svwhilelt_b32_s64
#define SV_DUP svdup_f32
#endif

int CNAME(BLASLONG m, BLASLONG n, BLASLONG dummy1, FLOAT alpha, FLOAT *a,
          BLASLONG lda, FLOAT *x, BLASLONG inc_x, FLOAT *y, BLASLONG inc_y,
          FLOAT *buffer)
{
  BLASLONG i;
  BLASLONG ix, iy;
  BLASLONG j;
  FLOAT *a_ptr;
  FLOAT temp;

  a_ptr = a;
  ix = 0;

  if (inc_y == 1) {
    BLASLONG sve_size = SV_COUNT();
    svbool_t pg_true = SV_TRUE();

    /* Process 3 consecutive columns at a time: (j, j+1, j+2) */
    for (j = 0; j + 2 < n; j += 3) {
      SV_TYPE temp0_vec = SV_DUP(alpha * x[ix]);
      SV_TYPE temp1_vec = SV_DUP(alpha * x[ix + inc_x]);
      SV_TYPE temp2_vec = SV_DUP(alpha * x[ix + inc_x * 2]);

      FLOAT *a0 = a_ptr;
      FLOAT *a1 = a_ptr + lda;
      FLOAT *a2 = a_ptr + lda * 2;

      i = 0;
      while ((i + sve_size - 1) < m) {
        SV_TYPE y0_vec = svld1(pg_true, y + i);

        SV_TYPE a00_vec = svld1(pg_true, a0 + i);
        SV_TYPE a01_vec = svld1(pg_true, a1 + i);
        SV_TYPE a02_vec = svld1(pg_true, a2 + i);

        y0_vec = svmla_m(pg_true, y0_vec, temp0_vec, a00_vec);
        y0_vec = svmla_m(pg_true, y0_vec, temp1_vec, a01_vec);
        y0_vec = svmla_m(pg_true, y0_vec, temp2_vec, a02_vec);

        svst1(pg_true, y + i, y0_vec);
        i += sve_size;
      }

      if (i < m) {
        svbool_t pg = SV_WHILE(i, m);

        SV_TYPE y0_vec = svld1(pg, y + i);

        SV_TYPE a00_vec = svld1(pg, a0 + i);
        SV_TYPE a01_vec = svld1(pg, a1 + i);
        SV_TYPE a02_vec = svld1(pg, a2 + i);

        y0_vec = svmla_m(pg, y0_vec, temp0_vec, a00_vec);
        y0_vec = svmla_m(pg, y0_vec, temp1_vec, a01_vec);
        y0_vec = svmla_m(pg, y0_vec, temp2_vec, a02_vec);

        svst1(pg, y + i, y0_vec);
      }

      a_ptr += lda * 3;
      ix += inc_x * 3;
    }

    /* Cleanup: remaining 1 or 2 columns */
    for (; j < n; j++) {
      SV_TYPE temp_vec = SV_DUP(alpha * x[ix]);

      i = 0;
      while ((i + sve_size - 1) < m) {
        SV_TYPE y_vec = svld1(pg_true, y + i);
        SV_TYPE a_vec = svld1(pg_true, a_ptr + i);
        y_vec = svmla_m(pg_true, y_vec, temp_vec, a_vec);
        svst1(pg_true, y + i, y_vec);
        i += sve_size;
      }

      if (i < m) {
        svbool_t pg = SV_WHILE(i, m);
        SV_TYPE y_vec = svld1(pg, y + i);
        SV_TYPE a_vec = svld1(pg, a_ptr + i);
        y_vec = svmla_m(pg, y_vec, temp_vec, a_vec);
        svst1(pg, y + i, y_vec);
      }

      a_ptr += lda;
      ix += inc_x;
    }

    return (0);
  }

  for (j = 0; j < n; j++) {
    temp = alpha * x[ix];
    iy = 0;
    for (i = 0; i < m; i++) {
      y[iy] += temp * a_ptr[i];
      iy += inc_y;
    }
    a_ptr += lda;
    ix += inc_x;
  }
  return (0);
}

[Sbozzolo]

This could also be a reflection of the fact that our set up might be very sensitive to numerical instabilities.

That said, with the above fix, our solution matches ARMPL extremely well. Otherwise, the error is very big.

[martin-frbg]

Hmm. Maybe if your test case handles a large amount of very small numbers, or numbers of widely different magnitudes in single precision , order of operations plays a role through intermediate rounding ??
@iha-taisei or @annop-w can you help with an educated guess ?

[Sbozzolo]

Hmm. Maybe if your test case handles a large amount of very small numbers, or numbers of widely different magnitudes in single precision , order of operations plays a role through intermediate rounding ??

I checked the following:

• Reverted only DGEMVNKERNEL (SGEMVN uses interleaved kernel) -> PASS
• Reverted only SGEMVNKERNEL (DGEMVN uses interleaved kernel) -> FAIL

So, I think the problem is in the double precision kernel.

[annop-w]

@Sbozzolo would you be able to extract M and N for the failing case ? That would be easier for me to reproduce. Thank you.

[Sbozzolo]

@Sbozzolo would you be able to extract M and N for the failing case ? That would be easier for me to reproduce. Thank you.

I am attaching a reproducer that I extracted from my failing application. The reproducer passes with the non-interleaved version.

dgemv_reproducer.txt
test_reproducer.c

[martin-frbg]

I'm a bit unsure what to make of this - running your reproducer on a machine without SVE, I get the same output of

Total mismatches: 128 / 128
Max absolute error: 8.701839e-11 at i=106
Max |y_ref|:        2.182807e+01
Error / (n*eps*max_ref): 270.8

Sanity: max|y_test - y_kernel_dumped| = 4.048495e-11
Sanity: max|y_naive - y_ref_dumped|    = 8.701839e-11

for the NEON gemv_n.S and the simple generic C implementation from kernel/arm/gemv_n.c (and also with 0.3.28)

and running on Zen5 with AVX512:

Total mismatches: 128 / 128
Max absolute error: 5.252421e-11 at i=106
Max |y_ref|:        2.182807e+01
Error / (n*eps*max_ref): 163.5

Sanity: max|y_test - y_kernel_dumped| = 4.678036e-11
Sanity: max|y_naive - y_ref_dumped|    = 8.701839e-11

which seems to be in the same ballpack as the "failed" result from the interleaved SVE kernel. I'm wondering if this might suggest that your non-interleaved "correction" is essentially just reordering intermediate rounding operations to make an inherently fragile calculation yield the desired result in one particular case.

[Sbozzolo]

I tend to agree. The largest Error / (n*eps*max_ref) I see on Neoverse V1 is 768.8, but given the similar ratios you report on other platforms, it's unclear whether that's evidence of an actual problem. It may indeed just be a different floating-point accumulation order.

That said, I should note that this is the only case where we see disagreement. OpenBLAS 0.3.29, MKL, AMDBlis, and ARMPL all produce consistent results for our test suite across a variety of architectures. But, I would understand it if you decided that this is not worth changing (and we will change our test instead).

[martin-frbg]

As far as I can tell the only change since 0.3.29 was further loop unrolling, which I assume would both change accumulation order and be expected to happen automatically with more advanced compilers at some point. But I'm still hoping for the author of the relevant PRs to comment

[martin-frbg]

On Cortex-X2 (ARMV8SVE target, which uses the same DGEMV kernel), I'm getting

Total mismatches: 128 / 128
Max absolute error: 4.048495e-11 at i=90
Max |y_ref|:        2.182807e+01
Error / (n*eps*max_ref): 126.0

Sanity: max|y_test - y_kernel_dumped| = 0.000000e+00
Sanity: max|y_naive - y_ref_dumped|    = 0.000000e+00

(difference to the 768.8 you say you're seeing may be different compiler (clang) and
different libc, this being a phone running Android (and termux)

It does pass with the Reference BLAS:

Total mismatches: 0 / 128
Max absolute error: 0.000000e+00 at i=-1
Max |y_ref|:        2.182807e+01
Error / (n*eps*max_ref): 0.0

Sanity: max|y_test - y_kernel_dumped| = 5.553069e-11
Sanity: max|y_naive - y_ref_dumped|    = 8.701839e-11

but also fails for OpenBLAS with TARGET=NEHALEM

Total mismatches: 128 / 128
Max absolute error: 2.393374e-11 at i=72
Max |y_ref|:        2.182807e+01
Error / (n*eps*max_ref): 74.5

Sanity: max|y_test - y_kernel_dumped| = 3.527401e-11
Sanity: max|y_naive - y_ref_dumped|    = 8.701839e-11

and HASWELL (same results already seen with 0.3.29):

Total mismatches: 128 / 128
Max absolute error: 8.196199e-11 at i=106
Max |y_ref|:        2.182807e+01
Error / (n*eps*max_ref): 255.1

Sanity: max|y_test - y_kernel_dumped| = 4.159362e-11
Sanity: max|y_naive - y_ref_dumped|    = 8.701839e-11

and even OpenBLAS-0.3.20 built for the venerable Prescott:

Total mismatches: 128 / 128
Max absolute error: 4.335554e-11 at i=106
Max |y_ref|:        2.182807e+01
Error / (n*eps*max_ref): 134.9

Sanity: max|y_test - y_kernel_dumped| = 3.913048e-11
Sanity: max|y_naive - y_ref_dumped|    = 8.701839e-11

also on POWER8:

Total mismatches: 128 / 128
Max absolute error: 8.701839e-11 at i=106
Max |y_ref|:        2.182807e+01
Error / (n*eps*max_ref): 270.8

Sanity: max|y_test - y_kernel_dumped| = 4.048495e-11
Sanity: max|y_naive - y_ref_dumped|    = 8.701839e-11

[annop-w]

Sorry for the delay. This is simply round-off errors due to different order of summation and has nothing to do with SVE.
I have tried converting the SVE loop to just scalar

  if (inc_y == 1) {                                                                                                                              
    BLASLONG width = n / 3; // Only process full 3-column blocks                                                                                 
                                                                                                                                                 
    FLOAT *a0_ptr = a_ptr + lda * width * 0;                                                                                                     
    FLOAT *a1_ptr = a_ptr + lda * width * 1;                                                                                                     
    FLOAT *a2_ptr = a_ptr + lda * width * 2;                                                                                                     
                                                                                                                                                 
    FLOAT *x0_ptr = x + inc_x * width * 0;                                                                                                       
    FLOAT *x1_ptr = x + inc_x * width * 1;                                                                                                       
    FLOAT *x2_ptr = x + inc_x * width * 2;                                                                                                       
                                                                                                                                                 
    for (j = 0; j < width; j++) {                                                                                                                
      FLOAT temp0 = alpha * x0_ptr[ix];                                                                                                          
      FLOAT temp1 = alpha * x1_ptr[ix];                                                                                                          
      FLOAT temp2 = alpha * x2_ptr[ix];                                                                                                          
      for (i = 0; i < m; i++) {                                                                                                                  
        y[i] += a0_ptr[i] * temp0;                                                                                                               
        y[i] += a1_ptr[i] * temp1;                                                                                                               
        y[i] += a2_ptr[i] * temp2;                                                                                                               
      }                                                                                                                                          
      a0_ptr += lda;                                                                                                                             
      a1_ptr += lda;                                                                                                                             
      a2_ptr += lda;                                                                                                                             
      ix += inc_x;                                                                                                                               
    }                                                                                                                                            
    for (j = width*3; j < n; j++) {                                                                                                              
      FLOAT *a_col = a + j * lda;                                                                                                                
      FLOAT temp = temp = alpha * x[j * inc_x];                                                                                                  
      for (i = 0; i < m; i++) {                                                                                                                  
        y[i] += a_col[i] * temp;                                                                                                                 
      }                                                                                                                                          
    }                                                                                                                                            
                                                                                                                                                 
    return(0);                                                                                                                                   
  }

and this also fails (max. difference at different i). Of course, if I modify to width = n / 1 so that we process consecutive columns in the kernel above, the test passes. That is also why the SVE kernel processing consecutive columns above also passes.