Comment Numpy peut-il être plus rapide que ma routine Fortran?

Je reçois un tableau 512 ^ 3 représentant une dissortingbution de température provenant d’une simulation (écrite en Fortran). Le tableau est stocké dans un fichier binary d’environ 1 / 2G. J’ai besoin de connaître le minimum, le maximum et la moyenne de ce tableau et comme je vais bientôt avoir besoin de comprendre le code Fortran, j’ai décidé de l’essayer et j’ai trouvé la routine très facile suivante.

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. Cela m’a paru assez long et je suis donc allé de l’avant et j’ai fait ce qui 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 à ce que ce soit 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 celle que ma routine Fortran trouve (que j’ai aussi utilisée avec des flotteurs de 128 bits, donc je lui fais confiance en quelque sorte), mais seulement au 7ème chiffre significatif.

Comment numpy peut-il être si rapide? Je veux dire que vous devez regarder chaque entrée d’un tableau pour trouver ces valeurs, non? Est-ce que je fais quelque chose de très stupide dans ma routine de Fortran pour que cela prenne beaucoup plus de temps?

MODIFIER:

Pour répondre aux questions dans les commentaires:

  • Oui, j’ai également exécuté la routine Fortran avec des flotteurs 32 bits et 64 bits, mais cela n’a eu aucun impact sur les performances.
  • J’ai utilisé iso_fortran_env qui fournit des flottants de 128 bits.
  • En utilisant les flottants 32 bits, ma moyenne est un peu en retrait, donc la précision est vraiment un problème.
  • J’ai couru les deux routines sur différents fichiers dans un ordre différent, donc la mise en cache aurait dû être juste dans la comparaison, je suppose?
  • En fait, j’ai essayé d’ouvrir MP, mais de lire le fichier à différentes positions en même temps. Après avoir lu vos commentaires et réponses, cela semble vraiment stupide maintenant et cela a rendu la routine beaucoup plus longue aussi. Je vais peut-être essayer les opérations sur les tableaux, mais ce ne sera peut-être même pas nécessaire.
  • Les fichiers ont en fait une taille de 1 / 2G, c’était une faute de frappe, merci.
  • Je vais essayer la mise en œuvre du tableau maintenant.

EDIT 2:

J’ai implémenté ce que @Alexander Vogt et @casey ont suggéré dans leurs réponses, et c’est aussi rapide que numpy mais maintenant j’ai un problème de précision comme @Luaan l’a fait remarquer. En utilisant un tableau flottant de 32 bits, la moyenne calculée par sum est de 20%. 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 (pas beaucoup, mais sensiblement). Y a-t-il une meilleure façon de contourner ce problème? Je ne pouvais pas trouver un moyen de lire les singles du fichier directement en double. Et comment numpy évite-t-il cela?

Merci pour toute l’aide jusqu’ici.

    Votre implantation Fortran souffre de deux lacunes majeures:

    • Vous mélangez les entrées-sorties et les calculs (et lisez les entrées de fichier par entrée).
    • Vous n’utilisez pas les opérations vectorielles / masortingcielles.

    Cette implémentation effectue la même opération que la vôtre et 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 le fichier entier dans un tableau tmp en une seule fois. Ensuite, je peux utiliser directement les fonctions MAXVAL , MINVAL et SUM sur le tableau.


    Pour le problème de la précision: utiliser simplement des valeurs de double précision et effectuer la conversion à la volée

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

    n’augmente que marginalement le temps de calcul. J’ai essayé d’effectuer l’opération par éléments et par tranches, mais cela n’a fait qu’accroître le temps requirejs au niveau d’optimisation par défaut.

    Pour -O3 , l’addition d’éléments est environ 3% meilleure que celle du tableau. La différence entre les opérations de précision double et simple est inférieure à 2% sur ma machine – en moyenne (les courses individuelles s’écartent de beaucoup plus).


    Voici une implémentation très rapide utilisant 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 

    Cela utilise la masortingce de précision simple SLANGE 1 norme sur les colonnes de la masortingce. L’exécution est encore plus rapide que l’approche utilisant des fonctions de tableau de précision unique – et ne montre pas le problème de précision.

    Le numpy est plus rapide car vous avez écrit du code beaucoup plus efficace en python (et une grande partie du backpy numpy est écrite dans Fortran et C optimisés) et du code terriblement inefficace dans Fortran.

    Regardez votre code python. Vous chargez l’intégralité du tableau à la fois, puis vous appelez des fonctions pouvant fonctionner sur un tableau.

    Regardez votre code fortran. Vous lisez une valeur à la fois et vous faites une logique de twigment.

    La majorité de vos divergences sont les entrées-sorties fragmentées que vous avez écrites dans Fortran.

    Vous pouvez écrire le Fortran à peu près de la même manière que vous avez écrit le python et vous constaterez que son fonctionnement est beaucoup plus rapide.

     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