首页 文章

使用gfortran 4.8.5处理时出现SIGFPE错误

提问于
浏览 次
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))

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

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

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

EDIT 我看了一下这个答案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
fortran gfortran fortran90 floating-point-exceptions

1 回答

  • 2

    我们在Fortran中有min和max函数，所以我认为我们可以使用与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]并用 asin2 替换内置的 asin （作为宏或Fortran函数） .

    回复于

相关问题