实数相乘时如何避免下溢?
How to avoid underflow when multiplying real numbers?
我想知道在编译我的 Fortran 代码时是否有正确的方法将浮点数(或双精度数)相乘或平方而不出现下溢错误
gfortran -ffpe-trap=invalid,zero,overflow,underflow ...
我知道 underflow
选项并不总是一个好主意,但我想知道是否可以使用此选项进行乘法运算。事实上,在下面的示例中,我知道可能会发生下溢,但也许我不知道我的代码中的其他情况。这就是为什么我想尽可能保留此选项的原因。
这是一个示例,其中我为矩阵的每个 x、y 索引计算向量 u;组成这些向量的 2 个值介于 0 和 1 之间。然后我计算其范数的平方。
非常合乎逻辑,因为这个平方运算,我会有值下溢。因为,这些非常小的值对我来说可以视为零。有没有比使用 if
比较更好 underflow
的方法?
implicit none
double :: u(100,100,2), uSqr(100,100)
integer :: x,y
DO x= 1, 100
DO y = 1, 100
CALL Poisin( u(x,y,:), x, y )
ENDDO
ENDDO
uSqr = u(:,:,1)*u(:,:,1) + u(:,:,2) * u(:,:,2) ! where comes the underflow errors
如果只是为了避免下溢,也许你想使用 hypot
或者你真的需要斜边的平方吗?
hypot
的实施应避免 sqrt(x**2+y**2)
的 overflow/underflow 问题。
您的答案着眼于在某些情况下避免过度 下溢的具体方法。这使用 hypot
函数。这部分是一个答案:如果你想避免下溢,可能有一种方法可以重写算法来避免它。
对于更一般的情况(例如在这个问题中),需要对异常标志进行精细控制是不合适的。但是,编译器通常会提供异常处理例程的接口。
一种可移植的方法是使用 Fortran 2003 的 IEEE 工具。[如果使用 gfortran,您至少需要 5.0 版,但也有类似的特定于编译器的方法可用。]
Fortran 定义了 IEEE 异常和标志。标志可以是安静的或发信号的。您想要的是下溢不是有用诊断的部分,以便在该计算后不影响下溢标志状态。
该标志被称为 IEEE_UNDERFLOW
。我们可以通过子程序调用 IEEE_GET_FLAG(IEEE_UNDERFLOW, value)
和 IEEE_SET_FLAG(IEEE_UNDERFLOW, value)
来查询和设置它的状态。如果我们期待但不关心下溢,我们还希望确保异常不会停止。子程序 IEEE_SET_HALTING_MODE(IEEE_UNDERFLOW, value)
控制此模式。
所以,一个带注释的例子。
use, intrinsic :: ieee_arithmetic, only : IEEE_SELECTED_REAL_KIND
use, intrinsic :: ieee_exceptions
implicit none
! We want an IEEE kind, but this doesn't ensure support for underflow control
integer, parameter :: rk=IEEE_SELECTED_REAL_KIND(6, 70)
! State preservation/restoration
logical should_halt, was_flagged
real(rk) x
! Get the original halting mode and signal state
call ieee_get_halting_mode(IEEE_UNDERFLOW, should_halt)
call ieee_get_flag(IEEE_UNDERFLOW, was_flagged)
! Ensure we aren't going to halt on underflow
call ieee_set_halting_mode(IEEE_UNDERFLOW, .FALSE.)
! The irrelevant computation
x=TINY(x)
x=x**2
! And restore our old state
call ieee_set_halting_mode(IEEE_UNDERFLOW, should_halt)
call ieee_set_flag(IEEE_UNDERFLOW, was_flagged)
end program
我想知道在编译我的 Fortran 代码时是否有正确的方法将浮点数(或双精度数)相乘或平方而不出现下溢错误
gfortran -ffpe-trap=invalid,zero,overflow,underflow ...
我知道 underflow
选项并不总是一个好主意,但我想知道是否可以使用此选项进行乘法运算。事实上,在下面的示例中,我知道可能会发生下溢,但也许我不知道我的代码中的其他情况。这就是为什么我想尽可能保留此选项的原因。
这是一个示例,其中我为矩阵的每个 x、y 索引计算向量 u;组成这些向量的 2 个值介于 0 和 1 之间。然后我计算其范数的平方。
非常合乎逻辑,因为这个平方运算,我会有值下溢。因为,这些非常小的值对我来说可以视为零。有没有比使用 if
比较更好 underflow
的方法?
implicit none
double :: u(100,100,2), uSqr(100,100)
integer :: x,y
DO x= 1, 100
DO y = 1, 100
CALL Poisin( u(x,y,:), x, y )
ENDDO
ENDDO
uSqr = u(:,:,1)*u(:,:,1) + u(:,:,2) * u(:,:,2) ! where comes the underflow errors
如果只是为了避免下溢,也许你想使用 hypot
或者你真的需要斜边的平方吗?
hypot
的实施应避免 sqrt(x**2+y**2)
的 overflow/underflow 问题。
您的答案着眼于在某些情况下避免过度 下溢的具体方法。这使用 hypot
函数。这部分是一个答案:如果你想避免下溢,可能有一种方法可以重写算法来避免它。
对于更一般的情况(例如在这个问题中),需要对异常标志进行精细控制是不合适的。但是,编译器通常会提供异常处理例程的接口。
一种可移植的方法是使用 Fortran 2003 的 IEEE 工具。[如果使用 gfortran,您至少需要 5.0 版,但也有类似的特定于编译器的方法可用。]
Fortran 定义了 IEEE 异常和标志。标志可以是安静的或发信号的。您想要的是下溢不是有用诊断的部分,以便在该计算后不影响下溢标志状态。
该标志被称为 IEEE_UNDERFLOW
。我们可以通过子程序调用 IEEE_GET_FLAG(IEEE_UNDERFLOW, value)
和 IEEE_SET_FLAG(IEEE_UNDERFLOW, value)
来查询和设置它的状态。如果我们期待但不关心下溢,我们还希望确保异常不会停止。子程序 IEEE_SET_HALTING_MODE(IEEE_UNDERFLOW, value)
控制此模式。
所以,一个带注释的例子。
use, intrinsic :: ieee_arithmetic, only : IEEE_SELECTED_REAL_KIND
use, intrinsic :: ieee_exceptions
implicit none
! We want an IEEE kind, but this doesn't ensure support for underflow control
integer, parameter :: rk=IEEE_SELECTED_REAL_KIND(6, 70)
! State preservation/restoration
logical should_halt, was_flagged
real(rk) x
! Get the original halting mode and signal state
call ieee_get_halting_mode(IEEE_UNDERFLOW, should_halt)
call ieee_get_flag(IEEE_UNDERFLOW, was_flagged)
! Ensure we aren't going to halt on underflow
call ieee_set_halting_mode(IEEE_UNDERFLOW, .FALSE.)
! The irrelevant computation
x=TINY(x)
x=x**2
! And restore our old state
call ieee_set_halting_mode(IEEE_UNDERFLOW, should_halt)
call ieee_set_flag(IEEE_UNDERFLOW, was_flagged)
end program