Merge pull request #4729 from martin-frbg/issue4728

Fix handling of INF or NAN arguments in S/D/C SCAL

kernel/arm64/scal.S

@@ -168,8 +168,8 @@ USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 	cmp	N, xzr
 	ble	.Lscal_kernel_L999
 
-	fcmp	DA, #0.0
-	beq	.Lscal_kernel_zero
+	//fcmp	DA, #0.0
+	//beq	.Lscal_kernel_zero
 
 	cmp	INC_X, #1
 	bne	.Lscal_kernel_S_BEGIN

kernel/mips/dscal_msa.c

@@ -42,7 +42,7 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da, FLOAT *x,
 
     if (1 == inc_x)
     {
-        if (0.0 == da)
+        if (0)     //if (0.0 == da )
         {
             v2f64 zero_v = {0.0, 0.0};
 
@@ -243,9 +243,11 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da, FLOAT *x,
         if (da == 0.0)
         {
             for (i = n; i--;)
-            {
-                *x = 0.0;
-
+            {         
+                if (isfinite(*x))
+                          *x = 0.0;
+                else
+                          *x = NAN;
                 x += inc_x;
             }
         }

kernel/mips/scal.c

@@ -35,7 +35,12 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da, FLOAT *x, BLAS
 	{
 
 		if ( da == 0.0 )
-			x[i]=0.0;
+			if (isnan(x[i])||isinf(x[i]))
+				x[i]=NAN;
+			else
+				x[i]=0.0;
+		else if (isnan(da))
+			x[i]=NAN;
 		else
 			x[i] = da * x[i] ;
 

kernel/mips/sscal_msa.c

@@ -42,7 +42,7 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da, FLOAT *x,
 
     if (1 == inc_x)
     {
-        if (0.0 == da)
+        if (0) // if (0.0 == da)
         {
             v4f32 zero_v = {0.0, 0.0, 0.0, 0.0};
 
@@ -259,7 +259,10 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da, FLOAT *x,
         {
             for (i = n; i--;)
             {
-                *x = 0;
+                if (isfinite(*x))
+                          *x = 0;
+                else
+                          *x = NAN;
                 x += inc_x;
             }
         }

kernel/mips64/scal.S

@@ -79,6 +79,9 @@
 	bc1f	$fcc0, .L50
 	NOP
 
+	bc1t	$fcc0, .L50
+	NOP
+
 	bne	INCX, TEMP, .L20
 	dsra	I, N, 3
 

kernel/power/dscal.c

@@ -73,14 +73,38 @@ static void dscal_kernel_8_zero (BLASLONG n, FLOAT *x)
 
         for( i=0; i<n; i+=8 )
         {
-                x[0] = alpha;
-                x[1] = alpha;
-                x[2] = alpha;
-                x[3] = alpha;
-                x[4] = alpha;
-                x[5] = alpha;
-                x[6] = alpha;
-                x[7] = alpha;
+                if(isfinite(x[0]))
+			x[0] = alpha;
+		else	
+			x[0] = NAN;
+                if(isfinite(x[1]))
+			x[1] = alpha;
+		else
+			x[1] = NAN;
+		if (isfinite(x[2]))
+			x[2] = alpha;
+		else
+			x[2] = NAN;
+		if (isfinite(x[3]))
+	                x[3] = alpha;
+		else
+			x[3] = NAN;
+		if (isfinite(x[4]))
+                	x[4] = alpha;
+		else
+			x[4] = NAN;
+		if (isfinite(x[5]))
+                	x[5] = alpha;
+		else
+			x[5] = NAN;
+		if (isfinite(x[6]))
+                	x[6] = alpha;
+		else
+			x[6] = NAN;
+		if (isfinite(x[7]))
+                	x[7] = alpha;
+		else
+			x[7] = NAN;
                 x+=8;
         }
 
@@ -107,7 +131,10 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da, FLOAT *x, BLAS
 			{
 				BLASLONG align = ((32 - ((uintptr_t)x & (uintptr_t)0x1F)) >> 3) & 0x3;
 				for (j = 0; j < align; j++) {
-					x[j] = 0.0;
+					if (isfinite(x[j]))
+						x[j] = 0.0;
+					else
+						x[j] = NAN;
 				}
 			}
 			BLASLONG n1 = (n-j) & -16;
@@ -127,8 +154,10 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da, FLOAT *x, BLAS
 
 			while(j < n)
 			{
-
-				x[j]=0.0;
+				if (!isfinite(x[j]))
+					x[j]=NAN;
+				else
+					x[j]=0.0;
 				j++;
 			}
 
@@ -176,8 +205,10 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da, FLOAT *x, BLAS
 
 			while(j < n)
 			{
-
-				x[i]=0.0;
+				if (!isfinite(x[i]))
+					x[i]=NAN;
+				else
+					x[i]=0.0;
 				i += inc_x ;
 				j++;
 			}

kernel/power/scal.S

@@ -84,8 +84,9 @@
 	cmpwi	cr0, N, 0
 	blelr-	cr0
 
-	fcmpu	cr0, FZERO, ALPHA
-	bne-	cr0, LL(A1I1)
+//	fcmpu	cr0, FZERO, ALPHA
+//	bne-	cr0, LL(A1I1)
+	b LL(A1I1)
 
 	cmpwi	cr0, INCX, SIZE
 	bne-	cr0, LL(A0IN)

kernel/power/sscal.c

@@ -74,14 +74,38 @@ static void sscal_kernel_16_zero( BLASLONG n, FLOAT *x )
 
         for( i=0; i<n; i+=8 )
         {
-                x[0] = alpha;
-                x[1] = alpha;
-                x[2] = alpha;
-                x[3] = alpha;
-                x[4] = alpha;
-                x[5] = alpha;
-                x[6] = alpha;
-                x[7] = alpha;
+                if (isfinite(x[0]))
+			x[0] = alpha;
+		else
+			x[0] = NAN;
+		if (isfinite(x[1]))
+                	x[1] = alpha;
+		else
+			x[1] = NAN;
+		if (isfinite(x[2]))
+                	x[2] = alpha;
+		else
+			x[2] = NAN;
+		if (isfinite(x[3]))
+                	x[3] = alpha;
+		else
+			x[3] = NAN;
+		if (isfinite(x[4]))
+                	x[4] = alpha;
+		else
+			x[4] = NAN;
+		if (isfinite(x[5]))
+                	x[5] = alpha;
+		else
+			x[5] = NAN;
+		if (isfinite(x[6]))
+                	x[6] = alpha;
+		else
+			x[6] = NAN;
+		if (isfinite(x[7]))
+                	x[7] = alpha;
+		else
+			x[7] = NAN;
                 x+=8;
         }
 
@@ -109,7 +133,10 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da, FLOAT *x, BLAS
 			{
 				BLASLONG align = ((32 - ((uintptr_t)x & (uintptr_t)0x1F)) >> 2) & 0x7;
 				for (j = 0; j < align; j++) {
-					x[j] = 0.0;
+					if (isfinite(x[j]))
+						x[j] = 0.0;
+					else
+						x[j] = NAN;
 				}
 			}
 			BLASLONG n1 = (n-j) & -32;
@@ -129,8 +156,10 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da, FLOAT *x, BLAS
 
 			while(j < n)
 			{
-
-				x[j]=0.0;
+				if (isfinite(x[j]))
+					x[j]=0.0;
+				else
+					x[j]=NAN;
 				j++;
 			}
 
@@ -178,8 +207,10 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da, FLOAT *x, BLAS
 
 			while(j < n)
 			{
-
-				x[i]=0.0;
+				if (isfinite(x[i]))
+					x[i]=0.0;
+				else
+					x[i]=NAN;
 				i += inc_x ;
 				j++;
 			}

kernel/riscv64/scal.c

@@ -48,7 +48,10 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da, FLOAT *x, BLAS
 	{
 
 		if ( da == 0.0 )
-			x[i]=0.0;
+			if (isfinite(x[i]))
+				x[i]=0.0;
+			else
+				x[i]=NAN;
 		else
 			x[i] = da * x[i] ;
 

kernel/riscv64/scal_vector.c

@@ -71,7 +71,7 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da, FLOAT *x, BLAS
         FLOAT_V_T v0, v1;
         unsigned int gvl = 0;
         if(inc_x == 1){
-                if(da == 0.0){
+                if (0){ //if(da == 0.0){
                         memset(&x[0], 0, n * sizeof(FLOAT));
                 }else{
                         gvl = VSETVL(n);
@@ -96,7 +96,7 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da, FLOAT *x, BLAS
                         }
                 }
         }else{
-                if(da == 0.0){
+                if (0) { //if(da == 0.0){
                         BLASLONG stride_x = inc_x * sizeof(FLOAT);
                         BLASLONG ix = 0;
                         gvl = VSETVL(n);

kernel/riscv64/zscal.c

@@ -61,9 +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;
-				if (!isinf(x[ip+1])) 
+				if (!isinf(x[ip+1]))
 					x[ip+1] = da_i * x[ip]  ;
-				else 	x[ip+1] = NAN;
+				else	x[ip+1] = NAN;
 			}
 		}
 		else

kernel/sparc/scal.S

@@ -120,8 +120,10 @@
 
 	FCLR(29)
 
-	FCMP	ALPHA, FZERO
-	fbne	.LL100
+//	FCMP	ALPHA, FZERO
+//	fbne	.LL100
+	b	.LL100
+
 	sll	INCX, BASE_SHIFT, INCX
 
 	cmp	INCX, SIZE

kernel/x86/scal.S

@@ -68,7 +68,8 @@
 	ftst
 	fnstsw	%ax
 	andb	$68, %ah
-	je	.L300		# Alpha != ZERO
+//	je	.L300		# Alpha != ZERO
+	jmp	.L300
 
 /* Alpha == ZERO */
 	cmpl	$1,%esi

kernel/x86_64/cscal.c

@@ -259,11 +259,22 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da_r, FLOAT da_i,
 				while(j < n1)
 				{
 			
+					if (isnan(x[i]) || isinf(x[i]))
+						temp0	= NAN;
+					else
 					temp0        = -da_i * x[i+1];
+					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;
@@ -272,9 +283,14 @@ 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];
+				        if (!isinf(x[i+1]))
 					x[i+1]       =  da_i * x[i];
+					else x[i+1] = NAN;
 					x[i]         =  temp0;
 					i += inc_x ;
 					j++;
@@ -365,42 +381,51 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da_r, FLOAT da_i,
 			else
 				cscal_kernel_16_zero_r(n1 , alpha , x);
 		else
-			if ( da_i == 0 )
-				cscal_kernel_16_zero_i(n1 , alpha , x);
-			else
 				cscal_kernel_16(n1 , alpha , x);
 
 		i = n1 << 1;
 		j = n1;
 	}
 
-
-	if ( da_r == 0.0 )
+	if ( da_r == 0.0 || isnan(da_r) )
 	{
-
 		if ( da_i == 0.0 )
 		{
-
+		FLOAT res=0.0;
+		if (isnan(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 (isinf(da_r)) {
+			while(j < n)
+			{
+					x[i]= NAN;
+					x[i+1] = da_r;
+					i += 2 ;
+					j++;
+
+			}
+
+		} else
 		{
 
 			while(j < n)
 			{
-			
 				temp0        = -da_i * x[i+1];
-				x[i+1]       =  da_i * x[i];
-				x[i]         =  temp0;
+				if (isinf(x[i]))
+					temp0 = NAN;
+				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 ;
 				j++;
 

kernel/x86_64/dscal.c

@@ -169,8 +169,12 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da, FLOAT *x, BLAS
 			while(j < n1)
 			{
 			
-				x[i]=0.0;
-				x[i+inc_x]=0.0;
+				if (isinf(x[i])||isnan(x[i]))
+					x[i]=NAN;
+				else x[i]=0.0;
+				if (isinf(x[i+inc_x])||isnan(x[i+inc_x]))
+					x[i+inc_x]=NAN;
+				else x[i+inc_x]=0.0;
 				i += 2*inc_x ;
 				j+=2;
 
@@ -179,7 +183,9 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da, FLOAT *x, BLAS
 			while(j < n)
 			{
 			
-				x[i]=0.0;
+				if (isinf(x[i])||isnan(x[i]))
+					x[i]=NAN;
+				else x[i]=0.0;
 				i += inc_x ;
 				j++;
 
@@ -213,9 +219,9 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da, FLOAT *x, BLAS
 	BLASLONG n1 = n & -8;
 	if ( n1 > 0 )
 	{
-		if ( da == 0.0 )
-			dscal_kernel_8_zero(n1 , &da , x);
-		else
+//		if ( da == 0.0 )
+//			dscal_kernel_8_zero(n1 , &da , x);
+//		else
 			dscal_kernel_8(n1 , &da , x);
 	}
 
@@ -223,15 +229,24 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da, FLOAT *x, BLAS
 	{
 		for ( i=n1 ; i<n; i++ )
 		{
-			x[i] = 0.0;
+			if(isinf(x[i])||isnan(x[i]))
+				x[i]=NAN;
+			else x[i] = 0.0;
 		}
 	}
+	else if (isinf(da)){
+		for ( i=n1 ; i<n; i++)
+			if (x[i]==0.) x[i]=NAN;
+			else x[i] *=da;
+	}	
 	else
 	{
 
 		for ( i=n1 ; i<n; i++ )
 		{
-			x[i] *= da;
+			if(isinf(x[i]))
+				x[i]=NAN;
+			else x[i] *= da;
 		}
 	}
 	return(0);

kernel/x86_64/scal_sse.S

@@ -76,7 +76,7 @@
 	shufps	$0, %xmm0, %xmm0
 
 	jne	.L100		# Alpha != ZERO
-
+	je	.L100
 /* Alpha == ZERO */
 	cmpq	$SIZE, INCX
 	jne	.L50

kernel/x86_64/scal_sse2.S

@@ -75,7 +75,7 @@
 	comisd	%xmm0, %xmm1
 	jne	.L100		# Alpha != ZERO
 	jp	.L100		# For Alpha = NaN
-
+	je 	.L100		# disable the Alpha=zero path as it does not handle x=inf or nan
 /* Alpha == ZERO */
 	cmpq	$SIZE, INCX
 	jne	.L50

kernel/x86_64/sscal.c

@@ -119,14 +119,16 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da, FLOAT *x, BLAS
 
 		if ( da == 0.0 )
 		{
-
 			BLASLONG n1 = n & -2;
 
 			while(j < n1)
 			{
-
-				x[i]=0.0;
-				x[i+inc_x]=0.0;
+                                if (isinf(x[i])||isnan(x[i]))
+                                        x[i]=NAN;
+                                else x[i]=0.0;
+                                if (isinf(x[i+inc_x])||isnan(x[i+inc_x]))
+                                        x[i+inc_x]=NAN;
+				else	x[i+inc_x]=0.0;
 				i += 2*inc_x ;
 				j+=2;
 
@@ -134,8 +136,9 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da, FLOAT *x, BLAS
 
 			while(j < n)
 			{
-
-				x[i]=0.0;
+                                if (isinf(x[i])||isnan(x[i]))
+                                        x[i]=NAN;
+				else x[i]=0.0;
 				i += inc_x ;
 				j++;
 
@@ -143,7 +146,7 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da, FLOAT *x, BLAS
 		}
 		else
 		{
-
+#if 1
 			BLASLONG n1 = n & -8;
 			if ( n1 > 0 )
 			{
@@ -151,10 +154,9 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da, FLOAT *x, BLAS
 				i = n1 * inc_x;
 				j = n1;
 		        }
-
+#endif
 			while(j < n)
 			{
-
 				x[i] *= da;
 				i += inc_x ;
 				j++;
@@ -162,16 +164,15 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da, FLOAT *x, BLAS
 			}
 
 		}
-
 		return(0);
 	}
 
 	BLASLONG n1 = n & -16;
 	if ( n1 > 0 )
 	{
-		if ( da == 0.0 )
-			sscal_kernel_16_zero(n1 , &da , x);
-		else
+		//if ( da == 0.0 )
+		//	sscal_kernel_16_zero(n1 , &da , x);
+		//else
 			sscal_kernel_16(n1 , &da , x);
 	}
 
@@ -179,7 +180,18 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da, FLOAT *x, BLAS
 	{
 		for ( i=n1 ; i<n; i++ )
 		{
-			x[i] = 0.0;
+                        if (isinf(x[i])||isnan(x[i]))
+                                x[i]=NAN;
+			else x[i]=0.0;
+		}
+	}
+	else if ( isinf(da) )
+	{
+		for ( i=n1 ; i<n; i++ )
+		{
+                        if (x[i] == 0.0)
+                                x[i]=NAN;
+			else x[i] *= da;
 		}
 	}
 	else
@@ -187,7 +199,9 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da, FLOAT *x, BLAS
 
 		for ( i=n1 ; i<n; i++ )
 		{
-			x[i] *= da;
+                        if (isinf(x[i]))
+                                x[i]=NAN;
+			else x[i] *= da;
 		}
 	}
 	return(0);

utest/test_zscal.c

@@ -1,5 +1,449 @@
 #include "openblas_utest.h"
 #include <cblas.h>
+#ifdef BUILD_SINGLE
+
+#ifndef NAN
+#define NAN 0.0/0.0
+#endif
+#ifndef INFINITY
+#define INFINITY 1.0/0.0
+#endif
+
+CTEST(sscal, 0_nan)
+{
+    blasint N=9;
+    blasint incX=1;
+    float i = 0.0;
+    float x[] = {NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN};
+    BLASFUNC(sscal)(&N, &i, x, &incX);
+    ASSERT_TRUE(isnan(x[0]));
+    ASSERT_TRUE(isnan(x[8]));
+}
+
+CTEST(sscal, 0_nan_inc_2)
+{
+    blasint N=9;
+    blasint incX=2;
+    float i = 0.0;
+    float x[] = {NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN,
+                    NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN};
+    BLASFUNC(sscal)(&N, &i, x, &incX);
+    ASSERT_TRUE(isnan(x[0]));
+    ASSERT_TRUE(isnan(x[8]));
+}
+
+CTEST(sscal, nan_0)
+{
+    blasint N=9;
+    blasint incX=1;
+    float i = NAN;
+    float x[] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
+    BLASFUNC(sscal)(&N, &i, x, &incX);
+    ASSERT_TRUE(isnan(x[0]));
+    ASSERT_TRUE(isnan(x[8]));
+}
+
+CTEST(sscal, nan_0_inc_2)
+{
+    blasint N=9;
+    blasint incX=2;
+    float i = NAN;
+    float x[] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+                    0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
+    BLASFUNC(sscal)(&N, &i, x, &incX);
+    ASSERT_TRUE(isnan(x[0]));
+    ASSERT_TRUE(isnan(x[8]));
+}
+
+CTEST(sscal, 0_inf)
+{
+    blasint N=9;
+    blasint incX=1;
+    float i = 0.0;
+    float x[] = {INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY};
+    BLASFUNC(sscal)(&N, &i, x, &incX);
+    ASSERT_TRUE(isnan(x[0]));
+    ASSERT_TRUE(isnan(x[8]));
+}
+
+CTEST(sscal, 0_inf_inc_2)
+{
+    blasint N=9;
+    blasint incX=2;
+    float i = 0.0;
+    float x[] = {INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY,
+                    INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY};
+    BLASFUNC(sscal)(&N, &i, x, &incX);
+    ASSERT_TRUE(isnan(x[0]));
+    ASSERT_TRUE(isnan(x[8]));
+}
+
+CTEST(sscal, inf_0)
+{
+    blasint N=9;
+    blasint incX=1;
+    float i = INFINITY;
+    float x[] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
+    BLASFUNC(sscal)(&N, &i, x, &incX);
+    ASSERT_TRUE(isnan(x[0]));
+    ASSERT_TRUE(isnan(x[8]));
+}
+
+CTEST(sscal, inf_0_inc_2)
+{
+    blasint N=9;
+    blasint incX=2;
+    float i = INFINITY;
+    float x[] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+                    0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
+    BLASFUNC(sscal)(&N, &i, x, &incX);
+    ASSERT_TRUE(isnan(x[0]));
+    ASSERT_TRUE(isnan(x[8]));
+}
+
+CTEST(sscal, nan_inf)
+{
+    blasint N=9;
+    blasint incX=1;
+    float i = NAN;
+    float x[] = {INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY};
+    BLASFUNC(sscal)(&N, &i, x, &incX);
+    ASSERT_TRUE(isnan(x[0]));
+    ASSERT_TRUE(isnan(x[8]));
+}
+
+CTEST(sscal, nan_inf_inc_2)
+{
+    blasint N=9;
+    blasint incX=2;
+    float i = NAN;
+    float x[] = {INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY,
+                    INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY};
+    BLASFUNC(sscal)(&N, &i, x, &incX);
+    ASSERT_TRUE(isnan(x[0]));
+    ASSERT_TRUE(isnan(x[8]));
+}
+
+CTEST(sscal, inf_nan)
+{
+    blasint N=9;
+    blasint incX=1;
+    float i = INFINITY;
+    float x[] = {NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN};
+    BLASFUNC(sscal)(&N, &i, x, &incX);
+    ASSERT_TRUE(isnan(x[0]));
+    ASSERT_TRUE(isnan(x[8]));
+}
+
+CTEST(sscal, inf_nan_inc_2)
+{
+    blasint N=9;
+    blasint incX=2;
+    float i = INFINITY;
+    float x[] = {NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN,
+                    NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN};
+    BLASFUNC(sscal)(&N, &i, x, &incX);
+    ASSERT_TRUE(isnan(x[0]));
+    ASSERT_TRUE(isnan(x[8]));
+}
+
+#endif
+
+#ifdef BUILD_DOUBLE
+
+#ifndef NAN
+#define NAN 0.0/0.0
+#endif
+#ifndef INFINITY
+#define INFINITY 1.0/0.0
+#endif
+
+CTEST(dscal, 0_nan)
+{
+    blasint N=9;
+    blasint incX=1;
+    double i = 0.0;
+    double x[] = {NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN};
+    BLASFUNC(dscal)(&N, &i, x, &incX);
+    ASSERT_TRUE(isnan(x[0]));
+    ASSERT_TRUE(isnan(x[8]));
+}
+
+CTEST(dscal, 0_nan_inc_2)
+{
+    blasint N=9;
+    blasint incX=2;
+    double i = 0.0;
+    double x[] = {NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN,
+                    NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN};
+    BLASFUNC(dscal)(&N, &i, x, &incX);
+    ASSERT_TRUE(isnan(x[0]));
+    ASSERT_TRUE(isnan(x[8]));
+}
+
+CTEST(dscal, nan_0)
+{
+    blasint N=9;
+    blasint incX=1;
+    double i = NAN;
+    double x[] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
+    BLASFUNC(dscal)(&N, &i, x, &incX);
+    ASSERT_TRUE(isnan(x[0]));
+    ASSERT_TRUE(isnan(x[8]));
+}
+
+CTEST(dscal, nan_0_inc_2)
+{
+    blasint N=9;
+    blasint incX=2;
+    double i = NAN;
+    double x[] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+                    0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
+    BLASFUNC(dscal)(&N, &i, x, &incX);
+    ASSERT_TRUE(isnan(x[0]));
+    ASSERT_TRUE(isnan(x[8]));
+}
+
+CTEST(dscal, 0_inf)
+{
+    blasint N=9;
+    blasint incX=1;
+    double i = 0.0;
+    double x[] = {INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY};
+    BLASFUNC(dscal)(&N, &i, x, &incX);
+    ASSERT_TRUE(isnan(x[0]));
+    ASSERT_TRUE(isnan(x[8]));
+}
+
+CTEST(dscal, 0_inf_inc_2)
+{
+    blasint N=9;
+    blasint incX=2;
+    double i = 0.0;
+    double x[] = {INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY,
+                    INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY};
+    BLASFUNC(dscal)(&N, &i, x, &incX);
+    ASSERT_TRUE(isnan(x[0]));
+    ASSERT_TRUE(isnan(x[8]));
+}
+
+CTEST(dscal, inf_0)
+{
+    blasint N=9;
+    blasint incX=1;
+    double i = INFINITY;
+    double x[] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
+    BLASFUNC(dscal)(&N, &i, x, &incX);
+    ASSERT_TRUE(isnan(x[0]));
+    ASSERT_TRUE(isnan(x[8]));
+}
+
+CTEST(dscal, inf_0_inc_2)
+{
+    blasint N=9;
+    blasint incX=2;
+    double i = INFINITY;
+    double x[] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+                    0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
+    BLASFUNC(dscal)(&N, &i, x, &incX);
+    ASSERT_TRUE(isnan(x[0]));
+    ASSERT_TRUE(isnan(x[8]));
+}
+
+CTEST(dscal, nan_inf)
+{
+    blasint N=9;
+    blasint incX=1;
+    double i = NAN;
+    double x[] = {INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY};
+    BLASFUNC(dscal)(&N, &i, x, &incX);
+    ASSERT_TRUE(isnan(x[0]));
+    ASSERT_TRUE(isnan(x[8]));
+}
+
+CTEST(dscal, nan_inf_inc_2)
+{
+    blasint N=9;
+    blasint incX=2;
+    double i = NAN;
+    double x[] = {INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY,
+                    INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY};
+    BLASFUNC(dscal)(&N, &i, x, &incX);
+    ASSERT_TRUE(isnan(x[0]));
+    ASSERT_TRUE(isnan(x[8]));
+}
+
+CTEST(dscal, inf_nan)
+{
+    blasint N=9;
+    blasint incX=1;
+    double i = INFINITY;
+    double x[] = {NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN};
+    BLASFUNC(dscal)(&N, &i, x, &incX);
+    ASSERT_TRUE(isnan(x[0]));
+    ASSERT_TRUE(isnan(x[8]));
+}
+
+CTEST(dscal, inf_nan_inc_2)
+{
+    blasint N=9;
+    blasint incX=2;
+    double i = INFINITY;
+    double x[] = {NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN,
+                    NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN};
+    BLASFUNC(dscal)(&N, &i, x, &incX);
+    ASSERT_TRUE(isnan(x[0]));
+    ASSERT_TRUE(isnan(x[8]));
+}
+
+#endif
+
+#ifdef BUILD_COMPLEX
+
+#ifndef NAN
+#define NAN 0.0/0.0
+#endif
+#ifndef INFINITY
+#define INFINITY 1.0/0.0
+#endif
+
+CTEST(cscal, i_nan)
+{
+    blasint N=9;
+    blasint incX=1;
+    float i[] = {0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1 };
+    float nan[] = {NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0};
+    BLASFUNC(cscal)(&N, i, nan, &incX);
+    ASSERT_TRUE(isnan(nan[0]));
+    ASSERT_TRUE(isnan(nan[1]));
+    ASSERT_TRUE(isnan(nan[16]));
+    ASSERT_TRUE(isnan(nan[17]));
+}
+
+CTEST(cscal, i_nan_inc_2)
+{
+    blasint N=9;
+    blasint incX=2;
+    float i[] = {0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1 };
+    float nan[] = {NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0,
+                    NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0};
+    BLASFUNC(cscal)(&N, i, nan, &incX);
+    ASSERT_TRUE(isnan(nan[0]));
+    ASSERT_TRUE(isnan(nan[1]));
+    ASSERT_TRUE(isnan(nan[16]));
+    ASSERT_TRUE(isnan(nan[17]));
+}
+
+CTEST(cscal, nan_i)
+{
+    blasint N=9;
+    blasint incX=1;
+    float i[] = {0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1 };
+    float nan[] = {NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0};
+    BLASFUNC(cscal)(&N, nan, i, &incX);
+    ASSERT_TRUE(isnan(i[0]));
+    ASSERT_TRUE(isnan(i[1]));
+    ASSERT_TRUE(isnan(i[16]));
+    ASSERT_TRUE(isnan(i[17]));
+}
+
+CTEST(cscal, nan_i_inc_2)
+{
+    blasint N=9;
+    blasint incX=2;
+    float i[] = {0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1,
+                  0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1 };
+    float nan[] = {NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0};
+    BLASFUNC(cscal)(&N, nan, i, &incX);
+    ASSERT_TRUE(isnan(i[0]));
+    ASSERT_TRUE(isnan(i[1]));
+    ASSERT_TRUE(isnan(i[16]));
+    ASSERT_TRUE(isnan(i[17]));
+}
+
+CTEST(cscal, i_inf)
+{
+    blasint N=9;
+    blasint incX=1;
+    float i[] = {0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1 };
+    float inf[] = {INFINITY, 0, INFINITY,0, INFINITY,0, INFINITY,0, INFINITY,0, INFINITY,0, INFINITY,0, INFINITY,0, INFINITY,0};
+    BLASFUNC(cscal)(&N, i, inf, &incX);
+    ASSERT_TRUE(isnan(inf[0]));
+    ASSERT_TRUE(isinf(inf[1]));
+    ASSERT_TRUE(isnan(inf[16]));
+    ASSERT_TRUE(isinf(inf[17]));
+}
+
+CTEST(cscal, i_inf_inc_2)
+{
+    blasint N=9;
+    blasint incX=2;
+    float i[] = {0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1 };
+    float inf[] = {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,0};
+    BLASFUNC(cscal)(&N, i, inf, &incX);
+    ASSERT_TRUE(isnan(inf[0]));
+    ASSERT_TRUE(isinf(inf[1]));
+    ASSERT_TRUE(isnan(inf[16]));
+    ASSERT_TRUE(isinf(inf[17]));
+}
+
+CTEST(cscal, inf_i)
+{
+    blasint N=9;
+    blasint incX=1;
+    float i[] = {0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1 };
+    float inf[] = {INFINITY, 0, INFINITY,0, INFINITY,0, INFINITY,0, INFINITY,0, INFINITY,0, INFINITY,0, INFINITY,0, INFINITY,0};
+    BLASFUNC(cscal)(&N, inf, i, &incX);
+    ASSERT_TRUE(isnan(i[0]));
+    ASSERT_TRUE(isinf(i[1]));
+    ASSERT_TRUE(isnan(i[16]));
+    ASSERT_TRUE(isinf(i[17]));
+}
+
+CTEST(cscal, inf_i_inc_2)
+{
+    blasint N=9;
+    blasint incX=2;
+    float i[] = {0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1,
+                  0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1 };
+    float inf[] = {INFINITY, 0, INFINITY,0, INFINITY,0, INFINITY,0, INFINITY,0, INFINITY,0, INFINITY,0, INFINITY,0, INFINITY,0};
+    BLASFUNC(cscal)(&N, inf, i, &incX);
+    ASSERT_TRUE(isnan(i[0]));
+    ASSERT_TRUE(isinf(i[1]));
+    ASSERT_TRUE(isnan(i[16]));
+    ASSERT_TRUE(isinf(i[17]));
+}
+
+CTEST(cscal, i_0inf)
+{
+    blasint N=9;
+    blasint incX=1;
+    float i[] = {0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1 };
+    float inf[] = {0,INFINITY, 0, INFINITY,0, INFINITY,0, INFINITY,0, INFINITY,0, INFINITY,0, INFINITY,0, INFINITY,0, INFINITY};
+    BLASFUNC(cscal)(&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(cscal, i_0inf_inc_2)
+{
+    blasint N=9;
+    blasint incX=2;
+    float i[] = {0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1 };
+    float 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(cscal)(&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
+
 #ifdef BUILD_COMPLEX16
 
 #ifndef NAN
@@ -25,7 +469,7 @@ CTEST(zscal, i_nan)
 CTEST(zscal, i_nan_inc_2)
 {
     blasint N=9;
-    blasint incX=1;
+    blasint incX=2;
     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,
                     NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0};
@@ -52,7 +496,7 @@ CTEST(zscal, nan_i)
 CTEST(zscal, nan_i_inc_2)
 {
     blasint N=9;
-    blasint incX=1;
+    blasint incX=2;
     double i[] = {0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1,
                   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};