Функция sum возвращает ответы, отличные от явного цикла

Question

Я преобразовываю код f77 в код f90, и часть кода должна суммироваться по элементам трехмерной матрицы. В f77 это было достигнуто с помощью 3 петель (по внешнему, среднему, внутреннему индексам). Я решил использовать внутреннюю сумму f90 (3 раза) для достижения этой цели, и, к моему удивлению, ответы отличаются. Я использую компилятор ifort, у меня есть отладка, контрольные границы, все оптимизации не включены

Вот код в стиле f77

r1 = 0.0
do k=1,nz
  do j=1,ny
    do i=1,nx
      r1 = r1 + foo(i,j,k)
    end do
  end do
end do

, а вот код f90

r = SUM(SUM(SUM(foo, DIM=3), DIM=2), DIM=1)

Я испробовал все виды вариаций, такие как изменение порядка циклов для кода f77 или создание временных 2D-матриц и 1D-массивов, чтобы «уменьшить» размеры при использовании SUM, но явный стиль f77Циклы всегда дают разные ответы от функции f90 + SUM.

Буду признателен за любые предложения, которые помогут понять расхождение.

Кстати, при этом используется один последовательный процессор.

Отредактировано 12:13 вечера, чтобы показать полный пример

! ifort -check bounds -extend-source 132 -g -traceback -debug inline-debug-info -mkl -o verify  verify.f90
! ./verify

program verify

implicit none

integer :: nx,ny,nz

parameter(nx=131,ny=131,nz=131)

integer :: i,j,k
real :: foo(nx,ny,nz)
real :: r0,r1,r2
real :: s0,s1,s2
real :: r2Dfooxy(nx,ny),r1Dfoox(nx)

call random_seed
call random_number(foo)

r0 = 0.0
do k=1,nz
  do j=1,ny
    do i=1,nx
      r0 = r0 + foo(i,j,k)
    end do
  end do
end do

r1 = 0.0
do i=1,nx
  do j=1,ny
    do k=1,nz
      r1 = r1 + foo(i,j,k)
    end do
  end do
end do

r2 = 0.0
do j=1,ny
  do i=1,nx
    do k=1,nz
      r2 = r2 + foo(i,j,k)
    end do
  end do
end do

!*************************

s0 = 0.0
s0 = SUM(SUM(SUM(foo, DIM=3), DIM=2), DIM=1)

s1 = 0.0
r2Dfooxy = SUM(foo,   DIM = 3)
r1Dfoox  = SUM(r2Dfooxy, DIM = 2)
s1 = SUM(r1Dfoox)

s2 = SUM(foo)

!*************************

print *,'nx,ny,nz = ',nx,ny,nz
print *,'size(foo) = ',size(foo)

write(*,'(A,4(ES15.8))') 'r0,r1,r2          = ',r0,r1,r2
write(*,'(A,3(ES15.8))') 'r0-r1,r0-r2,r1-r2 = ',r0-r1,r0-r2,r1-r2

write(*,'(A,4(ES15.8))') 's0,s1,s2          = ',s0,s1,s2
write(*,'(A,3(ES15.8))') 's0-s1,s0-s2,s1-s2 = ',s0-s1,s0-s2,s1-s2

write(*,'(A,3(ES15.8))') 'r0-s1,r1-s1,r2-s1    = ',r0-s1,r1-s1,r2-s1

stop
end

!**********************************************

sample output

nx,ny,nz =          131         131         131
size(foo) =      2248091

r0,r1,r2          =  1.12398225E+06 1.12399525E+06 1.12397238E+06
r0-r1,r0-r2,r1-r2 = -1.30000000E+01 9.87500000E+00 2.28750000E+01
s0,s1,s2          =  1.12397975E+06 1.12397975E+06 1.12398225E+06
s0-s1,s0-s2,s1-s2 =  0.00000000E+00-2.50000000E+00-2.50000000E+00
r0-s1,r1-s1,r2-s1    =  2.50000000E+00 1.55000000E+01-7.37500000E+00

