Sun Studio 12:Fortran 编程指南

6.5 更深层次的数值风险

本部分解决一些运算操作无意间可能会产生无效、被零除、上溢、下溢或不精确异常的实际问题。

例如,在 IEEE 标准之前,如果在计算机中将两个非常小的数相乘,结果可能为零。多数大型机和小型机的行为亦是如此。使用 IEEE 运算,渐进下溢会扩大动态计算范围。

例如,假设某一 32 位处理器以 1.0E-38 作为机器中可表示的最小值 epsilon。将两个小数相乘:


      a = 1.0E-30
      b = 1.0E-15
      x = a * b

较早的运算会得到 0.0,但如果使用 IEEE 运算和相同的字长,却会得到 1.40130E-45。此时便出现了下溢,告诉您答案比机器自然表示的值小。该结果是通过从尾数中“窃取”一些位并将其移交给指数来完成的。得到的结果(即一个非规范化数)从某种意义上说精确度比较低,但从另外一种意义上说精确度又比较高。其深层含意已超出了本次讨论的范围。如果有兴趣,可以参考《Computer》1980 年 1 月,第 13 卷,第 1 期,尤其是 J. Coonen 的文章 "Underflow and the Denormalized Numbers"。

大多数科学程序都有对舍入很敏感的代码段,通常是在方程求解或矩阵因子分解中。若不采用渐进下溢,程序员就得自己实现检测接近不准确阈值的方法。否则,他们就必须放弃对实现强大、稳定算法的追求。

有关这些主题的详细信息,请参见《数值计算指南》。

6.5.1 避免简单下溢

有些应用程序实际上执行了许多非常接近零的计算。这在计算残数或微分修正的算法中很常见。为获得在数值上安全的最大性能,需要采用扩展精度运算来执行关键计算。如果应用程序是单精度应用程序,可采用双精度执行关键计算。

示例:采用单精度的简单点积计算:


      sum = 0
      DO i = 1, n
         sum = sum + a(i) * b(i)
      END DO

如果 a(i)b(i) 非常小,会出现很多下溢。通过强制计算采用双精度,计算点积时会具有更高的准确性,并且可避免出现下溢情况:


      DOUBLE PRECISION sum
      DO i = 1, n
         sum = sum + dble(a(i)) * dble(b(i))
      END DO
      result = sum

通过增加对库例程 nonstandard_arithmetic() 的调用,或通过使用 -fns 选项编译应用程序的主程序,可以强制 SPARC 处理器在涉及到下溢时(存储零)像较早的系统一样进行处理。

6.5.2 以错误答案继续

您可能会对在答案明显错误的情况下为什么能继续进行计算感到奇怪。IEEE 运算允许您区分可以忽略的错误答案的种类,如 NaNInf。然后,可以根据此种区分来作决定。

不妨考虑一个电路模拟的例子。在 50 行的特定计算中,(出于参数原因)唯一关注的变量是电压。进一步假设其值只可能是 +5v、0、-5v。

仔细安排计算的每一部分以强制每个子结果在正确范围内,这是很可能实现的:

此外,由于 Inf 不是允许值,所以需要特殊的逻辑来确保不会与较大的数相乘。

利用 IEEE 运算,此逻辑可以简化许多。计算可以用显而易见的方式编写,并且只需强制最终结果为正确的值-因为 Inf 可以出现并且可以很容易地测出。

此外,还可检测到 0/0 的特殊情况并按照您的意愿进行处理。结果更易读取且执行起来更快,因为无需进行不必要的比较。

6.5.3 过度下溢

如果将两个非常小的数字相乘,结果将会出现下溢。

如果预知乘法(或减法)中的操作数可能会很小并且很可能会下溢,可采用双精度进行计算,而后再将结果转换成单精度。

例如,类似下面的点积循环:


  real sum, a(maxn), b(maxn)
  ...
  do i =1, n
      sum = sum + a(i)*b(i)
  enddo

其中,已知 a(*)b(*) 具有小元素,为保持数值准确度,应采用双精度来运行:


  real a(maxn), b(maxn)
  double sum
  ...
  do i =1, n
      sum = sum + a(i)*dble(b(i))
  enddo

鉴于以软件方式解决了原始循环造成的过度下溢,这样做还有可能提高性能。但是,对此并无绝对而快速的法则;只能通过对计算代码进行高强度实验来确定最有利的解决方案。