dscal: 0 * nan ≠ nan * 0 and 0 * inf ≠ inf * 0

[davidcl]

Here another issue related to #4726 and #4413, impacting dscal this time.

With OpenBLAS 0.3.27:

$ LD_LIBRARY_PATH=lib/thirdparty:lib/thirdparty/redist ./openblas_dscal
result of  nan *  inf =  nan
result of  inf *  nan =  nan
result of    0 *  nan =    0
result of  nan *    0 =  nan
result of    0 *  inf =    0
result of  inf *    0 = -nan

With refBLAS 3.12.0 :

$ LD_LIBRARY_PATH=$HOME/work/dependencies/lapack/build/lib ./openblas_dscal
result of  nan *  inf =  nan
result of  inf *  nan =  nan
result of    0 *  nan =  nan
result of  nan *    0 =  nan
result of    0 *  inf = -nan
result of  inf *    0 = -nan

test code

#include <stdio.h>
#include <math.h>

extern void dscal_(int*, double*, double*, int*);

int main(int argc, char ** argv)
{
    int ONE = 1;

    double scalar = NAN;
    double inout = INFINITY;

    printf("result of %4.5g * %4.5g = ", scalar, inout);
    dscal_(&ONE, &scalar, &inout, &ONE);
    printf("%4.5g\n", inout);

    scalar = INFINITY;
    inout = NAN;

    printf("result of %4.5g * %4.5g = ", scalar, inout);
    dscal_(&ONE, &scalar, &inout, &ONE);
    printf("%4.5g\n", inout);

    scalar = 0;
    inout = NAN;

    printf("result of %4.5g * %4.5g = ", scalar, inout);
    dscal_(&ONE, &scalar, &inout, &ONE);
    printf("%4.5g\n", inout);

    scalar = NAN;
    inout = 0;

    printf("result of %4.5g * %4.5g = ", scalar, inout);
    dscal_(&ONE, &scalar, &inout, &ONE);
    printf("%4.5g\n", inout);

    scalar = 0;
    inout = INFINITY;

    printf("result of %4.5g * %4.5g = ", scalar, inout);
    dscal_(&ONE, &scalar, &inout, &ONE);
    printf("%4.5g\n", inout);

    scalar = INFINITY;
    inout = 0;

    printf("result of %4.5g * %4.5g = ", scalar, inout);
    dscal_(&ONE, &scalar, &inout, &ONE);
    printf("%4.5g\n", inout);

    return 0;
}

[martin-frbg]

Already suspected that, curious how nobody seems to have hit (or at least bothered reporting) these in all those years...