AArch64 NRM2 returns infinity for a norm below DBL_MAX #6094

Activity

1. BenKnill commented on Oct 5, 2026

   

   BenKnill

   on Oct 5, 2026

   Author

   More actions

   A candidate fix is available in this draft PR.
2. martin-frbg commented on Oct 9, 2026

   

   martin-frbg

   on Oct 9, 2026

   Collaborator

   More actions

   Waiting for a response on the related issue in Reference-LAPACK, but strictly speaking this is not a bug in OpenBLAS when it matches the behavior of the reference implementation. (Might be interesting to see if other implementations like BLIS, MKL or Accelerate fare better)
3. BenKnill commented on Oct 10, 2026

   

   BenKnill

   on Oct 10, 2026

   Author

   More actions

   Here is a quick comparison. It uses the two vectors from this report and the n=8 vector from Reference-LAPACK/lapack#1436, single-threaded, in round-to-nearest. In the table, "max" means DBL_MAX or FLT_MAX, and "max−1ulp" means one ulp below it.

   Implementation	dnrm2 n=2	dznrm2 n=1	snrm2 n=2	scnrm2 n=1	dnrm2 n=8
   OpenBLAS 0.3.32 arm64 (ARMV8, NEOVERSEN1)	inf	inf	max	inf	max
   OpenBLAS #6116 e8b92a8 (base b366de2), HASWELL	max	max	max	max	max
   Reference BLAS 3.12.1 (Ubuntu aarch64)	max−1ulp	max−1ulp	max	max	inf
   BLIS a8037cf (haswell)	inf	inf	max	max	max
   MKL 2026.1.0, default dispatch	max	max	max	max	max
   MKL 2026.1.0, MKL_CBWR=COMPATIBLE	inf	inf	max	max	max
   Accelerate (macOS 27.0.1, M5 Pro)	max	max−1ulp	max	max	max

   The x86 rows ran on an i7-10750H. Both MKL_CBWR=AVX2 and MKL_ENABLE_INSTRUCTIONS=AVX2 matched default dispatch on these five calls; both MKL_CBWR=SSE4_2 and MKL_ENABLE_INSTRUCTIONS=SSE4_2 matched COMPATIBLE.

   The tested OpenBLAS arm64 build and Reference BLAS return Inf on different inputs. In every case here, the exact norm is below DBL_MAX or FLT_MAX; rounding it to nearest gives that largest finite value. BLIS and MKL's compatible mode return Inf on the same double inputs as the tested OpenBLAS arm64 build. The OpenBLAS Haswell build, MKL's default dispatch and Accelerate return finite results for all five calls.

   On the Reference-scaling direction you mentioned in #6095: Reference 3.12.1 is finite on this report's inputs but returns Inf on the lapack#1436 vector, so its current scaling still leaves an avoidable-overflow case. I'm happy to leave #6095 parked while that upstream discussion proceeds.

   Test program

   #include <math.h>
   #include <stdio.h>
   #include <cblas.h>   /* <mkl.h> for MKL; <Accelerate/Accelerate.h> with ACCELERATE_NEW_LAPACK */
   
   static void show(const char *name, double r) {
     printf("  %-26s %-24a %s\n", name, r, isfinite(r) ? "finite" : "NONFINITE");
   }
   
   int main(void) {
     const double d[2] = {0x1.a681399499ad2p+1023, 0x1.2131e6d880433p+1023};
     const float f[2] = {0x1.800008p+127f, 0x1.52a7eep+127f};
     const double v8[8] = {0x1.bd2788e9e8b9ap+1023, 0x1.cb33a7ccd84f4p+1021,
                           0x1.52c1bd33cfe1fp+1021, 0x1.fb2a0fcb66de1p+1021,
                           0x1.6c278b240b52dp+1021, 0x1.1f0c45e85aa4dp+1021,
                           0x1.27cf1d741ec94p+1021, 0x1.73f734386dff4p+1021};
     show("dnrm2 n=2", cblas_dnrm2(2, d, 1));
     show("dznrm2 n=1", cblas_dznrm2(1, (const void *)d, 1));
     show("snrm2 n=2", (double)cblas_snrm2(2, f, 1));
     show("scnrm2 n=1", (double)cblas_scnrm2(1, (const void *)f, 1));
     show("dnrm2 n=8", cblas_dnrm2(8, v8, 1));
     return 0;
   }

   On Debian and Ubuntu, force the reference library at runtime (for example with LD_LIBRARY_PATH) and verify the loaded path with ldd. An explicit link path alone can still record the libblas.so.3 SONAME and load OpenBLAS through alternatives.