From f1248b849d11d6a436a6e19983c9c6695c143bc7 Mon Sep 17 00:00:00 2001
From: Martin Kroeker <martin@ruby.chemie.uni-freiburg.de>
Date: Sat, 22 Jun 2024 15:55:29 +0200
Subject: [PATCH] handle INF and NAN in input

---
 kernel/power/dscal.c | 57 ++++++++++++++++++++++++++++++++++----------
 1 file changed, 44 insertions(+), 13 deletions(-)

diff --git a/kernel/power/dscal.c b/kernel/power/dscal.c
index 96c4e51bc..2bbc1ea6d 100644
--- a/kernel/power/dscal.c
+++ b/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++;
 			}