Skip to content

Commit

Permalink
Merge pull request OpenMathLib#4419 from martin-frbg/issue4413
Browse files Browse the repository at this point in the history 
 [WIP] Add fixes and utests for ZSCAL with NaN or Inf arguments
  • Loading branch information
@martin-frbg
martin-frbg authored Jan 12, 2024
2 parents 3599f2d + 20413ee commit f31bea0
Show file tree
Hide file tree
Showing 15 changed files with 172 additions and 42 deletions.
2 changes: 1 addition & 1 deletion kernel/arm64/zscal.S
View file
Original file line number Diff line number Diff line change
Expand Up @@ -223,7 +223,7 @@ zscal_begin:
fcmp DA_I, #0.0
beq .Lzscal_kernel_RI_zero

b .Lzscal_kernel_R_zero
// b .Lzscal_kernel_R_zero

.Lzscal_kernel_R_non_zero:

Expand Down
6 changes: 4 additions & 2 deletions kernel/mips/KERNEL.P5600
View file
Original file line number Diff line number Diff line change
Expand Up @@ -103,8 +103,10 @@ endif
ifdef HAVE_MSA
SSCALKERNEL = ../mips/sscal_msa.c
DSCALKERNEL = ../mips/dscal_msa.c
CSCALKERNEL = ../mips/cscal_msa.c
ZSCALKERNEL = ../mips/zscal_msa.c
#CSCALKERNEL = ../mips/cscal_msa.c
#ZSCALKERNEL = ../mips/zscal_msa.c
CSCALKERNEL = ../mips/zscal.c
ZSCALKERNEL = ../mips/zscal.c
else
SSCALKERNEL = ../mips/scal.c
DSCALKERNEL = ../mips/scal.c
Expand Down
8 changes: 6 additions & 2 deletions kernel/mips/zscal.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,6 +47,7 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da_r,FLOAT da_i, F
else
{
temp = - da_i * x[ip+1] ;
if (isnan(x[ip]) || isinf(x[ip])) temp = NAN;
x[ip+1] = da_i * x[ip] ;
}
}
Expand All @@ -63,8 +64,11 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da_r,FLOAT da_i, F
x[ip+1] = da_r * x[ip+1] + da_i * x[ip] ;
}
}
x[ip] = temp;

if ( da_r != da_r )
x[ip] = da_r;
else
x[ip] = temp;

ip += inc_x2;
}

