Как numpy может быть намного быстрее моей процедуры на Фортране?

Я получаю массив 512 ^ 3, представляющий распределение температуры из моделирования (написанного на Фортране). Массив хранится в двоичном файле размером около 1/2 ГБ. Мне нужно знать минимум, максимум и среднее значение этого массива, и, поскольку мне скоро все равно понадобится разбираться в коде Фортрана, я решил попробовать и придумал следующую очень простую процедуру.

  integer gridsize,unit,j
  real mini,maxi
  double precision mean

  gridsize=512
  unit=40
  open(unit=unit,file='T.out',status='old',access='stream',&
       form='unformatted',action='read')
  read(unit=unit) tmp
  mini=tmp
  maxi=tmp
  mean=tmp
  do j=2,gridsize**3
      read(unit=unit) tmp
      if(tmp>maxi)then
          maxi=tmp
      elseif(tmp<mini)then
          mini=tmp
      end if
      mean=mean+tmp
  end do
  mean=mean/gridsize**3
  close(unit=unit)

Это занимает около 25 секунд на файл на машине, которую я использую. Это показалось мне довольно длинным, поэтому я сделал следующее на Python:

    import numpy

    mmap=numpy.memmap('T.out',dtype='float32',mode='r',offset=4,\
                                  shape=(512,512,512),order='F')
    mini=numpy.amin(mmap)
    maxi=numpy.amax(mmap)
    mean=numpy.mean(mmap)

Конечно, я ожидал, что это будет быстрее, но я действительно был потрясен. При идентичных условиях это занимает меньше секунды. Среднее значение отклоняется от того, которое находит моя подпрограмма Fortran (которую я также использовал со 128-битными числами с плавающей запятой, поэтому я как-то доверяю ему больше), но только на 7-й значащей цифре или около того.

Как Numpy может быть таким быстрым? Я имею в виду, что вам нужно просматривать каждую запись массива, чтобы найти эти значения, верно? Я делаю что-то очень глупое в моей программе Fortran, чтобы это заняло так много времени?

РЕДАКТИРОВАТЬ:

Чтобы ответить на вопросы в комментариях:

• Да, я также запускал процедуру Fortran с 32-битными и 64-битными числами с плавающей запятой, но это не повлияло на производительность.
• Я использовал iso_fortran_env128-битные числа с плавающей запятой.
• Используя 32-битные числа с плавающей запятой, мое среднее значение немного отличается, поэтому точность действительно является проблемой.
• Я запускал обе процедуры для разных файлов в разном порядке, так что кеширование должно было быть справедливым при сравнении, как я полагаю?
• Я действительно пробовал открыть MP, но читать из файла одновременно с разных позиций. После прочтения ваших комментариев и ответов это звучит действительно глупо, и это сделало рутину намного дольше. Я мог бы попробовать с ним операции с массивом, но, может быть, это даже не понадобится.
• Файлы на самом деле имеют размер 1/2 ГБ, это была опечатка, спасибо.
• Сейчас попробую реализовать массив.

РЕДАКТИРОВАТЬ 2:

Я реализовал то, что @Alexander Vogt и @casey предложили в своих ответах, и это так же быстро, как, numpyно теперь у меня проблема с точностью, как указал @Luaan, я мог бы получить. При использовании 32-битного массива с плавающей запятой среднее значение, вычисленное по нему, sumсоставляет 20%. Делать

...
real,allocatable :: tmp (:,:,:)
double precision,allocatable :: tmp2(:,:,:)
...
tmp2=tmp
mean=sum(tmp2)/size(tmp)
...

Решает проблему, но увеличивает время вычислений (не очень, но заметно). Есть ли лучший способ обойти эту проблему? Я не мог найти способ читать одиночные игры из файла прямо в парные. И как этого numpyизбежать?

Спасибо за помощь.

Ваша реализация Fortran страдает двумя основными недостатками:

• Вы смешиваете ввод-вывод и вычисления (и читаете из файла запись за записью).
• Вы не используете векторные / матричные операции.

