Проблема с суммой (массив, маска = ...)

Question

У меня проблема в моей программе. Внутренняя сумма функций, вызванная с помощью маски, приводит к подозрительным результатам: когда я выполняю среднее значение, я получаю значение из границ массива. Я подозреваю, что это связано с ошибками округления. Я работаю с большими массивами, и ошибки округления приводят к большим отклонениям (разница примерно на 40% по сравнению с ожидаемым значением в размере 40 000 элементов).

Ниже приведен минимальный пример для его воспроизведения и связанный с ним выход.

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 является:

test : 40019.9961  40011.9961 

теоретических результатов 40020.

Запуск GNU Fortran (GCC) 4.9.0

Answer 1

Вы работаете с одинарной точностью массивов. Компьютер в основном сохраняет реальные числа как расширение в степени 2. Это отлично работает для чисел, таких как 2 и 4 и 8 и т. Д., Которые легко могут быть записаны как целые степени двух с целыми коэффициентами, но менее хорошо для некоторых действительных чисел (например, 1.d0/2000.d0).

С одинарной точностью

real, allocatable, dimension(:) 

4 байта выделены. Это даст вам 8 цифр точности. И это то, что вы наблюдаете. Вторая сумма

sum(1.+real_array,mask=log_array) 

имеет только четыре цифры точности, но, хорошо, вы добавляете 1.0 и то, что в 1000 раз меньше. Это сужает его до четырех цифр (это то, что вы наблюдаете во втором случае).

Вы можете улучшить это, объявив все двойную точность (aka 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), которые снова можно было бы ожидать.