From 1abafcd9b2c32a07a67f9db7b57485c464828bad Mon Sep 17 00:00:00 2001
From: Martin Kroeker <martin@ruby.chemie.uni-freiburg.de>
Date: Thu, 6 Jun 2024 23:59:43 +0200
Subject: [PATCH] handle corner cases involving NAN and/or INF

---
 kernel/x86_64/cscal.c | 59 ++++++++++++++++++++++++++++++-------------
 kernel/x86_64/dscal.c |  5 ++++
 kernel/x86_64/sscal.c | 44 +++++++++++++++++++++-----------
 3 files changed, 76 insertions(+), 32 deletions(-)

diff --git a/kernel/x86_64/cscal.c b/kernel/x86_64/cscal.c
index 95a99b8b9..212a21594 100644
--- a/kernel/x86_64/cscal.c
+++ b/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++;
 
diff --git a/kernel/x86_64/dscal.c b/kernel/x86_64/dscal.c
index 0b028b782..e7182c5ce 100644
--- a/kernel/x86_64/dscal.c
+++ b/kernel/x86_64/dscal.c
@@ -234,6 +234,11 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da, FLOAT *x, BLAS
 			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
 	{
 
diff --git a/kernel/x86_64/sscal.c b/kernel/x86_64/sscal.c
index af1220f1b..a85d20564 100644
--- a/kernel/x86_64/sscal.c
+++ b/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);