问题与sum（array，mask = ...）

我在我的程序中有一个问题。用掩码调用的内部函数总和会导致可疑的结果：当我执行平均值时，我从数组边界中获得一个值。我怀疑这与舍入错误有关。我正在处理大数组，并且舍入误差会导致较大的偏差（与40,000个元素大小的预期值相比，差异约为40％）。问题与sum（array，mask = ...）

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

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

来源

2016-05-23 user1824346

可能重复[浮点数学是否被破坏？]（http://stackoverflow.com/questions/588004/is-floating-point-math-broken） –

你用单精度数组工作。计算机基本上将实数存储为2的幂的展开式。这适用于像2和4和8等数字，对于某些实数可以很容易地用整数系数表示为2的整数幂，但对于某些实数（如1.d0/2000.d0）。

随着单精度

real, allocatable, dimension(:)

4个字节被分配。这会给你8位数的精度。这就是你观察到的。第二个总和

sum(1.+real_array,mask=log_array)

只有四位数字的精度，但是，好吧，您正在添加1.0和一些小1000倍的东西。将其缩小到有效四位数字（这是您在第二种情况下观察到的情况）。

你可以通过声明所有的双精度（也就是8个字节的变量，精度为16位）来改进它，而不是1.0，你将不得不编写1.d0，或者添加一个编译器标志，如-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丢失）这又是你所期望的。

来源

2016-05-23 15:37:32

虽然我不能将所有真正的数组转换为双精度，本地切换到双精度工作：真正（sum（dble（myarray），mask = mymask））。这是一个丑陋的技巧，但工作之一... Thx – user1824346

问题与sum（array，mask = ...）

回答

相关问题