I am using a computational fluid dynamics software that is compiled with gfortran version 4.8.5 on Ubuntu 16.04 LTS. The software can be compiled with either single precision or double precision and the -O3 optimization option. As I do not have the necessary computational resources to run the CFD software on double precision I am compiling it with single precision and the following options
ffpe-trap=invalid,zero,overflow
I am getting a SIGFPE error on a line of code that contains the asin function-
INTEGER, PARAMETER :: sp = SELECTED_REAL_KIND( 6, 37) !< single precision
INTEGER, PARAMETER :: wp = sp
REAL(KIND=wp) zsm(:,:)
ela(i,j) = ASIN(zsm(ip,jp))
In other words the inverse sin function and this code is part of a doubly nested FOR loop with jp and ip as the indices. Currently the software staff is unable to help me for various other reasons and so I am trying to debug this on my own. The SIGFPE error is only being observed in the single precision compilation not double precision compilation.
I have inserted the following print statements in my code prior to the line of code that is failing i.e. the asin function call. Would this help me with unraveling the problem that I am facing ? This piece of code is executed for every time step and it is occurring after a series of time steps. Alternatively what other steps can I do to help me fix this problem ? Would adding "precision" to the compiler flag help ?
if (zsm(ip,jp) >= 1.0 .or. zsm(ip,jp) <= -1.0) then
print *,zsm(ip,jp),ip,jp
end if
EDIT I took a look at this answer Unexpected behavior of asin in R and I am wondering whether I could do something similar in fortran i.e. by using the max function. If it goes below -1 or greater than 1 then round it off in the proper manner. How can I do it with gfortran using the max function ?
On my desktop the following program executes with no problems(i.e. it has the ability to handle signed zeros properly) and so I am guessing the SIGFPE error occurs with either the argument greater than 1 or less than -1.
program testa
real a,x
x = -0.0000
a = asin(x)
print *,a
end program testa
We have min and max functions in Fortran, so I think we can use the same method as in the linked page, i.e., asin( max(-1.0,min(1.0,x) ). I have tried the following test with gfortran-4.8 & 7.1:
program main
implicit none
integer, parameter :: sp = selected_real_kind( 6, 37 )
integer, parameter :: wp = sp
! integer, parameter :: wp = kind( 0.0 )
! integer, parameter :: wp = kind( 0.0d0 )
real(wp) :: x, a
print *, "Input x"
read(*,*) x
print *, "x =", x
print *, "equal to 1 ? :", x == 1.0_wp
print *, asin( x )
print *, asin( max( -1.0_wp, min( 1.0_wp, x ) ) )
end
which gives with wp = sp (or wp = kind(0.0) on my computer)
$ ./a.out
Input x
1.00000001
x = 1.00000000
equal to 1 ? : T
1.57079625 (<- 1.5707964 for gfortran-4.8)
1.57079625
$ ./a.out
Input x
1.0000001
x = 1.00000012
equal to 1 ? : F
NaN
1.57079625
and with wp = kind(0.0d0)
$ ./a.out
Input x
1.0000000000000001
x = 1.0000000000000000
equal to 1 ? : T
1.5707963267948966
1.5707963267948966
$ ./a.out
Input x
1.000000000000001
x = 1.0000000000000011
equal to 1 ? : F
NaN
1.5707963267948966
If it is necessary to modify a lot of asin(x) and the program relies on a C or Fortran preprocessor, it may be convenient to define some macro like
#define clamp(x) max(-1.0_wp,min(1.0_wp,x))
and use it as asin( clamp(x) ). If we want to remove such a modification, we can simply change the definition of clamp() as #define clamp(x) (x). Another approach may be to define some asin2(x) function that limits x to [-1,1] and replace the built-in asin by asin2 (either as a macro or a Fortran function).