Merge pull request #4727 from martin-frbg/issue4726

Fix another corner case of infinity handling in x86_64 ZSCAL

kernel/arm/zscal.c

@@ -61,7 +61,9 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da_r,FLOAT da_i, F
 			{
 				temp = - da_i * x[ip+1] ;
 				if (isnan(x[ip]) || isinf(x[ip])) temp = NAN;
-				x[ip+1] = da_i * x[ip]  ;
+				if (!isinf(x[ip+1])) 
+					x[ip+1] = da_i * x[ip]  ;
+				else 	x[ip+1] = NAN;
 			}
 		}
 		else

kernel/mips/zscal.c

@@ -48,7 +48,9 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da_r,FLOAT da_i, F
 			{
 				temp = - da_i * x[ip+1] ;
 				if (isnan(x[ip]) || isinf(x[ip])) temp = NAN;
-				x[ip+1] = da_i * x[ip]  ;
+				if (!isinf(x[ip+1]))
+					x[ip+1] = da_i * x[ip]  ;
+				else    x[ip+1] = NAN;
 			}
 		}
 		else
@@ -56,12 +58,16 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da_r,FLOAT da_i, F
 			if ( da_i == 0.0 )
 			{
 				temp    = da_r * x[ip]  ;
-				x[ip+1] = da_r * x[ip+1];
+				if (!isinf(x[ip+1]))
+					x[ip+1] = da_r * x[ip+1];
+				else    x[ip+1] = NAN;
 			}
 			else
 			{
 				temp    = da_r * x[ip]   - da_i * x[ip+1] ;
-				x[ip+1] = da_r * x[ip+1] + da_i * x[ip]   ;
+				if (!isinf(x[ip+1]))
+					x[ip+1] = da_r * x[ip+1] + da_i * x[ip]   ;
+				else    x[ip+1] = NAN;
 			}
 		}
 		if ( da_r != da_r ) 

kernel/riscv64/zscal.c

@@ -61,7 +61,9 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da_r,FLOAT da_i, F
 			{
 				temp = - da_i * x[ip+1] ;
 				if (isnan(x[ip]) || isinf(x[ip])) temp = NAN;
-				x[ip+1] = da_i * x[ip]  ;
+				if (!isinf(x[ip+1])) 
+					x[ip+1] = da_i * x[ip]  ;
+				else 	x[ip+1] = NAN;
 			}
 		}
 		else

kernel/x86_64/zscal.c

@@ -258,13 +258,17 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da_r, FLOAT da_i,
 						temp0	= NAN;
 					else
 						temp0   = -da_i * x[i+1];
-					x[i+1]       =  da_i * x[i];
+					if (!isinf(x[i+1]))
+						x[i+1]  = da_i * x[i];
+					else    x[i+1]	= NAN;     
 					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];
+					if (!isinf(x[i+1+inc_x]))
+						x[i+1+inc_x] = da_i * x[i+inc_x];
+					else    x[i+1+inc_x] = NAN;
 					x[i+inc_x]   =  temp1;
 					i += 2*inc_x ;
 					j+=2;
@@ -278,7 +282,9 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da_r, FLOAT da_i,
 						temp0	= NAN;
 					else
 						temp0   = -da_i * x[i+1];
-					x[i+1]       =  da_i * x[i];
+					if (!isinf(x[i+1]))
+						x[i+1]       =  da_i * x[i];
+					else x[i+1] = NAN;
 					x[i]         =  temp0;
 					i += inc_x ;
 					j++;
@@ -412,7 +418,9 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da_r, FLOAT da_i,
 				temp0        = -da_i * x[i+1];
 				if (isinf(x[i]))
 					temp0 = NAN;
-				x[i+1]       =  da_i * x[i];
+				if (!isinf(x[i+1]))
+					x[i+1]       =  da_i * x[i];
+				else x[i+1] = NAN;
 				if ( x[i] == x[i]) //preserve NaN
 				  x[i]         =  temp0;
 				i += 2 ;

utest/test_zscal.c

@@ -117,4 +117,31 @@ CTEST(zscal, inf_i_inc_2)
     ASSERT_TRUE(isinf(i[17]));
 }
 
+CTEST(zscal, i_0inf)
+{
+    blasint N=9;
+    blasint incX=1;
+    double i[] = {0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1 };
+    double inf[] = {0,INFINITY, 0, INFINITY,0, INFINITY,0, INFINITY,0, INFINITY,0, INFINITY,0, INFINITY,0, INFINITY,0, INFINITY};
+    BLASFUNC(zscal)(&N, i, inf, &incX);
+    ASSERT_TRUE(isinf(inf[0]));
+    ASSERT_TRUE(isnan(inf[1]));
+    ASSERT_TRUE(isinf(inf[16]));
+    ASSERT_TRUE(isnan(inf[17]));
+}
+
+CTEST(zscal, i_0inf_inc_2)
+{
+    blasint N=9;
+    blasint incX=2;
+    double i[] = {0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1 };
+    double inf[] = {0,INFINITY, 0,INFINITY, 0,INFINITY, 0,INFINITY, 0,INFINITY, 0,INFINITY, 0,INFINITY, 0,INFINITY, 0,INFINITY,
+                    0,INFINITY, 0,INFINITY, 0,INFINITY, 0,INFINITY, 0,INFINITY, 0,INFINITY, 0,INFINITY, 0,INFINITY, 0,INFINITY};
+    BLASFUNC(zscal)(&N, i, inf, &incX);
+    ASSERT_TRUE(isinf(inf[0]));
+    ASSERT_TRUE(isnan(inf[1]));
+    ASSERT_TRUE(isinf(inf[16]));
+    ASSERT_TRUE(isnan(inf[17]));
+}
+
 #endif