Avoid overflow from division in GETF2 potentially causing NaN

[martin-frbg]

for fault seen in numpy/numpy#22025

[mattip]

Would it help to close the numpy issue if I create an openblas-lib build of the develop branch or should I wait for a 0.3.22 release?

[angsch]

Is there a reason to use fabs(temp1) > 1.e-305 rather than SFMIN as used in LAPACK? SFMIN = 2**-1022 = DBL_MIN $\approx 2.22*10^{-308}$ . I can think of a badly scaled matrix (admittedly a pathological case) where LAPACK succeeds, but OpenBLAS fails with the fix. Also, I would use >= rather than >.

[martin-frbg]

No reason except imited testing once the problem/solution was identified. Also seems possible now that the problem only exists for clang/arm64 so could be ifdef'd to skip the fabs call on other platforms.

[angsch]

@martin-frbg Here is an example that illustrates the problem.

A = [ 0  a
      a  0 ], where |a| is a value below 1.e-305

The expected LU factorization with pivoting is

L = [ 1 0 ]  U = [ a 0 ] P = [ 0 1 ]
    [ 0 1 ],     [ 0 a ],    [ 1 0 ]

The expected memory representation is ipiv = [ 2, 2], A = [ a 0 0 a ].

With the patch, the computed result is

L = [ 1 0 ]  U = [ 0 0 ] P = [ 0 1 ]
    [ a 1 ],     [ 0 a ],    [ 1 0 ]

I consider this a bug. It will not be possible to solve a system with that LU factorization because U(1,1) == 0.0 will incur a division by zero.

Copy & paste reproducer:

// gcc test-getf2.c /path/to/libopenblas_skylakexp-r0.3.21.dev.a
#include <stdio.h>
#include <float.h>

extern void dgetf2_(int *m, int *n, double *A, int *ldA, int *ipiv, int *info);

int main() {
  int m = 2, n = 2, ldA = 2, info = 0;
  double a = DBL_MIN; // smallest positive normal number
  int ipiv[2];
  double A[4] = {0.0, a,
                 a  , 0.0};

  dgetf2_(&m, &n, A, &ldA, ipiv, &info);

  printf("U = [%.6e %.6e\n     %.6e %.6e]\n", A[0], A[2], 0.0, A[3]);
  printf("L = [%.6e %.6e\n     %.6e %.6e]\n", 1.0, 0.0, A[1], 1.0);
  printf("ipiv %d, %d", ipiv[0], ipiv[1]);
  return 0;
}

[martin-frbg]

@angsch thanks for reminding me