https://gcc.gnu.org/bugzilla/show_bug.cgi?id=91337
Bug ID: 91337
Summary: gfortran skips an if statement with some mathematical
optimisations with complex numbers.
Product: gcc
Version: 9.1.0
Status: UNCONFIRMED
Severity: normal
Priority: P3
Component: fortran
Assignee: unassigned at gcc dot gnu.org
Reporter: chinoune.mehdi at hotmail dot com
Target Milestone: ---
I have encountered some underflows/overflows in my code compiled with -Ofast,
and after investigations it seems like the complex abs gives zero with small
numbers. So I added a workaround. but it didn't work:
module m
implicit none
!
integer, parameter :: sp = selected_real_kind(6)
real(sp), parameter :: tiny_sp = tiny(1._sp), sqrt_tiny = sqrt( tiny_sp )
!
contains
subroutine sub(z,y)
complex(sp), intent(in) :: z
real(sp), intent(out) :: y
real(sp) :: az
!
az = abs(z)
if( az<tiny_sp ) az = abs(z/sqrt_tiny)*sqrt_tiny
! if( az<tiny_sp .or. az==0._sp ) az = abs(z/sqrt_tiny)*sqrt_tiny
print*, az<tiny_sp, az
y = 1._sp/az
!
end subroutine sub
end module m
!
program bug_skip_if
use m
implicit none
complex(sp) :: z
real(sp) :: y
!
z = (1.e-19_sp,0._sp)
call sub(z,y)
print*,"y = ",y
!
end program bug_skip_if
gfortran-9 -O1 -funsafe-math-optimizations -ffinite-math-only bug_skip_if.f90
-o test.x
./test.x
T 0.00000000
y = Infinity
I tried to write an equivalent c code, but it gives the right result:
#include <stdio.h>
#include <float.h>
#include <math.h>
#include <complex.h>
void foo(float complex z, float *y){
float const sqrt_flt_min = sqrtf( FLT_MIN );
float az;
az = cabsf(z);
if( az<FLT_MIN ) az = cabsf(z/sqrt_flt_min) * sqrt_flt_min;
printf("az = %.8e\n",az);
*y = 1.f/az;
}
int main(){
float complex z;
float y;
z = 1.e-19f + 0.*I;
foo(z,&y);
printf("y = %.7e\n",y);
return 0;
}
gcc-9 -O1 -funsafe-math-optimizations -ffinite-math-only cbug_skip_if.c -lm -o
test.x
./test.x
az = 9.99999968e-20
y = 1.0000000e+19
Q : Why does gfortran skip the if statement?
OS : Linux Ubuntu 18.04 x86_64
Compilers : gfortran 9.1.0 and 8.3.0
