Jak numpy może być o wiele szybszy niż mój Fortran?

Otrzymuję tablicę 512^3 przedstawiającą rozkład temperatury z symulacji (zapisanej w Fortran). Tablica jest przechowywana w pliku binarnym o rozmiarze około 1/2G. Muszę znać minimum, maksimum i średnią tej tablicy, a ponieważ wkrótce będę musiał zrozumieć kod Fortran i tak, postanowiłem dać mu szansę i wymyśliłem następującą bardzo łatwą rutynę.

  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)

To zajmuje około 25 sekund na plik na komputerze, którego używam. To wydawało mi się dość długie, więc poszedłem przed i zrobił 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)
Spodziewałam się, że będzie szybciej, ale byłam zdumiona. To trwa mniej niż sekundę w identycznych warunkach. Średnia odbiega od tej, którą znajduje mój program Fortran (który również uruchomiłem z pływakami 128-bitowymi, więc jakoś bardziej mu ufam), ale tylko na siódmej znaczącej cyfrze lub tak. 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łupota w mojej rutynie Fortran, że tak długo to trwa?

EDIT:

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_env który zapewnia 128-bitowe pływaki.
  • Korzystanie z 32-bitowych pływaków jest trochę wyłączone, więc precyzja jest naprawdę problemem.
  • uruchomiłem obie procedury na różnych plikach w różnych zamówienie, więc buforowanie powinno być uczciwe w porównaniu, jak sądzę ?
  • próbowałem otworzyć MP, ale czytać z pliku w różnych pozycjach w tym samym czasie. Po przeczytaniu komentarzy i odpowiedzi brzmi to naprawdę głupio teraz i to sprawiło, że rutyna trwa znacznie dłużej. Mogę spróbować z operacjami array, ale może to nawet nie będzie konieczne.
  • pliki mają rozmiar 1 / 2G, to była literówka, dzięki.
  • spróbuję tablicy wdrożenie teraz.

Edytuj 2:

Zaimplementowałem to, co @Alexander Vogt i @ casey zasugerowali w swoich odpowiedziach, i jest to tak szybkie, jak numpy, ale teraz mam problem z precyzją, jak zauważył @Luaan. Używając 32-bitowej tablicy zmiennoprzecinkowej średnia obliczona przez sum wynosi 20%. Doing 

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

Rozwiązuje problem, ale wydłuża czas obliczeniowy (nie bardzo, ale zauważalnie). Czy istnieje lepszy sposób obejścia tego problemu? I couldn ' t find a way to read singles z pliku bezpośrednio do podwaja. Jak tego uniknąć?

Dzięki za pomoc.
Author: Mysticial, 2015-11-15
Source

2 answers

Twoja implementacja Fortrana ma dwie główne niedociągnięcia:

  • mieszasz IO i obliczenia (i odczytywasz z pliku wpis po wpisie).
  • nie używasz operacji wektorowych / macierzowych.

Ta implementacja wykonuje tę samą operację co twoja i jest szybsza o współczynnik 20 na mojej Maszynie: 

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 tmp za jednym zamachem. Następnie mogę korzystać z funkcji MAXVAL, MINVAL, oraz SUM bezpośrednio na tablicy.

Dla kwestii dokładności: po prostu używając podwójnych wartości dokładności i wykonując konwersję w locie jako 

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

Tylko nieznacznie zwiększa czas obliczeń. Próbowałem wykonać element operacji w plasterkach, ale to tylko wydłużyło wymagany czas na domyślnym poziomie optymalizacji.

Przy -O3 dodawanie elementów wykonuje ~3% lepiej niż operacja tablicy. Różnica między podwójnym a pojedyncze operacje precyzyjne to mniej niż 2% na mojej maszynie - średnio (poszczególne przebiegi różnią się znacznie bardziej).

Oto bardzo szybka implementacja przy użyciu LAPACKA: 

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 wykorzystuje macierz jednopunktową 1-normę SLANGE na kolumnach macierzy. Czas wykonania jest jeszcze szybszy niż podejście przy użyciu funkcji single precision array - i nie pokazuje problemu precyzji.

 109
Author: Alexander Vogt,
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/doraprojects.net/template/agent.layouts/content.php on line 54
2015-11-16 18:25:39

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

Spójrz na swój kod Pythona. Ładujesz całą tablicę na raz, a następnie wywołujesz funkcje, które mogą działać na tablicy.

Spójrz na swój kod fortran. Odczytujesz jedną wartość na raz i robisz z nią logikę rozgałęzień.

Większość twojej rozbieżności to fragmentaryczne IO, które masz napisane w Fortran.

Możesz napisać Fortran dokładnie tak samo jak napisałeś Pythona i przekonasz się, że działa 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
 54
Author: casey,
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/doraprojects.net/template/agent.layouts/content.php on line 54
2015-11-16 02:09:15