Comment peut-numpy être beaucoup plus rapide que celle de ma routine Fortran?

Question

Je reçois un 512^3 tableau représentant une distribution de la Température à partir d'une simulation (écrit en Fortran). Le tableau est stocké dans un fichier binaire, c'est environ 1/2G dans la taille. J'ai besoin de savoir le minimum, le maximum et la moyenne de ce tableau et que je vais bientôt besoin de comprendre le code Fortran de toute façon, j'ai décidé de lui donner un aller et venu avec la suite très facile de routine.

  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)

Cela prend environ 25 secondes par fichier sur la machine que j'utilise. Ce qui m'a frappé comme étant plutôt long et donc je suis allé de l'avant et ne le suit en 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)

Maintenant, je m'attendais à cette à être plus rapide, bien sûr, mais j'ai été vraiment époustouflé. Il faut moins d'une seconde dans des conditions identiques. La moyenne s'écarte de celui de ma routine Fortran trouve (j'ai aussi couru 128 bits, les flotteurs, donc j'ai un peu de confiance plus), mais seulement à la 7ème chiffre significatif.

Comment peut-numpy être si rapide? Je veux dire, vous avez à regarder à chaque entrée d'une table pour trouver ces valeurs, droit? Suis-je en train de faire quelque chose de très stupide dans ma routine Fortran pour qu'il prenne beaucoup plus de temps?

EDIT:

Pour répondre aux questions dans les commentaires:

• Oui, aussi, j'ai couru la routine Fortran avec 32-bit et 64-bit flotte mais il n'a eu aucun impact sur les performances.
• J'ai utilisé iso_fortran_env qui fournit un cryptage 128 bits flotteurs.
• En utilisant 32 bits flotteurs ma moyenne est éteint tout à fait un peu, donc la précision est vraiment un problème.
• J'ai couru les deux routines sur les différents fichiers dans un ordre différent, de sorte que la mise en cache doit avoir été juste dans la comparaison, je suppose ?
• J'ai effectivement essayé de l'ouvrir MP, mais à le lire à partir du fichier à différentes positions en même temps. Après avoir lu vos commentaires et de vos réponses cela semble vraiment stupide maintenant et il fait de la routine de prendre beaucoup plus de temps. Je pourrais faire un essai sur le tableau des opérations, mais peut-être de ne pas même être nécessaire.
• Les fichiers sont en fait 1/2G en taille, c'était une faute de frappe, Merci.
• Je vais essayer de la matrice de mise en œuvre maintenant.

EDIT 2:

J'ai mis ce que @Alexandre Vogt et @casey suggéré dans leurs réponses, et il est aussi rapide que numpy , mais maintenant j'ai un problème de précision de @Luaan souligné que je pourrais obtenir. À l'aide d'une virgule flottante de 32 bits de la matrice de la moyenne calculée en sum 20% off. Faire

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

Résout le problème, mais augmente le temps de calcul (ce n'est pas beaucoup, mais sensiblement). Est-il une meilleure façon de contourner ce problème? Je ne pouvais pas trouver un moyen de lire les singles de directement le fichier en double. Et comment est - numpy éviter cela?

Merci pour toute l'aide jusqu'à présent.

Answer 1

Votre Fortran mise en œuvre souffre de deux défauts majeurs:

• Vous mélangez les IO et les calculs (et lire à partir du fichier d'entrée par entrée).
• Vous n'avez pas l'utilisation de vecteurs/matrices opérations.

Cette mise en œuvre n'effectuer la même opération que le vôtre et il est plus rapide d'un facteur 20 sur ma machine:

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

L'idée est de lire dans l'ensemble du fichier dans un tableau tmp d'un seul coup. Alors, je peux utiliser les fonctions MAXVAL, MINVAL, et SUM sur le tableau directement.

---

Pour la question de l'exactitude: il vous suffit d'utiliser le double des valeurs de précision et de faire la conversion à la volée sur le

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

que de façon marginale augmente le temps de calcul. J'ai essayé d'effectuer l'opération de l'élément de sage et en tranches, mais qui ne fait qu'accroître le temps nécessaire à l'optimisation par défaut.

En -O3, l'élément-sage addition effectue ~3 % de mieux que la matrice de l'opération. La différence entre le double et le simple précision des opérations est inférieur à 2% sur ma machine - en moyenne (la personne s'exécute s'en écarter beaucoup plus).

---

Voici une très rapide mise en œuvre à l'aide de 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

Il utilise la simple précision de la matrice 1-norme SLANGE sur les colonnes de la matrice. Le temps d'exécution est plus rapide que l'approche à l'aide de simple précision les fonctions de tableau - et ne montre pas la question de précision.

Answer 2

Le numpy est plus rapide parce que vous avez écrit beaucoup plus efficace de code en python (et une grande partie de la numpy backend est écrit dans optimisé Fortran et C) et terriblement inefficace code en Fortran.

Regardez votre code python. Vous chargez le tableau d'ensemble à la fois et ensuite appeler des fonctions qui peuvent fonctionner sur un tableau.

Regardez votre code fortran. Vous lire une valeur à la fois et faire un peu de ramification logique avec elle.

La majorité de vos écart est le manque de IO que vous avez écrit en Fortran.

Vous pouvez écrire le Fortran peu près de la même manière que vous avez écrit en python et vous verrez qu'il fonctionne beaucoup plus rapidement de cette façon.

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