Jak odrętwienie może być o wiele szybsze niż moja rutyna Fortran?

Otrzymuję tablicę 512 ^ 3 reprezentującą rozkład temperatury z symulacji (napisanej w języku Fortran). Tablica jest przechowywana w pliku binarnym o rozmiarze około 1/2 GB. Muszę znać minimum, maksimum i średnią tej tablicy, a ponieważ i tak wkrótce będę musiał zrozumieć kod Fortran, postanowiłem spróbować i wymyśliłem następującą bardzo prostą procedurę.

  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)

Zajmuje to około 25 sekund na plik na używanym komputerze. Wydało mi się to dość długie, więc poszedłem do przodu i wykonałem następujące czynności w Pythonie:

    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)

Oczywiście spodziewałem się, że to będzie szybsze, ale naprawdę mnie zaskoczyło. W identycznych warunkach zajmuje to mniej niż sekundę. Średnia odbiega od tej znalezionej przez moją rutynę Fortran (którą również uruchomiłem z 128-bitowymi liczbami zmiennoprzecinkowymi, więc w jakiś sposób bardziej jej ufam), ale tylko na 7. cyfrze znaczącej.

Jak numpy może być tak szybki? Chodzi mi o to, że musisz spojrzeć na każdy wpis tablicy, aby znaleźć te wartości, prawda? Czy robię coś bardzo głupiego w mojej rutynie Fortran, żeby zajęło to dużo więcej czasu?

EDYTOWAĆ:

Aby odpowiedzieć na pytania w komentarzach:

• Tak, również uruchomiłem procedurę Fortran z 32-bitowymi i 64-bitowymi pływakami, ale nie miało to wpływu na wydajność.
• Użyłem, iso_fortran_envktóry zapewnia 128-bitowe liczby zmiennoprzecinkowe.
• Używając 32-bitowych liczb zmiennoprzecinkowych, moja średnia jest jednak nieco wyłączona, więc precyzja jest naprawdę problemem.
• Uruchomiłem obie procedury na różnych plikach w różnej kolejności, więc buforowanie powinno być sprawiedliwe w porównaniu, jak sądzę?
• Właściwie próbowałem otworzyć MP, ale czytać z pliku w różnych pozycjach w tym samym czasie. Po przeczytaniu twoich komentarzy i odpowiedzi brzmi to teraz naprawdę głupio i sprawiło, że rutyna również zajęła znacznie więcej czasu. Mógłbym spróbować operacji na macierzy, ale może to nawet nie będzie konieczne.
• Pliki mają rozmiar 1 / 2G, to była literówka, dzięki.
• Spróbuję teraz implementacji tablicy.

EDYCJA 2:

Zaimplementowałem to, co sugerowali @Alexander Vogt i @casey w ich odpowiedziach i jest to równie szybkie, jak, numpyale teraz mam problem z precyzją, jak wskazał @Luaan. Używając 32-bitowej tablicy zmiennoprzecinkowej, średnia obliczona przez sumto 20% zniżki. Robić

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

Rozwiązuje problem, ale wydłuża czas przetwarzania (nie bardzo, ale zauważalnie). Czy jest lepszy sposób na obejście tego problemu? Nie mogłem znaleźć sposobu na odczytanie singli bezpośrednio z pliku do podwójnych. Jak tego numpyuniknąć?

Dziękuję za całą dotychczasową pomoc.

Twoja implementacja Fortran ma dwie poważne wady:

• Możesz mieszać IO i obliczenia (i czytać z pliku wpis po wpisie).
• Nie używasz operacji wektor / macierz.

Ta implementacja wykonuje tę samą operację co twoja i jest 20-krotnie szybsza na moim komputerze:

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

Chodzi o to, aby wczytać cały plik do jednej tablicy tmpza jednym razem. Następnie można użyć funkcji MAXVAL, MINVALi SUMna wprost tablicy.

W kwestii dokładności: wystarczy użyć wartości podwójnej precyzji i wykonać konwersję w locie jako

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

tylko nieznacznie wydłuża czas obliczeń. Próbowałem wykonać operację z podziałem na elementy i w plasterkach, ale to tylko zwiększyło wymagany czas na domyślnym poziomie optymalizacji.

W -O3przypadku dodawanie elementarne działa o ~ 3% lepiej niż operacja na tablicy. Na mojej maszynie różnica między podwójną i pojedynczą precyzją jest mniejsza niż 2% - średnio (poszczególne przebiegi różnią się znacznie bardziej).

Oto bardzo szybka implementacja przy użyciu 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

To używa normy pojedynczej precyzji macierzy 1 SLANGEw kolumnach macierzy. Czas wykonywania jest nawet szybszy niż podejście wykorzystujące funkcje tablicowe o pojedynczej precyzji - i nie wykazuje problemu z precyzją.

Numpy jest szybsze, ponieważ napisałeś znacznie wydajniejszy kod w Pythonie (a większość zaplecza Numpy jest napisana w zoptymalizowanym Fortranie i C) i strasznie nieefektywny kod w Fortranie.

Spójrz na swój kod w Pythonie. Ładujesz całą tablicę naraz, a następnie wywołujesz funkcje, które mogą operować na tablicy.

Spójrz na swój kod Fortran. Odczytujesz jedną wartość na raz i wykonujesz z nią logikę rozgałęziania.

Większość twoich rozbieżności to fragmentaryczne IO, które napisałeś w Fortranie.

Możesz napisać Fortrana w taki sam sposób, jak napisałeś Pythona, a przekonasz się, że działa on znacznie szybciej.

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