发出sum(数组，mask=.)

Question

我的节目有个问题。带掩码调用的内禀函数和导致可疑结果:当我执行一个平均值时，我获得一个超出数组界限的值。我怀疑这与舍入错误有关。我正在处理大数组，舍入误差会导致很大的偏差(与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，或者添加一个编译器标志，比如-fdefaul-real-8-fdefaul-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时丢失的4位数字)，这也是您所期望的。