floraachy
/
OpenBLAS
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

handle corner cases involving NAN and/or INF

This commit is contained in:
Martin Kroeker 2024-06-06 23:59:43 +02:00 committed by GitHub
parent ffc1ab3f6e
commit 1abafcd9b2
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
3 changed files with 76 additions and 32 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

59
kernel/x86_64/cscal.c
View File

@ -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++;

5
kernel/x86_64/dscal.c
View File

@ -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
{

44
kernel/x86_64/sscal.c
View File

@ -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);