sum(array,mask=...) 问题

Question

我的程序有问题。使用掩码调用的内在函数 sum 会导致可疑结果：当我执行平均时，我获得了数组边界之外的值。我怀疑这与舍入误差有关。我正在处理大型数组，舍入误差导致较大偏差（与 40,000 个元素大小的预期值相比，差异约为 40%）。

下面是重现它的最小示例，以及相关的输出。

program main

  implicit none

  integer :: nelem
  real, allocatable, dimension(:) :: real_array
  logical, allocatable, dimension(:) :: log_array

  ! init
  nelem=40000
  allocate(real_array(nelem))
  allocate(log_array(nelem))
  real_array=0.
  log_array=.true.

  ! Fill arrays
  real_array=1./2000.
  log_array = real_array.le.(0.5)

  ! Test
  print *, ' test : ', &
           count(log_array)+sum(real_array, mask=log_array), &
           sum(1.+real_array,mask=log_array)

end program main

输出为：

test :    40019.9961       40011.9961

理论结果为40020。

运行 GNU Fortran (GCC) 4.9.0

Answer 1

您正在使用单精度数组。计算机基本上将实数存储为 2 的幂的扩展。这适用于 2、4 和 8 等数字，它们可以很容易地写成具有整数系数的 2 的整数幂，但对于某些实数则不太好（比如 1.d0/2000.d0）。

单精度

real, allocatable, dimension(:)

分配了 4 个字节。这将为您提供 8 位数的精度。这就是你观察到的。第二个总和

sum(1.+real_array,mask=log_array)

只有四位数的精度，但是，好吧，你要加上 1.0 和小 1000 倍的东西。这进一步缩小到有效的四位数（这是您在第二种情况下观察到的）。

您可以通过声明所有双精度（也就是具有 16 位精度的 8 字节变量）来改进这一点，并且您必须编写 1.d0 而不是 1.0，或者添加一个编译器标志，如 -fdefault-real -8 -fdefault-double-8.

如果你的舍入误差在你的操作过程中积累了那么多，我建议你重新考虑做事的顺序。添加范围大不相同的变量会显着降低精度。

如果这不是一个选项并且双精度不够，我可以向您指出四倍精度

quad precision in gfortran

但我个人没有使用过，因为这通常通过软件解决层预计会有巨大的性能损失。

编辑：尝试双精度：

变化：

double precision, allocatable, dimension(:) :: real_array

保留其余部分并使用提到的编译器选项进行编译。我获得

 test :    40020.000000000000        40019.999999987354

第一个结果很好，第二个是 12 位精度（最初是 16 位加上 1.0 和 1.0/2000.0 相加后丢失的四位），这也是您所期望的。

sum(array,mask=...) 问题

Issue with sum(array,mask=...)

arrays

fortran

sum

gfortran