Wie kann Numpy so viel schneller sein als meine Fortran-Routine?

Ich erhalte ein 512 ^ 3-Array, das eine Temperaturverteilung aus einer Simulation darstellt (geschrieben in Fortran). Das Array wird in einer Binärdatei mit einer Größe von etwa 1 / 2G gespeichert. Ich muss das Minimum, Maximum und den Mittelwert dieses Arrays kennen und da ich den Fortran-Code sowieso bald verstehen muss, habe ich beschlossen, es auszuprobieren, und mir die folgende sehr einfache Routine ausgedacht.

  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)

Dies dauert ungefähr 25 Sekunden pro Datei auf dem von mir verwendeten Computer. Das kam mir ziemlich lang vor und so machte ich in Python Folgendes:

    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)

Nun, ich hatte erwartet, dass dies natürlich schneller sein würde, aber ich war wirklich überwältigt. Unter identischen Bedingungen dauert es weniger als eine Sekunde. Der Mittelwert weicht von dem ab, den meine Fortran-Routine findet (den ich auch mit 128-Bit-Floats ausgeführt habe, also vertraue ich ihm irgendwie mehr), aber nur auf der 7. signifikanten Ziffer oder so.

Wie kann Numpy so schnell sein? Ich meine, Sie müssen sich jeden Eintrag eines Arrays ansehen, um diese Werte zu finden, oder? Mache ich in meiner Fortran-Routine etwas sehr Dummes, damit es so viel länger dauert?

BEARBEITEN:

So beantworten Sie die Fragen in den Kommentaren:

• Ja, ich habe auch die Fortran-Routine mit 32-Bit- und 64-Bit-Floats ausgeführt, aber sie hatte keinen Einfluss auf die Leistung.
• Ich habe verwendet, iso_fortran_envdie 128-Bit-Floats bietet.
• Bei Verwendung von 32-Bit-Floats ist mein Mittelwert jedoch ziemlich niedrig, sodass Präzision wirklich ein Problem darstellt.
• Ich habe beide Routinen auf verschiedenen Dateien in unterschiedlicher Reihenfolge ausgeführt, also sollte das Caching im Vergleich fair sein, denke ich?
• Ich habe tatsächlich versucht, MP zu öffnen, aber gleichzeitig an verschiedenen Positionen aus der Datei zu lesen. Nachdem Sie Ihre Kommentare und Antworten gelesen haben, klingt dies jetzt wirklich dumm und die Routine hat auch viel länger gedauert. Ich könnte es mit den Array-Operationen versuchen, aber vielleicht ist das gar nicht nötig.
• Die Dateien sind tatsächlich 1 / 2G groß, das war ein Tippfehler, danke.
• Ich werde jetzt die Array-Implementierung versuchen.

EDIT 2:

Ich habe implementiert, was @Alexander Vogt und @casey in ihren Antworten vorgeschlagen haben, und es ist so schnell wie, numpyaber jetzt habe ich ein Präzisionsproblem, wie @Luaan darauf hingewiesen hat, dass ich es bekommen könnte. Bei Verwendung eines 32-Bit-Float-Arrays sumbeträgt der berechnete Mittelwert 20%. Tun

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

Behebt das Problem, erhöht aber die Rechenzeit (nicht sehr stark, aber spürbar). Gibt es einen besseren Weg, um dieses Problem zu umgehen? Ich konnte keinen Weg finden, Singles aus der Datei direkt in Doppel zu lesen. Und wie numpyvermeidet man das?

Vielen Dank für die bisherige Hilfe.

Ihre Fortran-Implementierung weist zwei Hauptmängel auf:

• Sie mischen E / A und Berechnungen (und lesen Eintrag für Eintrag aus der Datei).
• Sie verwenden keine Vektor- / Matrixoperationen.

Diese Implementierung führt den gleichen Vorgang wie Ihre aus und ist auf meinem Computer um den Faktor 20 schneller:

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

Die Idee ist, die gesamte Datei auf einmal in ein Array einzulesen tmp. Dann kann ich die Funktionen nutzen MAXVAL, MINVALund SUMauf dem Array direkt.

Für das Genauigkeitsproblem: Verwenden Sie einfach Werte mit doppelter Genauigkeit und führen Sie die Konvertierung im laufenden Betrieb durch

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

erhöht die Rechenzeit nur unwesentlich. Ich habe versucht, die Operation elementweise und in Slices auszuführen, aber das hat nur die erforderliche Zeit auf der Standardoptimierungsstufe erhöht.

Bei -O3ist die elementweise Addition ~ 3% besser als die Array-Operation. Der Unterschied zwischen Operationen mit doppelter und einfacher Genauigkeit beträgt auf meiner Maschine weniger als 2% - im Durchschnitt (die einzelnen Läufe weichen weitaus stärker ab).

Hier ist eine sehr schnelle Implementierung mit 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

Dies verwendet die Matrix 1-Norm SLANGEmit einfacher Genauigkeit für Matrixspalten. Die Laufzeit ist sogar schneller als der Ansatz mit Array-Funktionen mit einfacher Genauigkeit - und zeigt das Problem der Genauigkeit nicht.

Das Numpy ist schneller, weil Sie viel effizienteren Code in Python geschrieben haben (und ein Großteil des Numpy-Backends in optimiertem Fortran und C geschrieben ist) und schrecklich ineffizienten Code in Fortran.

Schauen Sie sich Ihren Python-Code an. Sie laden das gesamte Array auf einmal und rufen dann Funktionen auf, die auf einem Array ausgeführt werden können.

Sehen Sie sich Ihren fortran-Code an. Sie lesen jeweils einen Wert und führen damit eine Verzweigungslogik durch.

Der größte Teil Ihrer Diskrepanz ist das fragmentierte IO, das Sie in Fortran geschrieben haben.

Sie können den Fortran genauso schreiben, wie Sie den Python geschrieben haben, und Sie werden feststellen, dass er auf diese Weise viel schneller läuft.

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