Answer 1

Теперь я вижу разницу. Это типичные ошибки округления при добавлении небольших чисел к большой сумме. Процессору разрешено использовать любой порядок суммирования, который он хочет. нет "правильного" заказа . Вы можете не действительно сказать, что исходные циклы дают "правильный" ответ, а другие нет.

Что вы можете сделать, это использовать двойную точность . В экстремальных обстоятельствах есть уловки, такие как суммирование по Кахану, но это редко требуется.

Добавление небольшого числа к большой сумме является неточным, особенно с одинарной точностью. В вашем результате все еще есть четыре значащие цифры.

---

Обычно один из них не использует аргумент DIM=, который используется в определенных особых обстоятельствах.

Если вы хотите сложить всеэлементы foo, используйте только

s0 = SUM(foo)

Этого достаточно.

Что

s0 = SUM(SUM(SUM(foo, DIM=3), DIM=2), DIM=1)

делает то, что он создает временные двумерные массивы с каждым элементомбыть суммой соответствующей строки в измерении z, затем одномерный массив с каждым элементом суммой по последнему измерению двумерного массива и, наконец, суммой этого одномерного массива. Если все сделано хорошо, конечный результат будет таким же, но он потребляет много циклов ЦП.

Answer 2

Встроенная функция sum возвращает зависящее от процессора приближение к сумме элементов аргумента массива. Это не то же самое, что последовательное добавление всех элементов.

Просто найти массив x, где

summation = x(1) + x(2) + x(3)

(выполняется строго слева направо) - не лучшее приближение для суммы, рассматривающей значения как "математические вещественные числа", ачем числа с плавающей точкой.

---

В качестве конкретного примера, чтобы взглянуть на природу аппроксимации с помощью ifort, мы можем взглянуть на следующую программу. Нам нужно включить оптимизацию, чтобы увидеть эффекты;важность порядка суммирования очевидна даже при отключенных оптимизациях (с -O0 или -debug).

  implicit none

  integer i
  real x(50)
  real total

  x = [1.,(EPSILON(0.)/2, i=1, SIZE(x)-1)]
  total = 0
  do i=1, SIZE(x)
     total = total+x(i)
     print '(4F17.14)', total, SUM(x(:i)), SUM(DBLE(x(:i))), REAL(SUM(DBLE(x(:i))))
  end do
end program

Если сложить в строгом порядке, мы получим 1.,видя, что все, что меньше величины epsilon(0.), не влияет на сумму.

Вы можете поэкспериментировать с размером массива и порядком его элементов, масштабированием небольших чисел и компиляцией с плавающей запятой ifort. параметры (такие как -fp-model strict, -mieee-fp, -pc32). Вы также можете попытаться найти пример, подобный приведенному выше, используя двойную точность вместо действительного по умолчанию.

Answer 3

Во-первых, добро пожаловать в StackOverflow. Пожалуйста, возьмите тур! Есть причина, по которой мы ожидаем Минимальный, Полный и Проверяемый пример , потому что мы смотрим на ваш код и можем только догадываться, в чем может быть дело, и этоне слишком полезен для сообщества.

Я надеюсь, что следующие предложения помогут вам выяснить, что происходит.

Используйте функцию size () и выведите то, что, по мнению Фортрана, является размерами измерений. а также печать nx, ny и nz. Насколько нам известно, массив объявлен больше, чем nx, ny и nz, и эти переменные установлены в соответствии с набором данных. Fortran не обязательно инициализирует массивы нулями в зависимости от того, является ли это статическим или размещаемым массивом.

Вы также можете попробовать указать экстенты массива в функции sum:

r = Sum(foo(1:nx,1:ny,1:nz))

Если все сделано так, по крайней мере, мы знаем, что функция sum работает с тем же фрагментом foo, что и цикл loop.

Если это так, вы получите неправильный ответ, даже если в этом нет ничего «неправильного»с кодом. Вот почему особенно важно привести этот пример Минимальный, Полный и Проверяемый.