实数相乘时如何避免下溢？答案

【问题标题】：How to avoid underflow when multiplying real numbers?实数相乘时如何避免下溢？
【发布时间】：2017-04-11 17:37:44
【问题描述】：

我想知道在编译我的 Fortran 代码时是否有适当的方法将浮点数（或双精度数）相乘或平方而不出现下溢错误

gfortran -ffpe-trap=invalid,zero,overflow,underflow ...

我知道underflow 选项并不总是一个好主意，但我想知道是否可以使用此选项进行乘法运算。事实上，在下面的示例中，我知道可能会发生下溢，但也许我不知道我的代码中有其他情况。这就是为什么我想尽可能保留这个选项。

这是一个例子，我为矩阵的每个 x,y 索引计算一个向量 u；组成这些向量的 2 个值介于 0 和 1 之间。然后我计算其范数的平方。

非常合乎逻辑，由于这个平方运算，我的值会下溢。因为，这些非常小的值对我来说可以被认为是零。有没有办法不使用underflow 而不是使用if 比较？

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

【问题讨论】：

明确地说，您是要避免下溢，还是要避免下溢成为错误条件？
你不能在每个文件的基础上设置-ffpe-trap 选项吗？如果是这样，也许您可以将上面引用的部分分离到一个源文件中，并从-ffpe-trap 列表中省略underflow 以编译此文件，同时保留其他源文件的原始异常陷阱列表。我假设您希望通过将 uSqr 值视为零来允许程序继续运行。
我必须支持@francescalus 的评论。你想保持计算而不是退出错误吗？为什么要或需要捕获underflow？认真考虑不要诱捕underflow。在我看来，这是非常无用的，甚至是有害的。
@norio 这就是我所知道的。但是有了这个解决方案，我不知道他们在这个文件中是否还有其他underflow 错误。
@VladimirF 和@francescalus 我想捕获underflow 以确保它们在我的其余代码中没有underflow 错误。所以我想知道他们是否可以更好地编写避免underflow 问题的示例。也许我的问题很愚蠢，答案是因为我进行平方运算underflow 和overflow 可能会发生，所以没有其他方法，除非我使用条件测试来阈值“小”和“大”值。

标签： floating-point fortran underflow

【解决方案1】：

您有一个答案，该答案着眼于在某些情况下避免过度下溢的特定方法。这使用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

【讨论】：

您的回答似乎很合适。这个周末我会在家里试试，我会验证你的答案。

【解决方案2】：

如果只是为了避免下溢，也许你想使用hypot 或者你真的需要斜边的平方吗？

hypot的实现应该避免sqrt(x**2+y**2)的上溢/下溢问题。

【讨论】：

你说得对，我一直认为它只是为了溢出问题，但它应该以同样的方式避免下溢。
感谢hypot 的提示，但我真的需要斜边的平方。
@Romain 所以计算hypot的平方，问题出在哪里......如果在此期间发生下溢，这是不可避免的。