2

我正在使用在 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))

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

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

  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，则以适当的方式四舍五入。如何使用 gfortran 使用 max 函数来做到这一点？

在我的桌面上，以下程序可以毫无问题地执行（即它能够正确处理带符号的零），因此我猜测 SIGFPE 错误发生在参数大于 1 或小于 -1 的情况下。

 program testa

 real a,x

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

1 回答 1

2

我们在 Fortran 中有min和max函数，所以我认为我们可以使用与链接页面中相同的方法，即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 函数）。

于 2017-06-24T16:03:20.420 回答