Expand Down
1 change: 1 addition & 0 deletions kernel/riscv64/zscal.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,6 +60,7 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da_r,FLOAT da_i, F
else
{
temp = - da_i * x[ip+1] ;
if (isnan(x[ip]) || isinf(x[ip])) temp = NAN;
x[ip+1] = da_i * x[ip] ;
}
}
Expand Down
2 changes: 2 additions & 0 deletions kernel/riscv64/zscal_vector.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -80,6 +80,7 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da_r,FLOAT da_i, F
j += gvl;
ix += inc_x * 2 * gvl;
}
#if 0
}else if(da_r == 0.0){
gvl = VSETVL(n);
BLASLONG stride_x = inc_x * 2 * sizeof(FLOAT);
Expand All @@ -97,6 +98,7 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da_r,FLOAT da_i, F
j += gvl;
ix += inc_xv;
}
#endif
if(j < n){
gvl = VSETVL(n-j);
v0 = VLSEV_FLOAT(&x[ix], stride_x, gvl);
Expand Down
2 changes: 1 addition & 1 deletion kernel/x86/zscal.S
View file
Original file line number Diff line number Diff line change
Expand Up @@ -98,7 +98,7 @@
fcomip %st(1), %st
ffreep %st(0)
jne .L30

jp .L30
EMMS

pxor %mm0, %mm0
Expand Down
1 change: 1 addition & 0 deletions kernel/x86/zscal_sse.S
View file
Original file line number Diff line number Diff line change
Expand Up @@ -87,6 +87,7 @@
xorps %xmm7, %xmm7
comiss %xmm0, %xmm7
jne .L100 # Alpha_r != ZERO
jp .L100 # Alpha_r NaN

comiss %xmm1, %xmm7
jne .L100 # Alpha_i != ZERO
Expand Down
1 change: 1 addition & 0 deletions kernel/x86/zscal_sse2.S
View file
Original file line number Diff line number Diff line change
Expand Up @@ -98,6 +98,7 @@
xorps %xmm7, %xmm7
comisd %xmm0, %xmm7
jne .L100
jp .L100

comisd %xmm1, %xmm7
jne .L100
Expand Down
61 changes: 31 additions & 30 deletions kernel/x86_64/zscal.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,7 @@ USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#endif

#include "common.h"

#include <float.h>

#if defined (SKYLAKEX) || defined (COOPERLAKE) || defined (SAPPHIRERAPIDS)
#include "zscal_microk_skylakex-2.c"
Expand Down Expand Up @@ -222,12 +222,10 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da_r, FLOAT da_i,

if ( da_r == 0.0 )
{

BLASLONG n1 = n & -2;

if ( da_i == 0.0 )
{

while(j < n1)
{

Expand All @@ -253,7 +251,6 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da_r, FLOAT da_i,
}
else
{

while(j < n1)
{

Expand Down Expand Up @@ -356,49 +353,59 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da_r, FLOAT da_i,

alpha[0] = da_r;
alpha[1] = da_i;

if ( da_r == 0.0 )
if ( da_i == 0 )
zscal_kernel_8_zero(n1 , alpha , x);
else
zscal_kernel_8_zero_r(n1 , alpha , x);
// zscal_kernel_8_zero_r(n1 , alpha , x);
zscal_kernel_8(n1 , alpha , x);
else
if ( da_i == 0 )
if ( da_i == 0 && da_r == da_r)
zscal_kernel_8_zero_i(n1 , alpha , x);
else
zscal_kernel_8(n1 , alpha , x);

}
i = n1 << 1;
j = n1;
}


if ( da_r == 0.0 )

if ( da_r == 0.0 || da_r != da_r )
{

if ( da_i == 0.0 )
{

FLOAT res=0.0;
if (da_r != da_r) res= da_r;
while(j < n)
{

x[i]=0.0;
x[i+1]=0.0;
x[i]=res;
x[i+1]=res;
i += 2 ;
j++;

}

}
else
else if (da_r < -FLT_MAX || da_r > FLT_MAX) {
while(j < n)
{
x[i]= NAN;
x[i+1] = da_r;
i += 2 ;
j++;

}

} else
{

while(j < n)
{

temp0 = -da_i * x[i+1];
if (x[i] < -FLT_MAX || x[i] > FLT_MAX)
temp0 = NAN;
x[i+1] = da_i * x[i];
x[i] = temp0;
if ( x[i] == x[i]) //preserve NaN
x[i] = temp0;
i += 2 ;
j++;

Expand All @@ -409,28 +416,24 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da_r, FLOAT da_i,
}
else
{

if ( da_i == 0.0 )
if (da_i == 0.0)
{

while(j < n)
{
while(j < n)
{

temp0 = da_r * x[i];
x[i+1] = da_r * x[i+1];
x[i] = temp0;
i += 2 ;
j++;

}

}
}
else
{

while(j < n)
{

temp0 = da_r * x[i] - da_i * x[i+1];
x[i+1] = da_r * x[i+1] + da_i * x[i];
x[i] = temp0;
Expand All @@ -445,5 +448,3 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da_r, FLOAT da_i,

return(0);
}


1 change: 1 addition & 0 deletions kernel/x86_64/zscal_sse2.S
View file
Original file line number Diff line number Diff line change
Expand Up @@ -82,6 +82,7 @@
pxor %xmm15, %xmm15
comisd %xmm0, %xmm15
jne .L100
jp .L100

comisd %xmm1, %xmm15
jne .L100
Expand Down
20 changes: 16 additions & 4 deletions kernel/zarch/zscal.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -233,9 +233,15 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da_r, FLOAT da_i,

while (j < n1) {

temp0 = -da_i * x[i + 1];
if (isnan(x[i]) || isinf(x[i]))
temp0 = NAN;
else
temp0 = -da_i * x[i + 1];
x[i + 1] = da_i * x[i];
x[i] = temp0;
if (isnan(x[i + inc_x]) || isinf(x[i + inc_x]))
temp1 = NAN;
else
temp1 = -da_i * x[i + 1 + inc_x];
x[i + 1 + inc_x] = da_i * x[i + inc_x];
x[i + inc_x] = temp1;
Expand All @@ -246,7 +252,10 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da_r, FLOAT da_i,

while (j < n) {

temp0 = -da_i * x[i + 1];
if (isnan(x[i]) || isinf(x[i]))
temp0 = NAN;
else
temp0 = -da_i * x[i + 1];
x[i + 1] = da_i * x[i];
x[i] = temp0;
i += inc_x;
Expand Down Expand Up @@ -320,7 +329,7 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da_r, FLOAT da_i,
if (da_i == 0)
zscal_kernel_8_zero(n1, x);
else
zscal_kernel_8_zero_r(n1, alpha, x);
zscal_kernel_8(n1, da_r, da_i, x);
else if (da_i == 0)
zscal_kernel_8_zero_i(n1, alpha, x);
else
Expand All @@ -347,7 +356,10 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da_r, FLOAT da_i,

while (j < n) {

temp0 = -da_i * x[i + 1];
if (isnan(x[i]) || isinf(x[i]))
temp0 = NAN;
else
temp0 = -da_i * x[i + 1];
x[i + 1] = da_i * x[i];
x[i] = temp0;
i += 2;
Expand Down
1 change: 1 addition & 0 deletions utest/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@ else ()
test_dsdot.c
test_dnrm2.c
test_swap.c
test_zscal.c
)
endif ()

Expand Down
2 changes: 1 addition & 1 deletion utest/Makefile
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@ UTESTBIN=openblas_utest

include $(TOPDIR)/Makefile.system

OBJS=utest_main.o test_min.o test_amax.o test_ismin.o test_rotmg.o test_axpy.o test_dotu.o test_dsdot.o test_swap.o test_rot.o test_dnrm2.o
OBJS=utest_main.o test_min.o test_amax.o test_ismin.o test_rotmg.o test_axpy.o test_dotu.o test_dsdot.o test_swap.o test_rot.o test_dnrm2.o test_zscal.o
#test_rot.o test_swap.o test_axpy.o test_dotu.o test_dsdot.o test_fork.o

ifneq ($(NO_LAPACK), 1)
Expand Down
56 changes: 56 additions & 0 deletions utest/test_zscal.c
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,56 @@
#include "openblas_utest.h"
#include <cblas.h>
#ifdef BUILD_COMPLEX16

#ifndef NAN
#define NAN 0.0/0.0
#endif
#ifndef INFINITY
#define INFINITY 1.0/0.0
#endif

CTEST(zscal, i_nan)
{
double i[] = {0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1 };
double nan[] = {NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0};
cblas_zscal(9, i, &nan, 1);
ASSERT_TRUE(isnan(nan[0]));
ASSERT_TRUE(isnan(nan[1]));
ASSERT_TRUE(isnan(nan[16]));
ASSERT_TRUE(isnan(nan[17]));
}

CTEST(zscal, nan_i)
{
double i[] = {0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1 };
double nan[] = {NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0};
cblas_zscal(9, &nan, &i, 1);
ASSERT_TRUE(isnan(i[0]));
ASSERT_TRUE(isnan(i[1]));
ASSERT_TRUE(isnan(i[16]));
ASSERT_TRUE(isnan(i[17]));
}

CTEST(zscal, i_inf)
{
double i[] = {0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1 };
double inf[] = {INFINITY, 0, INFINITY,0, INFINITY,0, INFINITY,0, INFINITY,0, INFINITY,0, INFINITY,0, INFINITY,0, INFINITY,0};
cblas_zscal(9, i, &inf, 1);
ASSERT_TRUE(isnan(inf[0]));
ASSERT_TRUE(isinf(inf[1]));
ASSERT_TRUE(isnan(inf[16]));
ASSERT_TRUE(isinf(inf[17]));
}

CTEST(zscal, inf_i)
{
double i[] = {0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1 };
double inf[] = {INFINITY, 0, INFINITY,0, INFINITY,0, INFINITY,0, INFINITY,0, INFINITY,0, INFINITY,0, INFINITY,0, INFINITY,0};
cblas_zscal(9, &inf, &i, 1);
ASSERT_TRUE(isnan(i[0]));
ASSERT_TRUE(isinf(i[1]));
ASSERT_TRUE(isnan(i[16]));
ASSERT_TRUE(isinf(i[17]));
}

#endif
Loading

0 comments on commit f31bea0

Please sign in to comment.