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

Question

I have an issue in my program. The intrinsic function sum called with a mask leads to suspicious results: when I perform an average, I obtain a value out of the array bounds. I suspect this is related to rounding errors. I am working with large arrays and the rounding errors lead to large deviations (difference around 40% compared with the expected value for a size of 40,000 elements).

Below is a minimal example to reproduce it, and the associated output.

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

Ouput is :

test :    40019.9961       40011.9961

Theoretical results is 40020.

Running GNU Fortran (GCC) 4.9.0

Possible duplicate of [Is floating point math broken?](http://stackoverflow.com/questions/588004/is-floating-point-math-broken) — Alexander Vogt, May 23 '16 at 17:28

score 1 · Accepted Answer · edited May 23 '17 at 12:31

You are working with single precision arrays. The computer basically stores real numbers as an expansion in the power of 2. This works fine for numbers like 2 and 4 and 8 and so on, which can easily be written as integer powers of two with integer coefficients but less well for some real numbers (like 1.d0/2000.d0).

With single precision

real, allocatable, dimension(:)

4 bytes are allocated. Which will give you 8 digits of precision. And this is what you observe. The second sum

sum(1.+real_array,mask=log_array)

has only four digits of precision, but, well, you are adding 1.0 and something that is 1000 times smaller. Which narrows it down further to effectively four digits (which is what you observe in the second case).

You can improve upon this by declaring everything double precision (aka 8 byte variables with 16 digits precision) and either instead of 1.0 you will have to write 1.d0, or you add a compiler flag like -fdefault-real-8 -fdefault-double-8.

If your rounding errors accumulate that much during your operations, I recommend rethinking the order of doing things. Adding variables of vastly different scopes will bring your precision down significantly.

If that is not an option and double precision is not enough, I can point you to quadruple precision

quad precision in gfortran

but I personally have not used it and since this usually solved through a software layer expect a huge performance loss.

edit: tried with double precision:

change :

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

keep the rest and compile with the mentioned compiler options. I obtain

 test :    40020.000000000000        40019.999999987354

The first result is fine, the second one is 12 digits precision (originally 16 digits plus four digits lost by adding 1.0 and 1.0/2000.0) which again is what you would expect.

Although I can not convert all real arrays to double precision, locally switching to double precision worked : real(sum(dble(myarray),mask=mymask)). Which is an ugly trick but a working one... Thx — user1824346, May 23 '16 at 16:01

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

1 Answers1