gfortran 4.8.5 处理的 SIGFPE 错误

Question

我正在使用在 Ubuntu 16.04 LTS 上使用 gfortran 4.8.5 版编译的计算流体动力学软件。该软件可以使用单精度或双精度和 -O3 优化选项进行编译。由于我没有运行双精度 CFD 软件所需的计算资源，因此我使用单精度和以下选项编译它

ffpe-trap=invalid,zero,overflow

我在包含 asin 函数的代码行上收到 SIGFPE 错误 -

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

换句话说，反sin函数和这段代码是一个以jp和ip为索引的双重嵌套FOR循环的一部分。目前，由于各种其他原因，软件人员无法帮助我，因此我正在尝试自行调试。 SIGFPE 错误仅在单精度编译而非双精度编译中观察到。

我在我的代码中的失败代码行（即 asin 函数调用）之前插入了以下打印语句。这会帮助我解决我面临的问题吗？这段代码针对每个时间步执行，并且在一系列时间步之后发生。或者我可以采取哪些其他步骤来帮助我解决此问题？将 "precision" 添加到编译器标志有帮助吗？

  if (zsm(ip,jp) >= 1.0 .or. zsm(ip,jp) <= -1.0) then
       print *,zsm(ip,jp),ip,jp
  end if

编辑我看了一下这个答案 Unexpected behavior of asin in R，我想知道我是否可以在 fortran 中做类似的事情，即使用 max 函数。如果它低于 -1 或大于 1，则以适当的方式对其进行四舍五入。我如何使用 max 函数使用 gfortran 来做到这一点？

在我的桌面上，以下程序执行没有问题（即它能够正确处理带符号的零），所以我猜测 SIGFPE 错误发生在参数大于 1 或小于 -1 时。

 program testa

 real a,x

 x = -0.0000
 a = asin(x)
 print *,a
 end program testa

Answer 1

我们有 min and max functions in Fortran, so I think we can use the same method as in the linked page，即 asin( max(-1.0,min(1.0,x) )。我用 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

给出 wp = sp（或在我的计算机上为 wp = kind(0.0)）

$ ./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

和 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

如果需要大量修改 asin(x) 并且程序依赖于 C 或 Fortran 预处理器，定义一些宏可能会很方便

#define clamp(x) max(-1.0_wp,min(1.0_wp,x))

并将其用作asin( clamp(x) )。如果我们想删除这样的修改，我们可以简单地将 clamp() 的定义更改为 #define clamp(x) (x)。另一种方法可能是定义一些 asin2(x) 函数，将 x 限制为 [-1,1] 并将内置的 asin 替换为 asin2 （作为宏或Fortran 函数）。

gfortran 4.8.5 处理的 SIGFPE 错误

SIGFPE error with gfortran 4.8.5 handling

fortran

gfortran

fortran90

floating-point-exceptions