Эта реализация выполняет ту же операцию, что и ваша, и работает на моей машине в 20 раз быстрее:

program test
  integer gridsize,unit
  real mini,maxi,mean
  real, allocatable :: tmp (:,:,:)

  gridsize=512
  unit=40

  allocate( tmp(gridsize, gridsize, gridsize))

  open(unit=unit,file='T.out',status='old',access='stream',&
       form='unformatted',action='read')
  read(unit=unit) tmp

  close(unit=unit)

  mini = minval(tmp)
  maxi = maxval(tmp)
  mean = sum(tmp)/gridsize**3
  print *, mini, maxi, mean

end program

Идея состоит в том, чтобы за один раз прочитать весь файл в один массив tmp. Тогда я могу использовать функции MAXVAL, MINVALи SUMна массиве непосредственно.

Для проблемы точности: просто используя значения двойной точности и выполняя преобразование на лету как

mean = sum(real(tmp, kind=kind(1.d0)))/real(gridsize**3, kind=kind(1.d0))

лишь незначительно увеличивает время расчета. Я пробовал выполнять операцию поэлементно и по частям, но это только увеличило необходимое время на уровне оптимизации по умолчанию.

При -O3поэлементное сложение на ~ 3% лучше, чем операция с массивом. Разница между операциями с двойной и одинарной точностью на моей машине в среднем составляет менее 2% (отдельные прогоны отклоняются намного больше).

Вот очень быстрая реализация с использованием LAPACK:

program test
  integer gridsize,unit, i, j
  real mini,maxi
  integer  :: t1, t2, rate
  real, allocatable :: tmp (:,:,:)
  real, allocatable :: work(:)
!  double precision :: mean
  real :: mean
  real :: slange

  call system_clock(count_rate=rate)
  call system_clock(t1)
  gridsize=512
  unit=40

  allocate( tmp(gridsize, gridsize, gridsize), work(gridsize))

  open(unit=unit,file='T.out',status='old',access='stream',&
       form='unformatted',action='read')
  read(unit=unit) tmp

  close(unit=unit)

  mini = minval(tmp)
  maxi = maxval(tmp)

!  mean = sum(tmp)/gridsize**3
!  mean = sum(real(tmp, kind=kind(1.d0)))/real(gridsize**3, kind=kind(1.d0))
  mean = 0.d0
  do j=1,gridsize
    do i=1,gridsize
      mean = mean + slange('1', gridsize, 1, tmp(:,i,j),gridsize, work)
    enddo !i
  enddo !j
  mean = mean / gridsize**3

  print *, mini, maxi, mean
  call system_clock(t2)
  print *,real(t2-t1)/real(rate)

end program

Это использует матрицу одинарной точности 1-норму SLANGEдля столбцов матрицы. Время выполнения даже быстрее, чем подход с использованием функций массива с одинарной точностью, и не показывает проблемы с точностью.

Numpy работает быстрее, потому что вы написали гораздо более эффективный код на python (а большая часть numpy backend написана на оптимизированном Fortran и C) и ужасно неэффективный код на Fortran.

Посмотрите на свой код на Python. Вы загружаете сразу весь массив, а затем вызываете функции, которые могут работать с массивом.

Посмотрите на свой код фортрана. Вы читаете одно значение за раз и выполняете с ним некоторую логику ветвления.

Большая часть вашего несоответствия - это фрагментированный ввод-вывод, который вы написали на Фортране.

Вы можете написать Fortran примерно так же, как вы писали питон, и вы обнаружите, что так он работает намного быстрее.

program test
  implicit none
  integer :: gridsize, unit
  real :: mini, maxi, mean
  real, allocatable :: array(:,:,:)

  gridsize=512
  allocate(array(gridsize,gridsize,gridsize))
  unit=40
  open(unit=unit, file='T.out', status='old', access='stream',&
       form='unformatted', action='read')
  read(unit) array    
  maxi = maxval(array)
  mini = minval(array)
  mean = sum(array)/size(array)
  close(unit)
end program test